Pers.narod.ru. Алгоритмы. Краевая задача для дифференциального уравнения 2 порядка |
Вот здесь уже есть немало примеров на реализацию численных методов в Excel.
Однако, решать в Excel краевую задачу для дифференциального уравнения второго порядка, по-моему, извращение. Гораздо проще написать программку на Паскале, чем мы сейчас и займёмся. Но сначала немного теории, только самые нужные формулы.
Для дифференциального уравнения второго порядка
поставлена краевая задача:
Здесь u(t) -искомое
решение, A, B - краевые условия,заданные функции коэффициентов. Введем на
отрезке [a,b] разностную
сетку (t0, t1, t2, ..., tM,),
ti=a+τ*i, i=0,1,...,M,
τ=(b-a)/M, M –
параметр задачи (число шагов сетки).
Заменим уравнение разностной схемой
Собирая коэффициенты при yi, yi+1 , yi-1, получим следующую СЛАУ на вектор неизвестных (y0, y1, y2, ..., yM):
(k1—k2/τ)y0 + k2 y1/τ =A
AAi yi-1 - CCiyi + BBiyi+1 = FFi, i=1,2,...,M-1
(l1+l2/τ) yM - l2 yM-1 /τ=B,
где CCi = 2 - g(ti) τ2 , AAi = 1 - p(ti) τ/2, BBi = 1 + p(ti) τ/2, FFi = τ2f(ti).
Система решается методом прогонки, в котором решение отыскивается в виде yi=αi+1yi+1+βi+1, i=0,1,…,M-1. Алгоритм состоит из следующих шагов.
1. Поскольку y0=α1y1+β1, и из первого уравнения y0 = (A - k2 y1/τ) / (k1—k2/τ), то α1 = - k2/ [ τ (k1—k2/τ) ], β1 = A/(k1—k2/τ)
2.
Следующие прогоночные
коэффициенты ai, bi определяются по рекурентным формулам (прямой ход
метода прогонки с трёхдиагональной матрицей):
αi+1=,
β i+1= i=1,2,...,M-1.
3. Для определения решения на правой границе отрезка воспользуемся последним уравнением системы (l1+l2/τ) yM - l2yM-1 /τ=B и соотношением yM-1=αM yM + βM, откуда yM определяется как решение системы двух линейных алгебраических уравнений.
4. Обратный ход метода прогонки для определения решения.
yi=αi+1yi+1+βi+1, i = M-1,M-2, ..., 0.
Теперь
напишем программу решения краевой задачи для ОДУ второго порядка с помощью
конечно-разностного метода. Это программа будет тестовой, то есть, содержит
функцию yt(t)
, определяющую точное решение.
Уравнение:
Точное решение:
ut(t) = cos(sin(t))
+ C1 sin(sin(t)), C1=(10 - cos(sin(1)))/sin(sin(1)).
Программа
для заданного константой M будет выдавать на экран таблицу со столбцами
i
ti
yi
u(ti)
|u(ti)-yi|,
также по последнему столбцу определяется максимальная погрешность.
Метод показывает первый порядок аппроксимации, то есть, при увеличении
числа узлов в 2 раза точность решения также повышается примерно в 2 раза.
Возможно, здесь чуть-чуть иначе вычисляются α1
и
β1
, но с посчитанной ниже задачей в Excel совпало.
program odu2test; {Решение тестовой задачи с известным точным решением yt(t)} const n=10; {Число точек} a=0; {Границы интервала} b=1; {Краевые условия} k1=1; k2=0; a1=1; l1=1; l2=0; b1=10; type vec=array [0..n] of real; var h,eps,epsmax:real; al,bet,t,u,aa,bb,cc,ff:vec; i,j:integer; {Функция коэффициентов при 1-й производной u'(t)} function p(t:real):real; begin p:=sin(t)/cos(t); end; {Функция коэффициентов при u(t)} function g(t:real):real; begin g:=sqr(cos(t)); end; {Функция правой части уравнения} function f(t:real):real; begin f:=0; end; {Точное решение yt(t)} function yt(t:real):real; var c1:real; begin c1:=(10-cos(sin(1)))/sin(sin(1)); yt:=cos(sin(t))+c1*sin(sin(t)); end; begin {Вычисление шага сетки} h:=(b-a)/n; {Формирование векторов коэффициентов для 3-диагональной СЛАУ} for i:=0 to n do begin t[i]:=a+h*i; aa[i]:=1-p(t[i])*h/2; bb[i]:=1+p(t[i])*h/2; cc[i]:=2-g(t[i])*h*h; ff[i]:=h*h*f(t[i]); end; {Краевое условие при t=a} al[1]:=k2/(k2-k1*h); bet[1]:=-(a1*h)/(k2-k1*h); {Прямой ход метода прогонки} for i:=1 to n-1 do begin al[i+1]:=bb[i]/(CC[i]-al[i]*aa[i]); bet[i+1]:=(aa[i]*bet[i]-ff[i])/(cc[i]-al[i]*aa[i]); end; {Решение на правой границе отрезка} u[n]:=(l2*bet[n]+b1*h)/(l2+h*l1-l2*al[n]); {Обратный ход метода прогонки} for i:=n-1 downto 0 do u[i]:=al[i+1]*u[i+1]+bet[i+1]; {Вывод таблицы и определение макс.погрешности} writeln('i':4,'t[i]':10,'y[i]':10,'yt(t[i])':10,'Погрешн.':10); for i:=0 to n do begin eps:=abs(u[i]-yt(t[i])); if eps>epsmax then epsmax:=eps; writeln (i:4,t[i]:10:4,u[i]:10:4,yt(t[i]):10:4,eps); end; writeln ('Максимальная погрешность=',epsmax); reset (input); readln; end.
Используя
эту программу, можно решить и другую аналогичную краевую задачу, изменив в
нашей программке всего несколько строчек. Например, такую:
Уравнение
, интервал [1.1,1.4], краевое условие
при t=a имеет
вид u-2u ′ =-1,
краевое условие
при t=b имеет вид u=4. Вроде бы, должно получиться так:
program odu2my; const n=10; {Число точек} a=1.1; {Границы интервала} b=1.4; {Краевые условия} k1=1; k2=-2; a1=-1; l1=1; l2=0; b1=4; type vec=array [0..n] of real; var h:real; al,bet,t,u,aa,bb,cc,ff:vec; i,j:integer; {Функция коэффициентов при 1-й производной u'(t)} function p(t:real):real; begin p:=-1; end; {Функция коэффициентов при u(t)} function g(t:real):real; begin g:=2/t; end; {Функция правой части уравнения} function f(t:real):real; begin f:=t+0.4; end; begin {Вычисление шага сетки} h:=(b-a)/n; {Формирование векторов коэффициентов для 3-диагональной СЛАУ} for i:=0 to n do begin t[i]:=a+h*i; aa[i]:=1-p(t[i])*h/2; bb[i]:=1+p(t[i])*h/2; cc[i]:=2-g(t[i])*h*h; ff[i]:=h*h*f(t[i]); end; {Краевое условие при t=a} al[1]:=k2/(k2-k1*h); bet[1]:=-(a1*h)/(k2-k1*h); {Прямой ход метода прогонки} for i:=1 to n-1 do begin al[i+1]:=bb[i]/(CC[i]-al[i]*aa[i]); bet[i+1]:=(aa[i]*bet[i]-ff[i])/(cc[i]-al[i]*aa[i]); end; {Решение на правой границе отрезка} u[n]:=(l2*bet[n]+b1*h)/(l2+h*l1-l2*al[n]); {Обратный ход метода прогонки} for i:=n-1 downto 0 do u[i]:=al[i+1]*u[i+1]+bet[i+1]; {Вывод таблицы и определение макс.погрешности} writeln('i':4,'t[i]':10,'y[i]':10); for i:=0 to n do writeln (i:4,t[i]:10:4,u[i]:10:4); reset (input); readln; end.
Ниже то же самое посчитано в Excel для 10, 20 и 40 узлов сетки, в файле серым фоном обозначены ячейки с формулами и константами, которые являются входными данными, жёлтым - формулы, которые считаются иначе, чем соседние формулы выше и/или ниже.
odu2.xls (56 Кб, Excel XP/2003)
гостевая; E-mail |