Дифференциальное уравнение (ДУ) выражает зависимость между функцией и ее производными. Решением ДУ является функция. Численные методы позволяют получить функцию в виде таблицы.
В библиотеке рассмотрены только обыкновенные ДУ и их системы. Уравнения в частных производных, как правило, требуют более мощной вычислительной машины.
Численное решение дифференциальных уравнений и их систем может осуществляться как одношаговыми, так и многошаговыми методами. Наиболее удачный из одношаговых - метод Рунге - Кутта 4-го порядка.
В тех случаях, когда уравнение сложное, а решение необходимо получить на протяженном интервале, используют многошаговые методы, которые обеспечивают ту же точность при на порядок большем шаге. Для запуска многошагового используется одношаговый метод.
y'=f(x,y);
начальные условия: (x0,y0).
Программа 110. Метод Рунге - Кутта 4-го порядка.
k1=sf(xi, yi);
k2=sf(xi+s/2, yi+k1/2);
k3=sf(xi+s/2, yi+k2/2);
k4=sf(xi+s, yi+k3);
yi+1=yi+(k1+2k2+2k3+k4)/6;
где s - шаг интегрирования.
10 INPUT "s/2",A,"x",B,"y",C
20 D=C:GOSUB 50:E=F:B=B+A:D=C+F:GOSUB 50:E=E+F+F
30 D=C+F:GOSUB 50:E=E+F+F:B=B+A:D=C+F+F:GOSUB 50
40 C=C+(E+F)/3:PRINT "x";B,"y";C:GOTO 20
50 F=A*(...f{B,D}...):RETURN
Размер: 142+, F
Пример: y'=x/y; s/2=0,05; x0=0; y0=1.
В программе: 50 F=A*B/D:RETURN
Ответ: ... f(1)=1,41421356 ...
Программа 111. Метод Милна 6-го порядка. Решение устойчиво лишь при y'>1.
Предикатор: ypi+1=yi-5+3s/10{11(y'i+y'i-4)-14(y'i-1+y'i-3)+26y'i-2}.
Корректор: yi+1=yi-3+2s/45{7(y'i+1+y'i-3)+32(y'i+y'i-2)+12y'i-1}.
Для получения первых пяти отсчетов y и y' используется метод Рунге - Кутта 4-го порядка с шагом, уменьшенным в 10 раз.
10 INPUT "s",A,"x",B,"y",C:A=A/20:G=C:GOSUB
90
20 FOR F=.1 TO 5 STEP .1:D=C:E=Q:B=B+A:C=D+Q:GOSUB 90:E=E+Q+Q
30 C=D+Q:GOSUB 90:E=E+Q+Q:B=B+A:C=D+Q+Q:GOSUB 90:C=D+(E+Q)/3
40 GOSUB 90:IF FRAC F=0;PRINT "x";B,"y";C:G(F)=C:L(F)=20*Q
50 NEXT F:A=20*A
60 C=G+.3*(11*(Q+M)-14*(P+N)+26*O):FOR E=0 TO 9:G(E)=H(E)
70 NEXT E:B=B+A:GOSUB 90:L=H+2/45*(7*(Q+M)+32*(P+N)+12*O)
80 PRINT "x";B,"y";L:GOTO 60
90 Q=A*(...f{B,C}...):RETURN
Размер: 313+, Q
Пример: y'=x/y; s=0,1; x0=0; y0=1.
В программе: 90 Q=A*B/C:RETURN
Ответ: ... f(2)=2,23606797 ...
Программа 112. Метод Хемминга 4-го порядка.
Предикатор: ypi+1=yi+3+4s/(32y'i-y'i-1+2y'i-2).
Модификатор: ymi+1=ypi+1+112/121(yci-ypi).
Корректор: yci+1=(9yi-yi-2)/8+3s/8(y'mi+1+2y'i-y'i-1).
yi+1=yci+1-9/121(yci+1-ypi+1).
Для получения первых трех отсчетов y и y' используется метод Рунге - Кутта 4-го порядка с шагом, уменьшенным в 10 раз.
10 INPUT "s",A,"x",B,"y",C:A=A/20:G=C:GOSUB
90
20 FOR F=.1 TO 3 STEP .1:D=C:E=M:B=B+A:C=D+M:GOSUB 90:E=E+M+M
30 C=D+M:GOSUB 90:E=E+M+M:B=B+A:C=D+M+M:GOSUB 90:C=D+(E+M)/3
40 GOSUB 90:IF FRAC F=0;PRINT "x";B,"y";C:G(F)=C:J(F)=20*M
50 NEXT F:A=20*A:D=0
60 F=G+4/3*(M+M-L+K+K):C=F+112/121*D:FOR E=0 TO 5:G(E)=H(E)
70 NEXT E:B=B+A:GOSUB 90:J=(9*I-G)/8+3/8*(M+L+L-K):D=J-F
80 J=J-9/121*D:PRINT "x";B,"y";J:GOSUB 90:GOTO 60
90 M=A*(...f{B,C}...):RETURN
Размер: 334+, M
Пример: y'=x/y; s=0,1; x0=0; y0=1.
В программе: 90 M=A*B/C:RETURN
Ответ: ... f(2)=2,23606797 ...
y"=f(x,y,y');
начальные условия: (x0,y0,y'0).
Программа 113. Метод Рунге - Кутта 4-го порядка.
k1=sf(xi, yi,
y'i);
k2=sf(xi+s/2, yi+y'is/2,
y'i+k1/2);
k3=sf(xi+s/2, yi+y'is/2+k1s/2,
y'i+k2/2);
k4=sf(xi+s, yi+y'is+k2s/2,
y'i+k3);
y'i+1=y'i+(k1+2k2+2k3+k4)/6
yi+1=yi+s(y'i+(k1+k2+k3)/6)
где s - шаг интегрирования.
10 INPUT "s/2",A,"x",B,"y",C,"y'",D
20 E=C:F=D:GOSUB 60:G=J:B=B+A:C=E+F*A:D=F+G:GOSUB 60:H=J
30 C=C+G*A:D=F+H:GOSUB 60:I=J:B=B+A:C=E+2*A*(F+H):D=F+I+I
40 GOSUB 60:C=E+2*A*(F+(G+H+I)/3):D=F+(G+H+H+I+I+J)/3
50 PRINT "x";B,"y";C,"y'";D:GOTO 20
60 J=A*(...f{B,C,D}...):RETURN
Размер: 211+, J
Пример: y"=-y'/x-y; s/2=0,05; x0=10-30; y0=1; y'0=0.
В программе: 60 J=-A*(D/B+C):RETURN
Ответ: Отсчеты функции Бесселя J0(x), в частности при x=1
y=0,764548247 и y'=-0,439685328.
При интегрировании с шагом s=0,01 y(1)=0,765191315,
при s=0,001 y(1)=0,76519768.
Точное значение J0(1)=0,7651976866.
y(n)=f(x, y, y', ... , y(n-1));
начальные условия: (x0, y0,y'0, ... , y(n-1)0).
Программа 114. Метод Рунге - Кутта 4-го порядка. Обыкновенное дифференциальное уравнение любого порядка сводится к системе ДУ 1-го порядка посредством замены:
y1=y; yj+1=y'j, j=1...n-1;
где n - порядок уравнения;
y - искомая функция;
yj - переменные системы (см. описание программы 115).
Программа 114 является модификацией 115 с учетом эквивалентности yj+1=y'j. Текст на бейсике записан для уравнения второго порядка. В остальных случаях пользуйтесь представленной ниже таблицей соответствия переменных.
|
n |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
|
y0j |
G(C) |
H(C) |
I(C) |
J(C) |
K(C) |
L(C) |
M(C) |
N(C) |
|
еkij |
I(C) |
K(C) |
M(C) |
O(C) |
Q(C) |
S(C) |
U(C) |
W(C) |
|
y(n) |
G |
H |
I |
J |
K |
L |
M |
N |
|
y...y(n-1) |
E...F |
E...G |
E...H |
E...I |
E...J |
E...K |
E...L |
E...M |
|
Всего |
K |
N |
Q |
T |
W |
Z |
28 |
31 |
Изменению подлежат имена и константы, выделенные в программе нежирным шрифтом.
10 INPUT "s",A,"x",B:FOR C=1 TO
2:PRINT C;:INPUT D(C):NEXT
C
20 FOR C=1 TO 2:G(C)=D(C):I(C)=0:NEXT C
30 FOR D=0 TO 1.5 STEP .5:B=B+A*FRAC D
40 G=...f{B,E,F,...}...
50 FOR C=1 TO 2:I(C)=I(C)+A*E(C)*(1+SGN FRAC (D/1.5))
60 D(C)=G(C)+A*E(C)*(1+INT D)/2:NEXT C:NEXT D
70 FOR C=1 TO 2:D(C)=G(C)+I(C)/6:NEXT C:PRINT B;E:GOTO
20
Размер: 206+, K+
Пример: y"=-y'/x-y; s=0,1; x0=10-30; y0=1; y'0=0.
В программе: 40 G=-F/B-E
Ответ: Отсчеты функции Бесселя J0(x), в частности при x=1 y=0,764548247.
y'j=fj(x, y1, y2, ... , yn);
начальные условия:(x0, y01 ,y02, ... , y0n).
Программа 115. Метод Рунге - Кутта 4-го порядка.
Алгоритм в векторной форме записи:
k1=sF(xi, yi);
k2=sF(xi+s/2, yi+k1/2);
k3=sF(xi+s/2, yi+k2/2);
k4=sF(xi+s, yi+k3);
yi+1=yi+(k1+2k2+2k3+k4)/6;
где s - шаг интегрирования.
Текст программы записан для системы из двух уравнений. В остальных случаях пользуйтесь представленной ниже таблицей соответствия переменных.
Имена переменных в зависимости от порядка системы
|
n |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
|
y'j |
F(C) |
G(C) |
H(C) |
I(C) |
J(C) |
K(C) |
L(C) |
M(C) |
|
y0j |
H(C) |
J(C) |
L(C) |
N(C) |
P(C) |
R(C) |
T(C) |
V(C) |
|
еkij |
J(C) |
M(C) |
P(C) |
S(C) |
V(C) |
Y(C) |
Z(2+C) |
Z(5+C) |
|
y'1...y'n |
G...H |
H...J |
I...L |
J...N |
K...P |
L...R |
M...T |
N...V |
|
y1...yn |
E...F |
E...G |
E...H |
E...I |
E...J |
E...K |
E...L |
E...M |
|
Всего |
L |
P |
T |
X |
27 |
31 |
35 |
39 |
Изменению подлежат имена и константы, выделенные в программе нежирным шрифтом.
10 INPUT "s",A,"x",B:FOR C=1 TO
2:PRINT C;:INPUT D(C):NEXT
C
20 FOR C=1 TO 2:H(C)=D(C):J(C)=0:NEXT C
30 FOR D=0 TO 1.5 STEP .5:B=B+A*FRAC D
40 G=...f1{B,E,F,...}...
41 H=...f2{B,E,F,...}...
50 FOR C=1 TO 2:J(C)=J(C)+A*F(C)*(1+SGN FRAC (D/1.5))
60 D(C)=H(C)+A*F(C)*(1+INT D)/2:NEXT C:NEXT D:PRINT "x";B
70 FOR C=1 TO 2:D(C)=H(C)+J(C)/6:PRINT C;D(C):NEXT C:GOTO
20
Размер: 219+, L+
Пример:
y'1=y2
y'2=-(y2/x+y1)
s=0,1; x=10-30; y0=(1, 0)T
В программе:
40 G=F
41 H=-F/B-E
Ответ: Отсчеты функции Бесселя J0(x), в частности при x=1 y1=0,764548247 и y2=-0,439685328.
Программа 116. Метод Рунге - Кутта - Мерсона. Одношаговый метод 4-го порядка с грубой оценкой локальной ошибки, благодаря которой алгоритм автоматически изменяет шаг интегрирования.
Расчетные формулы в векторной записи:
k0=sF(xi, yi);
k1=sF(xi+s/3, yi+k0/3);
k2=sF(xi+s/3, yi+k0/6+k1/6);
k3=sF(xi+s/2, yi+k0/8+3k2/8);
k4=sF(xi+s, yi+k0/2-3k2/2+2k3);
yi+1=yi+(k0+4k3+k4)/6
Управление шагом интегрирования s:
Локальная ошибка r=-2k0+9k2-8k3+k4.
Если существует j, для которого |rj|>30e,
тогда s=s/2;
иначе, если для каждого j |rj|<e,
тогда s=s*2;
i=i+1, x=x+s.
Здесь e - допустимая локальная ошибка.
Текст программы записан для системы из двух уравнений. В остальных случаях пользуйтесь представленной ниже таблицей соответствия переменных.
10 INPUT "e",E,"s",A,"x",B
20 FOR C=1 TO 2:PRINT
C;:INPUT J(C):NEXT C
30 FOR C=1 TO 2:F(C)=J(C):NEXT C:GOSUB 150:FOR C=1 TO 2
40 L(C)=A*H(C):F(C)=J(C)+L(C)/3:NEXT C:B=B+A/3:GOSUB
150
50 FOR C=1 TO 2:F(C)=J(C)+(L(C)+A*H(C))/6:NEXT
C:GOSUB 150
60 FOR C=1 TO 2:N(C)=A*H(C):F(C)=J(C)+L(C)/8+.375*N(C)
70 NEXT C:D=0:B=B+A/6:GOSUB 150:FOR C=1 TO 2:P(C)=A*H(C)
80 F(C)=J(C)+L(C)/2-1.5*N(C)+2*P(C):NEXT
C:B=B+A/2:GOSUB 150
90 FOR C=1 TO 2:R(C)=A*H(C):F(C)=J(C)+(L(C)+4*P(C)+R(C))/6
100 F=R(C)-2*L(C)+9*N(C)-8*P(C):IF ABS F<E;D=D+1
110 IF ABS F>30*E;B=B-A:A=A/2:GOTO 30
120 NEXT C:IF D=2;A=A+A
130 PRINT "x";B:FOR C=1 TO 2:J(C)=F(C)
140 PRINT C;J(C):NEXT C:GOTO 30
150 I=...f1{B,G,H,...}...
151 J=...f2{B,G,H,...}...
160 RETURN
Размер: 500+, T+
Имена переменных в зависимости от порядка системы
|
n |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
|
y'j |
H(C) |
I(C) |
J(C) |
K(C) |
L(C) |
M(C) |
N(C) |
|
y0j |
J(C) |
L(C) |
N(C) |
P(C) |
R(C) |
T(C) |
V(C) |
|
k0j |
L(C) |
O(C) |
R(C) |
U(C) |
X(C) |
Z(1+C) |
Z(4+C) |
|
k2j |
N(C) |
R(C) |
V(C) |
Z(C) |
Z(4+C) |
Z(8+C) |
Z(12+C) |
|
k3j |
P(C) |
U(C) |
Z(C) |
Z(5+C) |
Z(10+C) |
Z(15+C) |
Z(20+C) |
|
k4j |
R(C) |
X(C) |
Z(4+C) |
Z(10+C) |
Z(16+C) |
Z(22+C) |
Z(28+C) |
|
y'1...y'n |
I...J |
J...L |
K...N |
L...P |
M...R |
N...T |
O...V |
|
y1...yn |
G...H |
G...I |
G...J |
G...K |
G...L |
G...M |
G...N |
|
Всего |
T |
26 |
33 |
40 |
47 |
54 |
61 |
Изменению подлежат имена и константы, выделенные в программе нежирным шрифтом
Пример:
y'1=y2
y'2=-(y2/x+y1)
e=10-8, s=0,01; x=10-30;
y0=(1, 0)T
В программе: 150 I=H:J=-H/B-G:RETURN
Контроль управления шагом (значения x): 6,25 10-4; 1,25 10-3; 2,5 10-3; 5 10-3; 0,01; 0,02; 0,04; ...
Ответ: Отсчеты функции Бесселя J0(x), в частности при x=1 y1=0,76519768 и y2=-0,440050579.