马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
這是Runge kutta的程式碼
function rs= rk4sysnew(X0,a,b,m)
%to return the approximation values of X(t)
%by RK4 method, X0 is the initial
%t is in[a,b]
n=length(X0);
rs=zeros(n,m) ;
K1=zeros(n,1) ; K2=zeros(n,1) ;
K3=zeros(n,1) ; K4=zeros(n,1) ;
t=a ;
X=X0 ;X2=zeros(n,1) ;X3=X2; X4=X2 ;
h=(b-a)/m ;
for k=1:m
%compute K1,K2,K3and K4
K1=fxsys(t,X);
for i=1:n
X2(i)=X(i)+1/2*h*K1(i);
end
K2=fxsys(t+h/2,X2);
for i=1:n
X3(i)=X(i)+1/2*h*K2(i);
end
K3=fxsys(t+h/2,X3);
for i=1:n
X4(i)=X(i)+h*K3(i);
end
K4=fxsys(t+h,X4);
%compute X(t+h)
for i=1:n
X(i)=X(i)+h/6*(K1(i)+2*K2(i)+2*K3(i)+K4(i));
rs(i,k)=X(i) ;
end %for i
t=t+h ;
end % for k
以及我的高階非線性方程式我是降成一階常微分方程式
function fx=fxsys(t,X)
fx=zeros(length(X),1);
fx(1)=X(2);
fx(2)=(-6.15*10.^6*(X(3)-X(1)))/13.2;
fx(3)=X(4);
fx(4)=((8.1*10.^7*(X(5)-X(3)))-(6.15*10.^6*(X(3)-y(1))))/(1.3);
fx(5)=X(6);
fx(6)=((6.03*10.^6*(X(5)-X(3)))-(1.545*10.^4*0.176*(0.176*X(7)-0.176*X(5))))/0.76;
fx(7)=X(8);
fx(8)=((1.545*10.^4*0.176*(0.176*X(7)-0.176*X(5)))+(1.545*10.^4*0.176*(0.176*X(9)-0.176*X(7))))/0.76;
fx(9)=X(10);
fx(10)=((1.545*10.^4*0.176*(0.176*X(9)-0.176*X(7)))-(8.1*10.^7*(X(11)-X(9))))/0.76;
fx(11)=X(12);
fx(12)=((8.1*10.^7*(X(11)-X(9)))+(1.618*10.^4*0.150*(0.150*X(13)-0.150*X(11))))/0.56;
fx(13)=X(14);
fx(14)=((1.618*10.^4*0.150*(0.150*X(12)-0.150*X(11)))-(6.03*10.^6*(X(15)-X(13))))/0.56;
fx(15)=X(16);
fx(16)=(6.03*10.^6*(X(15)-y(13)))/3.66;
fx(17)=X(18);
fx(18)=(-1.159*10.^4*0.156*(0.084*X(19)-0.156*X(17)))/0.25;
fx(19)=X(20);
fx(20)=(-2.01*10.^6*(X(21)-X(19)))/0.03;
fx(21)=X(22);
fx(22)=(2.01*10.^6*(X(21)-X(19)))/1.5;
可是當我執行時,給定邊界條件
a=0; b= 1; m= 10;
>> X0=[0; 0] ;
>> Rs =rk4sysnew(X0,a,b,m);
出現這樣的錯誤
??? Index exceeds matrix dimensions.
Error in ==> C:\MATLAB6p5\work\rk4sysnew.m (fxsys)
On line 38 ==> fx(2)=(-6.15*10.^6*(X(3)-X(1)))/13.2;
Error in ==> C:\MATLAB6p5\work\rk4sysnew.m
On line 14 ==> K1=fxsys(t,X);
是那裡錯???麻煩大俠邦解 |