5 Метод конечных разностей

Предположим, что решается одномерная линейная краевая задача для уравнения 2-го порядка, т.е. требуется определить функцию y(x), удовлетворяющую заданному дифференциальному уравнению

(5.1)

на отрезке axb вместе с краевыми условиями при x=a и x=b

(5.2)

( ). Для решения этой задачи методом конечных разностей (МКР) прежде всего производится дискретизация независимой переменной x, т.е. строится множество (или сетка) n+1 точек (узлов) x_k (k=0,1,2,...,n) на отрезке axb, причем x₀=a, x_n=b. Шаг сетки h_k  x_k+1–x_k из-за удобства программирования чаще берется постоянным (h₁=h₂=…=h_nh), но может быть и переменным (особенно когда заранее известен характер поведения решения).

Следующий этап состоит в замене в дифференциальном уравнении членов, содержащих дифференцирование, разностными отношениями. На основе формулы Тейлора можно получить выражения для производных функции

, ,

где h=x_k+1–x_k=x_k–x_k–1, y_ky(x_k). Подставляя эти выражения в уравнение (5.1) для всех (может быть, за исключением граничных) узлов {x_k}, а также в уравнения граничных условий, получим систему линейных алгебраических уравнений относительно узловых значений {y_k}.

Для оценки точности получаемого решения необходимо провести расчеты для разных значений шага (не менее 3) и убедиться в том, что полученные значения функции в одних и тех же узлах близки между собой и разность их уменьшается, что говорит о стремлении решения к некоторому пределу при h ^ 0.

Следует отметить, что при определенных наборах коэффициентов уравнения (5.1) p(x), q(x), f(x) и граничных условий _i, _i, _i решение рассматриваемой краевой задачи может не существовать или быть не единственным. Кроме того, даже если оно и существует и единственно, то МКР может быть неустойчивым. Не касаясь деталей этого вопроса, ниже мы будем предполагать безусловный характер устойчивости как задачи, так и метода.

МКР может быть с успехом распространен на нелинейные дифференциальные уравнения, на системы уравнений, уравнения в частных производных, а также адаптирован к решению задач Коши для различных дифференциальных уравнений. В любом случае решение дифференциального уравнения конечно-разностным методом сводится к решению системы алгебраических уравнений относительно значений искомой функции на заданном множестве точек. Эта система является линейной или нелинейной в зависимости от того, линейно или нелинейно дифференциальное уравнение. Следует обратить внимание, что МКР непосредственно дает значения функции только на дискретном множестве узловых точек, а решение во всех остальных точках области задачи можно получить с помощью какого-нибудь интерполяционного метода. Из-за этого в точность решения будет вноситься дополнительная погрешность аппроксимации.

Предложенный ниже пример позволяет проследить процесс построения разностной схемы и выявить ключевые моменты реализации МКР.

Пример. Методом конечных разностей решить краевую задачу

; .

Решение. В соответствии с МКР в исходном уравнении

вместо первой и второй производной записываются аппроксимации:

.

Разобьем отрезок [1,2] на четыре равные части. Получим 5 узлов x_i=1+hi; i=0, 1, ..., 4; h=0,25. Запишем последовательно уравнения для всех узлов, начиная с x₁ (для узла x₀=1, которому соответствует индекс i=0, не пишем, поскольку для него значение известно из граничного условия: j₀ =j(1)=0).

Для i=1 (x=1,25) – первое уравнение:

;

для i=2 (x=1,5) – второе уравнение:

;

для i=3 (x=1,75) – третье уравнение:

;

для i=4 (x=2) – четвертое уравнение

.

Здесь j₅ – значение неизвестной функции в фиктивном узле x₅=2,25.

Имеем четыре уравнения при пяти неизвестных (j₁, j₂, j₃, j₄, j₅). Недостающее пятое уравнение получим, если аппроксимируем второе граничное условие, соответствующее узлу с номером i=4:

. (пятое уравнение)

Итак, решение полученной системы даст приближенные значения искомой функции в точках x=1,25; x=1,5; x=1,75; x=2. Найти значения в других точках отрезка [1, 2] можно с помощью методов интерполяции.

