Пример 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);