求助解常微分方程 四阶龙格库塔法

leadershuang

设：a1=y[0]; a2=y[1]; b1=y[2]; b2=y[3]; da1=dy[0]=y[4]; da2=dy[1]=y[5]; db1=dy[2]=y[6]; db2=dy[3]=y[7];
需要求解的四个式子我列在下边了，我用c语言编写的。其中c 、lanmta1、lanmta2、fai、pi为参数。我用c程序调用龙格库塔法求解出来数据绘图时总是不合适。哪位高手用matlab帮我编一下这个四阶龙格库塔法求解这四个式子！ 多谢！

dy[4]=16*c*y[5]/3-(lanmta1-fai)*pi*pi*y[0]
  -lanmta1*(3*pi*pi*pi*pi*y[0]*y[0]*y[0]/8+3*pi*pi*pi*y[0]*y[1]+3*pi*pi*pi*pi*y[0]*y[1]*y[1]
  +3*pi*pi*pi*pi*y[0]*y[2]*y[2]/8+pi*pi*pi*y[2]*y[3]+2*pi*pi*pi*pi*y[1]*y[2]*y[3]+pi*pi*pi*pi*y[0]*y[3]*y[3]);
dy[5]=-16*c*y[4]/3-4*(lanmta1-fai)*pi*pi*y[1]
  -lanmta1*(3*pi*pi*pi*y[0]*y[0]/2+3*pi*pi*pi*pi*y[0]*y[0]*y[1]+6*pi*pi*pi*pi*y[1]*y[1]*y[1]
  +pi*pi*pi*y[2]*y[2]/2+pi*pi*pi*pi*y[1]*y[2]*y[2]+2*pi*pi*pi*pi*y[0]*y[2]*y[3]+6*pi*pi*pi*pi*y[1]  *y[3]*y[3]);
dy[6]=16*c*y[7]/3-(pi*pi*lanmta2-fai)*pi*pi*y[2]-lanmta1*(3*pi*pi*pi*pi*y[2]*y[2]*y[2]/8+pi*pi*pi* y[0]*y[3]+3*pi*pi*pi*pi*y[2]*y[3]*y[3]+3*pi*pi*pi*pi*y[0]*y[0]*y[2]/8+pi*pi*pi*y[2]*y[1]+2*pi*pi*pi*pi*y[0]*y[1]*y[3]+pi*pi*pi*pi*y[1]*y[1]*y[2]);
dy[7]=-16*c*y[6]/3-4*(4*pi*pi*lanmta2-fai)*pi*pi*y[3]-lanmta1*(6*pi*pi*pi*pi*y[3]*y[3]*y[3]+pi*pi*pi*y[0]*y[2]+3*pi*pi*pi*pi*y[2]*y[2]*y[3]+pi*pi*pi*pi*y[0]*y[0]*y[3]+6*pi*pi*pi*pi*y[1]*y[1]*y[3]+2*pi*pi*pi*pi*y[0]*y[1]*y[2]);

无水1324

将你的方程编写成状态空间的形式保存为m函数然后就可以调用计算了。

function xxs=xxxs(t,x,flag,a0,a1)
a=a0+a1*sin(2*t);
% d2x/t2+a*x=0
xxs=[x(2);
    -a*x(1)];

%%%%%%%%%%上面既是mathieu方程状态方程的m函数

x0=[0.1;.1];
a0=6;a1=1;
[t,x]=ode45('xxxs',[0,100],x0,[],a0,a1)
plot(t,x)
%这是计算程序，你跟这一样编写程序就可以了。

还有可以在本版及matlab板块搜索到很多程序！

leadershuang

该如何把式子写成状态空间形式呢？

无水1324

就是简化为一阶的形式，然后把右边写在xxs=[；]里面

gghhjj

原帖由 leadershuang 于 2007-6-16 15:53 发表

                               
登录/注册后可看大图

该如何把式子写成状态空间形式呢？

x'=f(x)的形式就是