Порядок аппроксимации производных.
Численное решение краевых задач для обыкновенного дифференциального уравнения (ОДУ) второго порядка
Общий вид краевой задачи на отрезке для ОДУ второго порядка.
Краевая задача на отрезке для дифференциального уравнения второго порядка в общем случае представима в виде
(1.1)
Здесь все величины, кроме , предполагаются заданными.
Отметим, что существует некоторое бесконечное множество функций, удовлетворяющих приведенному дифференциальному уравнению, и лишь наличие области изменения ( – геометрическая характеристика конструкции) и условий на краях области обеспечивает единственность его решения. Однако и это решение (хотя оно и единственное), в общем случае, нельзя представить в аналитической форме. Поэтому, как правило, задача решается численным методом.
Метод конечных разностей.
Метод конечных разностей (МКР) – это наиболее простой и при этом достаточно эффективный способ численного решения краевой задачи. Его суть кратко опишем ниже.
Разобьем отрезок на отрезков (рис. 1.1).
Рис. 1.1. Схема аппроксимации по методу конечных разностей.
Введем обозначения: – координаты точек разбиения (узлов); – номер точки разбиения ( );
(1.2)
– длина -го отрезка (шаг разбиения);
(1.3)
– «средний» шаг;
, при этом , ; (1.4)
; ; . (1.5)
Производные в -ой точке заменим разностными соотношениями:
(1.6)
– правая разность;
(1.7)
– левая разность;
(1.8)
– центральная разность.
При в соответствии с известным определением производной все три величины будут стремиться к . Использование той или другой из них зависит от конкретной ситуации.
Вторая производная в -ой точке может быть приближенно представлена разностным отношением в виде разности от первых (вторая производная – это первая производная от первой производной), например
. (1.9)
При формулу второй разности можно упростить:
. (1.10)
Пользуясь приведенными обозначениями и формулами, можно представить задачу (1.1) в каждой -ой точке относительно величин следующим образом:
(1.11)
Это и есть разностный аналог краевой задачи (1.1). Здесь неизвестных и линейных уравнений. Приводя подобные члены, получим систему линейных алгебраических уравнений относительно неизвестных величин :
(1.12)
или в матричном виде
, (1.13)
где
; ; ,
(1.14)
при этом элементы матрицы и вектора правой части определяются формулами
, ; , ,
, ; , ; (1.15)
; , ; . (1.16)
Следует отметить, что матрица трехдиагональная.
После решения системы (1.13) получим приближенное решение задачи (1.1) – значения искомой функции в -х точках. Соединяя ломаной линией, получаем приближенное решение во всех точках области .
Точность решения задачи (1.1) зависит от точности аппроксимации значений производных разностями, которая, в свою очередь, зависит от величины шага разбиения.
Порядок аппроксимации производных.
Для упрощения оценок предположим, что . По формуле Тейлора получим:
,
где . (1.17)
Тогда для правой разности получим
т.е. окончательно
(1.18)
и имеем аппроксимацию первого порядка.
Аналогично для левой разности
,
т.е. также получаем результат (1.18) и первый порядок аппроксимации.
Рассмотрим также центральную разность.
т.е. (1.19)
и здесь имеем аппроксимацию второго порядка.
Следовательно, центральная разность является, как правило, более точной аппроксимацией первой производной.
Установим порядок аппроксимации для второй производной. Имеем:
т.е. . (1.20)
Итак, имеем аппроксимацию второго порядка. Аналогично можно установить порядок аппроксимации для любой разностной формулы.
1.4. Разностный аналог краевой задачи для случая .
В этом случае вид системы (1.11) упрощается, а именно
(1.21)
Приводя подобные члены, получим
(1.22)
Пример. Методом конечных разностей решить краевую задачу
В этом случае
, ,
; ; ,
, , , , , .
Пример М-файла
function krz_odu2
a=input('ввести координату левого конца отрезка a=');
b=input('ввести координату правого конца отрезка b=');
al1=input('ввести параметр alfa1='); bt1=input('ввести параметр beta1=');
gm1=input('ввести параметр gamma1=');
al2=input('ввести параметр alfa2='); bt2=input('ввести параметр beta2=');
gm2=input('ввести параметр gamma2=');
n=input('ввести число точек n=');
l=b-a; h=l/(n-1); xi=(a:h:b)'; x=xi(2:n-1);
P=3*x;
q=2+abs(x);
f=2-x;
tm=1/h^2-P/(2*h);
t0=q-2/h^2;
tp=1/h^2+P/(2*h);
A=diag([tm;-bt2/h],-1)+diag([al1-bt1/h;t0;al2+bt2/h])+diag([bt1/h;tp],1);
disp('матрица системы разностных уравнений'),disp(A)
b=[gm1;f;gm2];
disp('вектор правой части системы разностных уравнений'),disp(b)
y=A\b;
disp('решение:'),disp(y)
plot(xi,y),grid on,title('reshenie krajevoj zadachi y=y(x)')
Результаты решения при n=9
ввести координату левого конца отрезка a=-3
ввести координату правого конца отрезка b=3
ввести параметр alfa1=1
ввести параметр beta1=2
ввести параметр gamma1=0
ввести параметр alfa2=1
ввести параметр beta2=0
ввести параметр gamma2=2
ввести число точек n=9
матрица системы разностных уравнений
-1.67 2.67 0 0 0 0 0 0 0
6.28 0.69 -2.72 0 0 0 0 0 0
0 4.78 -0.06 -1.22 0 0 0 0 0
0 0 3.28 -0.81 0.28 0 0 0 0
0 0 0 1.78 -1.56 1.78 0 0 0
0 0 0 0 0.28 -0.81 3.28 0 0
0 0 0 0 0 -1.22 -0.06 4.78 0
0 0 0 0 0 0 -2.72 0.69 6.28
0 0 0 0 0 0 0 0 1.00
вектор правой части системы разностных уравнений
4.25
3.5
2.75
1.25
0.5
-0.25
решение:
-1.0629
-0.66432
-4.1819
-5.2704
43.962
44.862
7.6812
11.67
>>
Результаты решения при n=17
ввести координату левого конца отрезка a=-3
ввести координату правого конца отрезка b=3
ввести параметр alfa1=1
ввести параметр beta1=2
ввести параметр gamma1=0
ввести параметр alfa2=1
ввести параметр beta2=0
ввести параметр gamma2=2
ввести число точек n=17
матрица системы разностных уравнений
………………………………………………
вектор правой части системы разностных уравнений
…………………………………………………………..
решение:
0.58
0.47
0.32
0.09
-0.14
0.81
5.26
13.34
20.52
22.21
17.38
10.81
6.50
4.44
3.32
2.56
2.00
>>
Результаты решения при n=31
ввести координату левого конца отрезка a=-3
ввести координату правого конца отрезка b=3
ввести параметр alfa1=1
ввести параметр beta1=2
ввести параметр gamma1=0
ввести параметр alfa2=1
ввести параметр beta2=0
ввести параметр gamma2=2
ввести число точек n=31
матрица системы разностных уравнений
………………………………………………
вектор правой части системы разностных уравнений
…………………………………………………………..
решение:
0.57
0.51
0.44
0.36
0.27
0.16
0.07
0.06
0.30
1.06
2.65
5.33
9.06
13.41
17.60
20.76
22.34
21.95
19.83
16.65
13.24
10.21
7.85
6.15
4.96
4.12
3.50
3.01
2.62
2.28
2.00
>>
Варианты заданий
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.