马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
主要是利用四阶龙格——库塔公式求解的问题
利用龙格库塔方法计算ndetaz处的R的值
我编写的程序运行不出来
dudz=10;
dRdZ=dudz*0.12;
k=850;
detaz=1/k;
n=[0:1:850];
z=detaz*n;
z(1)=0;R(1)=0;A(1)=0;
z(2)=detaz;R(2)=dRdZ*detaz;
RR(2)=(R(2))^2+10*(R(2))^(4/3);
A(2)=2.54*detaz*RR(2)/2;
k1=sqrt(A(2))/detaz;
Rz(2)=R(2)+k1*0.5*detaz;
RRz(2)=(Rz(2))^2+10*(Rz(2))^(4/3);
B(2)=2.54*1.5*detaz*(RRz(2))/2;
k2=sqrt(B(2))/(1.5*detaz);
Rz1(2)=R(2)+k2*0.5*detaz;
RRz1(2)=(Rz1(2))^2+10*(Rz1(2))^(4/3);
C(2)=2.54*1.5*detaz*(RRz1(2))/2;
k3=sqrt(C(2))/(1.5*detaz);
R(3)=R(2)+k3*detaz;
RR(3)=(R(3))^2+10*(R(3))^(4/3);
D(3)=2.54*2*detaz*(RR(3))/2;
k4=sqrt(D(3))/(2*detaz);
R(3)=R(2)+detaz/6*(k1+2*k2+2*k3+k4);
while n>=2
i=2:1:849
zz=[0:detaz/2:detaz*i];
R=[R(1:(i+1))];
Re=spline(z,R,zz);
A(i+1)=quad(zz,Re);
k1=sqrt(A(i+1))/detaz;
Rz(i+1)=R(i+1)+k1*0.5*detaz;
zz1=[0:detaz/2:detaz*((i+1)+0.5)];
R1=[R(1:(i+1)),Rz(i+1)];
Re1=spline(z,R1,zz1);
B(i+1)=quad(zz1,Re1);
k2=sqrt(B(i+1))/(1.5*detaz);
Rz1(i+1)=R(i+1)+k2*0.5*detaz;
zz1=[0:detaz/2:detaz*((i+1)+0.5)];
R2=[R(1:(i+1)),Rz1(i+1)];
Re2=spline(z,R1,zz1);
C(i+1)=quad(zz1,Re2);
k3=sqrt(C(i+1))/1.5/detaz;
R(i+2)=R(i+1)+k3*detaz;
zz2=[0:detaz/2:detaz*(i+2)];
R3=[R(1:(i+1)),R(i+2)];
Re3=spline(z,R3,zz2);
D(i+1)=quad(zz2,Re3);
k4=sqrt(D(i+1))/2/detaz;
R(i+2)=R(i+1)+detaz/6*(k1+2*k2+2*k3+k4);
end
plot(z,R)
运行出现下面的提示:
Index exceeds matrix dimensions |