Порядок решения в системе Maple

> restart;

1. Введем переменную deq для хранения дифференциального уравнения:

> deq:=x^2*diff(y(x),x$2)+2*x*diff(y(x),x)-2*y(x)=-1-2/x;

(y(x)здесь обозначает (x)).

2. Определим конечно-разностные уравнения.

Первое уравнение (переменная eq1) получим, осуществляя подстановку в deq вместо (x) и (x) аппроксимации и соответственно, вместо y(x) узловую переменную y1 (т.е. j₁) и, наконец, вместо x значение 1,25:

> eq1:=subs(diff(y(x),x$2)=(y2-2*y1+y0)/0.25^2,diff(y(x),x)=(y2-y0)/(2*0.25),y(x)=y1,x=1.25,deq);

Аналогично зададим следующие три уравнения для точек x=1,5; x=1,75; x=2 – переменные eq2, eq3 и eq4 соответственно:

> eq2:=subs(diff(y(x),x$2)=(y3-2*y2+y1)/0.25^2, diff(y(x),x)=(y3-y1)/(2*0.25),y(x)=y2,x=1.5,deq);

> eq3:=subs(diff(y(x),x$2)=(y4-2*y3+y2)/0.25^2, diff(y(x),x)=(y4-y2)/(2*0.25),y(x)=y3,x=1.75,deq);

> eq4:=subs(diff(y(x),x$2)=(y5-2*y4+y3)/0.25^2, diff(y(x),x)=(y5-y3)/(2*0.25),y(x)=y4,x=2.0,deq);

Пятое уравнение – это аппроксимация граничного условия (2)=1:

> eq5:=(y5-y3)/(2*0.25)=1;

3. Решим эту систему с помощью функции solve. При этом следует положить ₀=0 (или y0⁼0), что следует из первого граничного условия:

> rez:=solve({y0=0,eq1,eq2,eq3,eq4,eq5},{y0,y1,y2, y3,y4,y5});

Получим результат:

(y0, y1, ... здесь обозначают j₀, j₁, ...). Это и есть конечно-разностное решение, которое непосредственно дает значения искомой функции лишь в отдельных точках – узлах сетки {x_k}.

4. Сравним с точным решением:

> st:=dsolve({deq,y(1)=0,D(y)(2)=1},y(x)):

> u:=unapply(subs(st,y(x)),x);

Выведем найденные числа {y_i} и значения функции u(x) в тех узлах, для которых найдено конечно-разностное решение:

> X:=[1,1.25,1.5,1.75,2]: Y:=subs(rez,[y0,y1,y2,y3,y4]);

> for j from 1 to 5 do printf(`x=%4.2g Y=%7.5g u=%7.5g\n`, X[j],Y[j],evalf(u(X[j]))); od;

x=1.00 Y=0.00000 u=0.00000

x=1.25 Y= .75290 u= .76700

x=1.50 Y=1.21836 u=1.23888

x=1.75 Y=1.55329 u=1.57806

x=2.00 Y=1.82171 u=1.85000

Совпадение есть. Погрешность ~1%.

5. Выполним построения графика точного решения и полученных приближенных узловых значений

> i:='i';

> g3:=plot([[X[i],Y[i]]$i=1..5],x=1..2,style=POINT, symbol=CIRCLE,color=black);

> plots[display](g3);

(Графическое построение не приводится).

6. Найдем приближенное решение в виде непрерывной функции с помощью интерполяционного метода

Сначала проведем глобальную интерполяцию

> Digits:=6: f:=interp(X,Y,a);

> fu:=proc(x) global a; a:=x; RETURN(f) end;

Построим графики интерполяционной функции fu(x) и точного решения вместе с узлами интерполяции

> g1:=plot(fu(x),x=1..2,color=blue):g2:=plot(u(x),

x=1..2,color=red):

> plots[display](g1,g2,g3);

Рис. 5.1.

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

> readlib(spline);

> f1:=spline(X,Y,x,linear);

> g4:=plot(f1,x=1..2,color=black,thickness=2):

> g2:=plot(u(x),x=1..2,color=red):

