Решение. Будем решать задачу в среде 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. Уравнение Эйлера для двойного интеграла.
Решить вариационную задачу для функционала, если на контуре области функция принимает заданные значения:
.
Решение. Будем решать задачу в среде Maple. Определим подынтегральную функцию. Для удобства обозначим производные функции 
по 
и по 
соответственно, через 
и 
:
> F:=proc(x, y,u, ux, uy) options operator, arrow;
ux^2+uy^2+2*y*u*cos(x)
end proc;
Уравнение Эйлера в этих обозначениях имеет вид:
.
Вычисляем последовательно производные в этом уравнении
> a1:=diff(F(x, y,u, ux, uy),u);
a2:=diff(F(x, y,u, ux, uy),ux);
a3:=diff(F(x, y,u, ux, uy),uy);
Составляем уравнение Эйлера для функционала
> EulerEq:=a1-(diff(subs(ux=diff(u(x, y),x),a2),x))-
(diff(subs(uy=diff(u(x, y),y),a3),y))=0;
Таким образом, получили следующее уравнение
> pde:=-(1/2)*op(2,lhs(EulerEq))-
(1/2)*op(3,lhs(EulerEq))=(1/2)*op(1,lhs(EulerEq));
Сформируем теперь граничные условия. Для этого определим граничную функцию
> f:=proc(x, y) options operator, arrow;
(1/10)*x+(1/50)*y^2
end proc;
|
Из за большого объема этот материал размещен на нескольких страницах:
1 2 3 4 5 6 7 8 9 10 |
