声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 4499|回复: 29

[编程技巧] 龙格库塔法求助(万急)

  [复制链接]
发表于 2011-5-7 10:09 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
本帖最后由 煜宸0922 于 2011-5-7 10:14 编辑

function dxdt=fun(t,x)
k1=2;k2=3;k3=0.05;k4=20;k5=1;k6=1;k7=3;
dxdt=zeros(4,1);
dxdt(1)=x(3);
dxdt(2)=x(4);
dxdt(3)=(1-k1)*sin(k2*t+k4)-2*k3*x(3)+2*k3*x(4)-x(1)+x(2);
dxdt(4)=(k1*sin(k2*t+k4)+2*k3*x(3)-2*k3*(1+k5)*x(4)+x(1)-(1+k6)*x(2))/k7;
end
这是所要求数值解的方程组程序,不用ode45,请教各位怎样编写龙格库塔法解这个方程组的程序?万急,求救!!!
我编的程序老是出错,请各位大侠指正,不胜感激,程序如下:
function R=rk41(f,a,b,xa,xb,M)
h=0.01;
x=zeros(1,M+1);
t=a:h:b;
x(3)=xa;
x(4)=xb;
k11=x(3);
k21=x(4);
k31=feval(f,t,x(1),x(3));
k41=feval(f,t,x(2),x(4));
k12=x(3)+h/2*k31;
k22=x(4)+h/2*k41;
k32=feval(f,t+h/2,x(1)+h/2*k11,x(3)+h/2*k31);
k42=feval(f,t+h/2,x(2)+h/2*k21,x(4)+h/2*k41);
k13=x(3)+h/2*k32;
k23=x(4)+h/2*k42;
k33=feval(f,t+h/2,x(1)+h/2*k12,x(3)+h/2*k32);
k43=feval(f,t+h/2,x(2)+h/2*k22,x(4)+h/2*k42);
k14=x(3)+h*k33;
k24=x(4)+h*k43;
k34=feval(f,t+h,x(1)+h*k13,x(3)+h*k33);
k44=feval(f,t+h,x(2)+h*k23,x(4)+h*k43);
x(12)=x(1)+h/6*(k11+2*k12+2*k13+k14);
x(22)=x(2)+h/6*(k21+2*k22+2*k23+k24);
x(32)=x(3)+h/6*(k31+2*k32+2*k33+k34);
x(42)=x(4)+h/6*(k41+2*k42+2*k43+k44);
end
这程序是想只求解一个点,等正确了再循环,可还是错的,在下matlab刚入门,十分幼稚的错误在所难免,只求大侠们给予纠错帮助。临表涕零,不胜感激。

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2011-5-7 10:16 | 显示全部楼层
给个出错信息!
 楼主| 发表于 2011-5-7 10:17 | 显示全部楼层
回复 2 # meiyongyuandeze 的帖子

??? Error using ==> fun
Too many input arguments.

Error in ==> rk41 at 9
    k31=feval(f,t,x(1),x(3));
发表于 2011-5-7 10:34 | 显示全部楼层
回复 3 # 煜宸0922 的帖子

首先你定义的fun函数只有两个输入量(t,x),所以k31=feval(f,t,x(1),x(3));这句是错的,

评分

1

查看全部评分

 楼主| 发表于 2011-5-7 10:36 | 显示全部楼层
回复 4 # meiyongyuandeze 的帖子

呵呵,哪我该怎么改一下才行?k31那一句要用的是方程组的第三个方程,大侠指点下吧!
发表于 2011-5-7 11:00 | 显示全部楼层
本帖最后由 meiyongyuandeze 于 2011-5-7 11:02 编辑
  1. function dxdt=fun(t,x)
  2. k1=2;k2=3;k3=0.05;k4=20;k5=1;k6=1;k7=3;
  3. dxdt=zeros(4,1);
  4. dxdt(1)=x(3);
  5. dxdt(2)=x(4);
  6. dxdt(3)=(1-k1)*sin(k2*t+k4)-2*k3*x(3)+2*k3*x(4)-x(1)+x(2);
  7. dxdt(4)=(k1*sin(k2*t+k4)+2*k3*x(3)-2*k3*(1+k5)*x(4)+x(1)-(1+k6)*x(2))/k7;
复制代码
  1. clc
  2. clear
  3. k=0.001;
  4. T=0:0.001:10;
  5. X(:,1)=[0.005 0 0 0]';
  6. for j=1:length(T)
  7. t(j,1)=T(j);
  8. k1=k*feval('fun',T(j),X(:,j));
  9. k2=k*feval('fun',T(j)+k/2,X(:,j)+k1/2);
  10. k3=k*feval('fun',T(j)+k/2,X(:,j)+k2/2);
  11. k4=k*feval('fun',T(j)+k,X(:,j)+k3);
  12. X(:,(j+1))=X(:,j)+(k1+2*k2+2*k3+k4)/6;
  13. end
  14. plot(T,X(1,2:end))
复制代码
比较忙,实在是不想找错误,就帮你直接编写了,结果:
untitled.jpg

评分

1

查看全部评分

发表于 2011-5-7 11:01 | 显示全部楼层
回复 5 # 煜宸0922 的帖子

帮你写了下程序,你看下吧,图不知道是不是你想要的!
 楼主| 发表于 2011-5-7 21:59 | 显示全部楼层
回复 7 # meiyongyuandeze 的帖子

十分感谢大侠,我是要用这个画相图最后根据初值和其他参数画分岔图和混沌。我本来觉得,要每个方程都有四个斜率值,四个方程得十六个斜率值,而对每个方程算对应的斜率我就不会调了,十分感谢大侠百忙之中给予的帮助,以后又问题还得向大侠请教呢,希望大侠到时不吝赐教。
 楼主| 发表于 2011-5-7 22:07 | 显示全部楼层
回复 8 # 煜宸0922 的帖子

对了,大侠,你能解释下你画图的那个命令里面的东西吗?
 楼主| 发表于 2011-5-8 08:47 | 显示全部楼层
回复 7 # meiyongyuandeze 的帖子

对了,大侠,如果精度不够,要采用变步长的方法,程序该怎么调一下?
发表于 2011-5-8 09:19 | 显示全部楼层
回复 10 # 煜宸0922 的帖子

那改动还是很大的,你可以参考数值计算相关的书来编写!
 楼主| 发表于 2011-5-9 09:29 | 显示全部楼层
回复 11 # meiyongyuandeze 的帖子

麻烦大侠给推荐些书
 楼主| 发表于 2011-5-10 16:54 | 显示全部楼层
回复 2 # meiyongyuandeze 的帖子

大侠,还有个小问题,要是加入了边界条件,该怎么变一下程序?边界条件是:x(1)=b,x(1)=-b,x'(a1)=-Rx'(b1),这里面,x'(a1)表示碰撞前的速度,x'(b1)表示碰撞后的速度。请大侠指点一二。
发表于 2011-5-10 17:26 | 显示全部楼层
回复 13 # 煜宸0922 的帖子

应该是在function dxdt=fun(t,x)中加入判定语句!
 楼主| 发表于 2011-5-10 20:19 | 显示全部楼层
回复 14 # meiyongyuandeze 的帖子

求大侠详细解释!万急!!!!
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-13 14:12 , Processed in 0.083564 second(s), 28 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表