Примеры решений Примеры реализации численных методов решения систем дифференциальных уравнений
Рассмотрим некоторые численные методы решения заданной системы.
Программная реализация (Turbo Pascal)
program Euler;
var T : Text;
a, b, h, x: real;
y1, y2, y3, y4, y1_1, y1_2, y1_3, y1_4 : real;
function F(x, y1, y2, y3, y4: real; n:byte):real;
begin
case n of
1: f:=-y2;
2: f:=y1;
3: f:=y1-y4;
4: f:=y2+y3;
end;
end;
begin
assign(T, 'c:\euler.txt'); rewrite(T);
a:=0; b:=2*pi; h:=pi/40;
y1:=1; y2:=0; y3:=0; y4:=0;
x:=a;
while x<=b+h do
begin
writeln(T, x:10:5, y1:10:5, y2:10:5, y3:10:5, y4:10:5);
y1_1:=y1+h*f(x, y1, y2, y3, y4, 1);
y1_2:=y2+h*f(x, y1, y2, y3, y4, 2);
y1_3:=y3+h*f(x, y1, y2, y3, y4, 3);
y1_4:=y4+h*f(x, y1, y2, y3, y4, 4);
y1:=y1_1; y2:=y1_2; y3:=y1_3; y4:=y1_4;
x:=x+h;
end;
flush(T); close(T);
end.
Программная реализация (Turbo Pascal)
program Euler_Koshi;
var T:Text;
x,h,a,b,y1, y2, y3, y4, y1_1, y1_2, y1_3,y1_4,
y11, y12, y13, y14:real;
function F(x, y1, y2, y3, y4: real; n:byte):real;
begin
case n of
1: f:=-y2;
2: f:=y1;
3: f:=y1-y4;
4: f:=y2+y3;
end;
end;
begin
assign(T, 'd:\euler_k.txt');
rewrite(T);
a:=0; b:=2*pi; h:=pi/160;
y1:=1; y2:=0; y3:=0; y4:=0;
x:=a;
while x<=b+h do
begin
writeln(T, x:10:5, y1:10:5, y2:10:5, y3:10:5, y4:10:5,
cos(x):10:5, sin(x):10:5, x*cos(x):10:5, x*sin(x):10:5);
y1_1:=y1+h*f(x, y1, y2, y3, y4, 1);
y1_2:=y2+h*f(x, y1, y2, y3, y4, 2);
y1_3:=y3+h*f(x, y1, y2, y3, y4, 3);
y1_4:=y4+h*f(x, y1, y2, y3, y4, 4);
y11:=y1+h/2*(f(x, y1, y2, y3, y4, 1)+f(x+h,y1_1, y1_2, y1_3, y1_4, 1));
y12:=y2+h/2*(f(x, y1, y2, y3, y4, 2)+f(x+h,y1_1, y1_2, y1_3, y1_4, 2));
y13:=y3+h/2*(f(x, y1, y2, y3, y4, 3)+f(x+h,y1_1, y1_2, y1_3, y1_4, 3));
y14:=y4+h/2*(f(x, y1, y2, y3, y4, 4)+f(x+h,y1_1, y1_2, y1_3, y1_4, 4));
y1:=y11; y2:=y12; y3:=y13; y4:=y14;
x:=x+h;
end;
flush(T); close(T);
end.
Программная реализация (Turbo Pascal)
program rk_III;
var T:Text;
x,h,a,b,y1, y2, y3, y4,
k1_1, k1_2, k1_3, k1_4,
k2_1, k2_2, k2_3, k2_4,
k3_1, k3_2, k3_3, k3_4:real;
function F(x, y1, y2, y3, y4: real; n:byte):real;
begin
case n of
1: f:=-y2;
2: f:=y1;
3: f:=y1-y4;
4: f:=y2+y3;
end;
end;
begin
assign(T, 'd:\rk_iii.txt');
rewrite(T);
a:=0; b:=2*pi; h:=pi/160;
y1:=1; y2:=0; y3:=0; y4:=0;
x:=a;
while x<=b+h do
begin
writeln(T, x:10:5, y1:10:5, y2:10:5, y3:10:5, y4:10:5,
cos(x):10:5, sin(x):10:5, x*cos(x):10:5, x*sin(x):10:5);
k1_1:=h*f(x, y1, y2, y3, y4, 1);
k1_2:=h*f(x, y1, y2, y3, y4, 2);
k1_3:=h*f(x, y1, y2, y3, y4, 3);
k1_4:=h*f(x, y1, y2, y3, y4, 4);
k2_1:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 1);
k2_2:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 2);
k2_3:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 3);
k2_4:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 4);
k3_1:=h*f(x+h,y1+k2_1, y2+k2_2, y3+k2_3, y4+k2_4, 1);
k3_2:=h*f(x+h,y1+k2_1, y2+k2_2, y3+k2_3, y4+k2_4, 2);
k3_3:=h*f(x+h,y1+k2_1, y2+k2_2, y3+k2_3, y4+k2_4, 3);
k3_4:=h*f(x+h,y1+k2_1, y2+k2_2, y3+k2_3, y4+k2_4, 4);
y1:=y1+1/4*(k1_1+2*k2_1+k3_1);
y2:=y2+1/4*(k1_2+2*k2_2+k3_2);
y3:=y3+1/4*(k1_3+2*k2_3+k3_3);
y4:=y4+1/4*(k1_4+2*k2_4+k3_4);
x:=x+h;
end;
flush(T); close(T);
end.
Программная реализация (Turbo Pascal)
program rk_iv;
var T:Text;
x,h,a,b,y1, y2, y3, y4,
k1_1, k1_2, k1_3, k1_4,
k2_1, k2_2, k2_3, k2_4,
k3_1, k3_2, k3_3, k3_4,
k4_1, k4_2, k4_3, k4_4:real;
function F(x, y1, y2, y3, y4: real; n:byte):real;
begin
case n of
1: f:=-y2;
2: f:=y1;
3: f:=y1-y4;
4: f:=y2+y3;
end;
end;
begin
assign(T, 'd:\rk_iv.txt');
rewrite(T);
a:=0; b:=2*pi; h:=pi/80;
y1:=1; y2:=0; y3:=0; y4:=0;
x:=a;
while x<=b+h do
begin
writeln(T, x:10:5, y1:10:5, y2:10:5, y3:10:5, y4:10:5,
cos(x):10:5, sin(x):10:5, x*cos(x):10:5, x*sin(x):10:5);
k1_1:=h*f(x, y1, y2, y3, y4, 1);
k1_2:=h*f(x, y1, y2, y3, y4, 2);
k1_3:=h*f(x, y1, y2, y3, y4, 3);
k1_4:=h*f(x, y1, y2, y3, y4, 4);
k2_1:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 1);
k2_2:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 2);
k2_3:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 3);
k2_4:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 4);
k3_1:=h*f(x+h/2, y1+k2_1/2, y2+k2_2/2, y3+k2_3/2, y4+k2_4/2, 1);
k3_2:=h*f(x+h/2, y1+k2_1/2, y2+k2_2/2, y3+k2_3/2, y4+k2_4/2, 2);
k3_3:=h*f(x+h/2, y1+k2_1/2, y2+k2_2/2, y3+k2_3/2, y4+k2_4/2, 3);
k3_4:=h*f(x+h/2, y1+k2_1/2, y2+k2_2/2, y3+k2_3/2, y4+k2_4/2, 4);
k4_1:=h*f(x+h,y1+k3_1, y2+k3_2, y3+k3_3, y4+k3_4, 1);
k4_2:=h*f(x+h,y1+k3_1, y2+k3_2, y3+k3_3, y4+k3_4, 2);
k4_3:=h*f(x+h,y1+k3_1, y2+k3_2, y3+k3_3, y4+k3_4, 3);
k4_4:=h*f(x+h,y1+k3_1, y2+k3_2, y3+k3_3, y4+k3_4, 4);
y1:=y1+1/6*(k1_1+2*k2_1+2*k3_1+k4_1);
y2:=y2+1/6*(k1_2+2*k2_2+2*k3_2+k4_2);
y3:=y3+1/6*(k1_3+2*k2_3+2*k3_3+k4_3);
y4:=y4+1/6*(k1_4+2*k2_4+2*k3_4+k4_4);
x:=x+h;
end;
flush(T); close(T);
end.