# Методы решения обыкновенных дифференциальных уравнений и систем обыкновенных дифференциальных уравнений, страница 2

0.0505

0.0551

0.0601

0.0653

0.0709

0.0768

0.0830

0.0896

0.0965

0.1039

0.1116

0.1198

0.1284

>> xx = 0:0.01:0.6;

>> pp =interp1(t,y,'cubic','pp');

>> yy = ppval(pp,xx);

>> plot(xx,yy,'-'),hold on;

>>

СОДУ:

function f1=f(x,y,z);

f1=z*x-y+1;

function f2=f(x,y,z);

f2=y-x^2;

function res = Euler(m,x0,y0,z0)

x=zeros(m,1);

y=zeros(m,1);

z=zeros(m,1);

y_m=zeros(m,1);

z_m=zeros(m,1);

x(1)=x0;

y(1)=y0;

z(1)=z0;

h=0.1;

for i=1:m-1

x(i+1)=x(1)+i*h;

end

'Method Euler'

for i=1:m-1

y(i+1)=y(i)+h*f1(x(i),y(i),z(i));

z(i+1)=z(i)+h*f2(x(i),y(i),z(i));

end

[x,y,z]

'Method Euler modified'

y_m(1)=y0;

z_m(1)=z0;

for i=1:m-1

y_m(i+1)=y_m(i)+h*f1(x(i)+h/2,y_m(i)+h/2*f1(x(i),y_m(i),z_m(i)),z_m(i)+h/2*(f2(x(i),y_m(i),z_m(i))));    z_m(i+1)=z_m(i)+h*f2(x(i)+h/2,y_m(i)+h/2*f1(x(i),y_m(i),z_m(i)),z_m(i)+h/2*(f2(x(i),y_m(i),z_m(i))));

end

[x,y_m,z_m]

function res = Runge(m,x0,y0,z0);

x=zeros(m,1);

y_RK=zeros(m,1);

z_RK=zeros(m,1);

y_RK(1)=y0;

z_RK(1)=z0;

x(1)=x0;

h=0.1;

for i=1:m-1

x(i+1)=x(1)+i*h;

end

for i=1:m-1

K1Y(i)=h*f1(x(i),y_RK(i),z_RK(i));

K1Z(i)=h*f2(x(i),y_RK(i),z_RK(i));

K2Y(i)=h*f1(x(i)+h/2,y_RK(i)+K1Y(i)/2,z_RK(i)+K1Z(i)/2);

K2Z(i)=h*f2(x(i)+h/2,y_RK(i)+K1Y(i)/2,z_RK(i)+K1Z(i)/2);

K3Y(i)=h*f1(x(i)+h/2,y_RK(i)+K2Y(i)/2,z_RK(i)+K2Z(i)/2);

K3Z(i)=h*f2(x(i)+h/2,y_RK(i)+K2Y(i)/2,z_RK(i)+K2Z(i)/2);

K4Y(i)=h*f1(x(i)+h,y_RK(i)+K3Y(i),z_RK(i)+K3Z(i));

K4Z(i)=h*f2(x(i)+h,y_RK(i)+K3Y(i),z_RK(i)+K3Z(i));

y_RK(i+1)=y_RK(i)+(K1Y(i)+2*K2Y(i)+2*K3Y(i)+K4Y(i))/6;

z_RK(i+1)=z_RK(i)+(K1Z(i)+2*K2Z(i)+2*K3Z(i)+K4Z(i))/6;

end

res = [x,y_RK,z_RK];

return

x=zeros(m,1);

y=zeros(m,1);

z=zeros(m,1);

h=0.1;

x(1)=x0;

y(1)=y0;

z(1)=z0;

for i=1:m-1

x(i+1)=x(1)+i*h;

end

v = Runge(8,0,1,2);

k=size(v);

k=k(1);

for i=2:k

y(i,1) = y(i,1)+v(i,2);

z(i,1) = z(i,1)+v(i,3);

end

for j=4:5

y(j+1)=y(j)+(h/24)*(55*f1(x(j),y(j),z(j))-59*f1(x(j-1),y(j-1),z(j-1))+37*f1(x(j-2),y(j-2),z(j-1))-9*f1(x(j-3),y(j-3),z(j-1)));

z(j+1)=z(j)+(h/24)*(55*f2(x(j),y(j),z(j))-59*f2(x(j-1),y(j-1),z(j-1))+37*f2(x(j-2),y(j-2),z(j-1))-9*f2(x(j-3),y(j-3),z(j-1)));

end

%'Mathod Miln'

for j=6:7

y12(j+1)=y(j-3)+(4*h/3)*(2*f1(x(j-2),y(j-2),z(j-2))-f1(x(j-1),y(j-1),z(j-1))+2*f1(x(j),y(j),z(j)));

y(j+1)=y(j-1)+(h/3)*(f1(x(j+1),y12(j+1),z(j+1))+4*f1(x(j),y(j),z(j))+f1(x(j-1),y(j-1),z(j-1)));

z12(j+1)=z(j-3)+(4*h/3)*(2*f2(x(j-2),y(j-2),z(j-2))-f2(x(j-1),y(j-1),z(j-1))+2*f2(x(j),y(j),z(j)));

z(j+1)=z(j-1)+(h/3)*(f2(x(j+1),y(j+1),z12(j+1))+4*f2(x(j),y(j),z(j))+f2(x(j-1),y(j-1),z(j-1)));

end

res = [x,y,z];

return

>> Euler(8,0,1,2)

ans =

Method Euler

ans =

0    1.0000    2.0000

0.1000    1.0000    2.1000

0.2000    1.0210    2.1990

0.3000    1.0629    2.2971

0.4000    1.1255    2.3944

0.5000    1.2087    2.4909

0.6000    1.3124    2.5868

0.7000    1.4364    2.6821

ans =

Method Euler modified

ans =

0    1.0000    2.0000

0.1000    1.0103    2.0998

0.2000    1.0405    2.1995

0.3000    1.0907    2.2993

0.4000    1.1608    2.3991

0.5000    1.2510    2.4990

0.6000    1.3611    2.5988

0.7000    1.4912    2.6986

v =

0    1.0000    2.0000

0.1000    1.0100    2.1000

0.2000    1.0400    2.2000

0.3000    1.0900    2.3000

0.4000    1.1615    2.4000

0.5000    1.2535    2.5004

0.6000    1.3611    2.6006

0.7000    1.4933    2.7007

Вывод: научилась находить решение дифференциальных уравнений с помощью методов Эйлера, модифицированного метода Эйлера, метода Адамса, метода Рунге-Кутта, метода Милна.