程序的目的是:实现ADMAS求二阶微分方程
程序如下:
f1=sym('Dy=z');
f2=sym('Dz=-2*z-36*y+36');
x=0:0.01:10;
h=0.01;
y(1)=0 ;z(1)=0;
for n=1:4
k1=-2*z(n)-36*y(n)+36;
k2=-2*[z(n)+h/2*k1]-36*[y(n)+h/2*z(n)]+36;
k3=-2*[z(n)+h/2*k2]-36*[y(n)+h/2*z(n)+h^2/4*k2]+36;
k4=-2*[z(n)+h*k3]-36*[y(n)+h*z(n)+h^2*k2]+36;
z(n+1)=z(n)+h*[k1+2*k2+2*k3+k4];
y(n+1)=y(n)+h*z(n)+h^2*[k1+k2+k3];
end
for k=1:4
x(k)=x(1)+k*h;
Dy(k)=z(k);
Dz(k)=-2*z(k)-36*y(k)+36;
end
for n=5:length(x)-1
x(5)=x(4)+h;
Dzp=z(4)+h*[55*Dz(4)-59*Dz(3)+37*Dz(2)-9*Dz(1)]/24;
Dyp=y(4)+h*[55*Dy(4)-59*Dy(3)+37*Dy(2)-9*Dy(1)]/24;
z(5)=z(4)+h*[9*Dzp+19*Dz(4)-5*Dz(3)+Dz(2)]/24;
y(5)=y(4)+h*[9*Dyp+19*Dy(4)-5*Dy(3)+Dy(2)]/24;
Dz(5)=-2*z(5)-36*y(5)+36;
Dy(5)=z(5);
x(4)=x(5);
y(4)=y(5);
z(4)=z(5);
for j=1:4
Dz(j)=Dz(j+1);
Dy(j)=Dy(j+1);
end
end
plot(x,y,'b')
图形结果:

请问怎么上传图片啊,我不会

[此贴子已经被作者于2006-12-3 11:07:33编辑过]