Методические указания и варианты курсовой работы по дисциплине «Вариационные методы в математической физике» (стр. 8 )

Отобразим решение на графиках

> plots[display]({pUp_0,pUp_1,pUp_2,py});

Видим, что и в этом случае удержание трех членов ряда вполне достаточно. Покажем полученные аппроксимации

> Up(x,0);Up(x,1);Up(x,2);

Пример 2. Конечно-разностный метод Эйлера.

Найти приближенное решение задачи о минимуме функционала

.

Показать сходимость полученного приближенного решения к точному при возрастании числа элементов, используемых в аппроксимации. Построить графики решений.

Решение. Будем решать задачу в среде Maple. Аппроксимация подынтегральной функции конечными разностями:

> restart; interface(displayprecision = 5):

> F:=proc(Y, m,h)

(Y[m+1]-Y[m])^2/h^2+Y[m]^2+2*X[m]*Y[m]

end proc;

Интеграл заменяем суммой по формуле прямоугольников

> S1:=proc(h, F,N) options operator, arrow;

h*(sum(F(Y, i,h),i = 0 .. N-1))

end proc;

Задаем пределы интегрирования:

> a := 0; b := 1;

Выбираем число узловых точек и определяем шаг интегрирования:

> N:=10: h:=(b-a)/N;

Вычисляем абсциссы вершин ломаной

> for j from 0 to N do X[j] := a+j*h end do:

Функционал как функция ординат вершин ломаной:

> Phi:=S1(h, F,N);

Учет граничных условий:

> Y[0]:=0; Y[N]:=0; Phi;

Составляем минимизирующую систему уравнений:

> for k to N-1 do

eq[k]:=evalf(diff(Phi, Y[k]))=0

end do:

var := {}: eqns := {}:

for k to N-1 do

var:=`union`(var,{Y[k]}):

eqns := `union`(eqns, {eq[k]}):

end do:

nops(var); nops(eqns);

Решаем систему:

> res:=solve(eqns, var);assign(res):

Сформируем список точек вершин ломаной:

> for j from 0 to N do P[j]:=[X[j],Y[j]] end do:

L:=[seq(P[k-1],k = 1 .. N+1)]:

Построим график решения:

> plot(L, x=0..1,title=cat("Число узлов N = ",

convert(N, string)),titlefont=[roman,15],

labelfont[Helvetica, roman,14],

legend=["метод Эйлера"],gridlines=true);

Для сравнения найдем точное решение задачи.

> with(VariationalCalculus):

> f:=(diff(y(x),x))^2+y(x)^2+2*x*y(x);

> ode:=EulerLagrange(f, x,y(x));

> problem := `union`(ode, {y(0) = 0, y(1) = 0});

> dsolve(problem, y(x));

> simplify(convert(%, trig));

> y := unapply(rhs(%), x);

Построим графики приближенного и точного решений:

> plot([L, y(x)],x=0..1,

title=cat("Число узлов N = ",

convert(N, string)),titlefont=[roman,15],

style=[point, line],labelfont[Helvetica, roman,14],

legend=["метод Эйлера","Точное решение"],

gridlines=true);

Пример 3. Метод Ритца для двойного интеграла.

Используя метод Ритца, найти экстремали функционала

Исследовать сходимость. Построить графики.

Решение. Будем решать задачу в среде Maple. Выбираем базисные функции и определяем аппроксимирующую функцию:

> restart;

> phi0:=(x, y)->x/10+y^2/50;

> phi:=(x, y,i, j)->

sin(i*Pi*(x-x1)/(x2-x1))*sin(j*Pi*(y-y1)/(y2-y1));

> U:=proc(x, y,M, N)option operator, arrow;

local n, m;

phi0(x, y)+

sum(sum(a[m, n]*phi(x, y,m, n),'m'=1..M),'n'=1..N);

end proc;

Здесь также удобно разработать процедуру для автоматического составления и решения уравнений метода Ритца. Например, такую

> Ritz:=proc(F, M,N, a)local Fu, eqns, var, eq, i,j, res;

Fu:=simplify(int(int(F, x=0..1),y=0..2));

eqns:={}:var:={}:

for i from 1 to M do

for j from 1 to N do

var:=var union {a[i, j]}:

eq[i, j]:=diff(Fu, a[i, j])=0:

eqns:=eqns union {eq[i, j]}:

od:

od:

res:=solve(eqns, var);

assign(res);

end proc:

Вводим теперь граничные точки

> x1:=0;x2:=1;y1:=0;y2:=2;

Определяем подынтегральную функцию

> F:=diff(u(x,y),x)^2-2*diff(u(x,y),y)^2+

2*y*u(x, y)*(sin(Pi*x)+x/5);

Задаемся числом аппроксимирующих функций и решаем задачу

> M:=2:N:=2:a:=array(1..M,1..N):

> F2:=subs(u(x, y)=U(x, y,M, N),F):Ritz(F2,M, N,a):

Посмотрим результат на графике

> p||M:=

> plot(U(0.5,y, M,N),y=0..2,color=black, linestyle=4,

legend=cat(`M=N=`,convert(M, string))):

plots[display]({p2});

Увеличим на единицу (по каждой переменной) число слагаемых в аппроксимации и повторим расчет

> M:=3:N:=3:a:=array(1..M,1..N):

> F3:=subs(u(x, y)=U(x, y,M, N),F):Ritz(F3,M, N,a):

> p||M:=

plot(U(0.5,y, M,N),y=0..2,color=black, linestyle=1,

legend=cat(`M=N=`,convert(> M, string))):

> plots[display]({p2,p3});

> M:=4;N:=4;a:=array(1..M,1..N):

> F4:=subs(u(x, y)=U(x, y,M, N),F):Ritz(F4,M, N,a);

> p||M:=

plot(U(0.5,y, M,N),y=0..2,color=black, linestyle=3,

legend=cat(`M=N=`,convert(M, string))):

> plots[display]({p2,p3,p4});

Видим, что приближенные решения хорошо сходятся, и можно ограничиться значениями M = N = 4.

Пример 4. Уравнение Эйлера для двойного интеграла.