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
Вывод: научилась находить решение дифференциальных уравнений с помощью методов Эйлера, модифицированного метода Эйлера, метода Адамса, метода Рунге-Кутта, метода Милна.
Уважаемый посетитель!
Чтобы распечатать файл, скачайте его (в формате Word).
Ссылка на скачивание - внизу страницы.