Методы решения обыкновенных дифференциальных уравнений и систем обыкновенных дифференциальных уравнений, страница 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

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

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

%'Method Adams'

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 = AdamsMilnSODE(8,0,1,2)

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

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