> g5:=plot([[X[i],Y[i]]$i=1..5],x=1..2.1, style=POINT,symbol=CIRCLE):

> plots[display]y(g2,g4,g5);

Рис. 5.2.

7. Проверка сходимости метода конечных разностей.

Сходимость означает, что с уменьшением расстояния h между ближайшими узлами (а это приводит к увеличению числа узлов) точность МКР должна возрастать, т.е. должна стремиться к нулю разность |y_i–y(x_i)| для любого узла с номером i при h^0. Таким образом, для проверки сходимости нужно провести серию расчетов по рассмотренной схеме МКР для разного числа узлов n.

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

или

Следует, однако, заметить, что тот способ реализации МКР в системе Maple, что приведен выше, достаточно нагляден, но с ростом числа узлов становится крайне громоздким. Действительно, если это число равно n, то придется вручную вводить столько же уравнений. Попробуйте оценить, на сколько возрастет текст программы при n =10, n =20, n =100.

Поэтому приведем другую, более универсальную, Maple-программу, дающую решение данной краевой задачи методом конечных разностей при любом n (т.е. n выступает как параметр):

> restart;

> Digits:=15;

> a:=1; b:=2; n:=10; h:=(b-a)/n;

Граничные условия зададим в общем виде (10), устанавливая для коэффициентов соответствующие значения (A0 сооветствует ₀, B0 сооветствует ₀, C0 сооветствует ₀ и т.д.):

> A0:=0; B0:=4; C0:=0; A1:=1.; B1:=0; C1:=1;

> diffeqn:=x^2*diff(y(x),x$2)+2*x*diff(y(x),x)-

2*y(x)=-1-2/x;

> Eqns_bas:={seq(subs(diff(y(x),x$2)=

(g[k+1]-2*g[k]+g[k-1])/h^2,diff(y(x),x)=(g[k+1]-g[k-1])/(2.*h),y(x)=g[k],

x=a+h*k,diffeqn),k=0..n)};

> eq1:={A0*(g[1]-g[-1])/2./h+B0*g[0]=C0}; # граничное условие в точке x=a

> eq2:={A1*(g[n+1]-g[n-1])/2./h+B1*g[n]=C1}; # граничное условие в точке x=b

> rez:=solve(Eqns_bas union eq1 union eq2, {g[i]$i=-1..n+1}):

> Y:=evalf([seq(subs(rez,g[k]),k=0..n)]):

> X:=evalf[5]([seq(a+k*h,k=0..n)]):

> gr1:=plot([[X[i],Y[i]]$i=1..n+1],x=a..b,

style=POINT,

symbol=CIRCLE,color=black):

> readlib(spline); # – получение непрерывного решения с помощью сплайн-интерполяции

> Digits:=10: f1:=spline(X,Y,x,linear):

> gr2:=plot(f1,x=a..b,color=blue):

> x:='x': y:='y':

> p:=dsolve({diffeqn,A0*D(y)(a)+B0*y(a)=C0,

A1*D(y)(b)+B1*y(b)=C1},y(x));

> u:=unapply(subs(p,y(x)),x); # u(x) – точное решение

> gr3:=plot(u(x),x=a..b,color=red):

> plots[display]([gr1,gr2,gr3])

Вычисление величин ₁ и ₂ – норм ошибки решения:

> delta1:=sum('( Y[i]-u(X[i]) )^2',

'i'=1..n+1)^(0.5);

> delta2:=max( seq(abs(Y[i]-u(X[i])),i=1..n+1) );

Следующие строки используются для табулирования:

> np:=trunc(0.1/h); if(np=0) then np:=1; fi;

> for j from 1 by np to n+1 do printf(`x=%4.2g Y=%8.5g u=%8.5g\n`,

X[j],Y[j],evalf(u(X[j]))); od;

Результаты вычисления величин ₁ и ₂ в зависимости от n приведены в таблице.

	n=4	n=8	n=16	n=50	n=100
1	0,04509	0,01527	0,005263	9,344910–4	3,288310–4
2	0,02829	0,007121	0,001784	1,182810–4	4,569910–5

Видно, что сходимостью данный метод конечных разностей обладает.