insects 发表于 2007-5-16 13:23

求教龙格库塔法解微分方程

有没有谁会编个二阶用龙格库塔法解非线性微分方程的程序给我看看，用matlab语言的，方程是
          xdot(2)=(-(1/tau1+K*tau2*(cosx(1))/tau1)*x(2)-K*(sinx(1))/tau1+Deltaomegao/tau1);
         xdot(1)=x(2);
其中xdot(1)就是x1求一阶导,xdot(2)就是x2求一阶导，Deltaomegao，tau1，tau2，K是要给定初值的常量
我把其中的cosx(1)，sinx(1)用泰勒公式展开了前三项然后用ode45编了程序，但是老师说这样不精确，也不专业，要我用龙格库塔法，不会，所以请大家帮帮忙，看看怎么改。
本人的程序如下：
function xdot=TwoPLL2(t,x,dummy,tau1,tau2,K,Deltaomegao)
         xdot=zeros(2,1);
         xdot(2)=(-(1/tau1+K*tau2*(1-x(1)^2/2+x(1)^4/24)/tau1)*x(2)-K*(x(1)-x(1)^3/6+x(1)^5/120)/tau1+Deltaomegao/tau1);
         xdot(1)=x(2);
主程序
tau1=5;tau2=1;K=12;Deltaomegao=;
tspan=linspace(0,15,1500);
options=odeset('Abstol',);
for i=1:4
   =ode45('TwoPLL2',tspan,[-pi,Deltaomegao(i)]',options,tau1,tau2,K,Deltaomegao(i));
   figure(1);plot(t,x(:,1));       hold on; title('Phase difference as a function of Time');ylabel('θe(t)/rad');xlabel('t/s');
   figure(2);plot(t,x(:,2));       hold on; title('Frequency difference as a function of Time');ylabel('(dθe(t)/dt)/(rad?sˉ1)');xlabel('t/s');
   figure(3);plot(x(:,1),x(:,2));hold on; title('Phase portrait');ylabel('(dθe(t)/dt)/(rad?sˉ1)');xlabel('/rad');
end:'(

insects 发表于 2007-5-16 13:26

上面不小心打错了一点内容，对不起了

:@L 不好意思，打错了，是四阶龙格库塔法的matlab程序求教，

咕噜噜 发表于 2007-5-16 14:42

回复 #2 insects 的帖子

你确定你上述给出的是四阶微分方程吗？不是二阶？你的m函数里面定义显示只有二阶
四阶龙格库塔法的matlab程序很多matlab书上都有，另ode45其实就是以龙格库塔法为基础的编程

gghhjj 发表于 2007-5-17 07:53

第一个问题
ode45并非4阶龙格库塔法，是四五阶变步长龙格库塔法
如果要用平常四阶龙格库塔法求解，请到matlab版搜索，有专门的matlab程序

当然这里用ode45求解应该是没有问题的，一般来说ode45的求解要比4阶龙格库塔法求解好

[ 本帖最后由 xinyuxf 于 2007-5-18 09:27 编辑 ]

gghhjj 发表于 2007-5-17 07:55

原帖由 咕噜噜 于 2007-5-16 14:42 发表 http://www.chinavib.com/forum/images/common/back.gif
你确定你上述给出的是四阶微分方程吗？不是二阶？你的m函数里面定义显示只有二阶
四阶龙格库塔法的matlab程序很多matlab书上都有，另ode45其实就是以龙格库塔法为基础的编程

搂主没说四阶微分方程，说的就是二阶
龙格库塔法的阶数和方程的价数无关

gghhjj 发表于 2007-5-17 07:57

第二个问题
类似的微分方程程序我在本版前面写过好多个，自己搜索以下看看怎么回事吧

另外doc ode45看帮助，里面有详细的例子

yufeng 发表于 2007-5-21 13:31

联系me：qq174456366，
页: [1]
查看完整版本: 求教龙格库塔法解微分方程