声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2566|回复: 7

[编程技巧] 大家帮我看看这个ODE45解微分方程程序哪里出错

[复制链接]
发表于 2012-9-18 21:13 | 显示全部楼层 |阅读模式

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

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

x

M文件如下
function dx=f2(t,x);
m1=1.8491 m2=1.3950 me=0.7951
a1=0.03   wg=0.46  we=0.02 Fm=0.04 Fin=0.2
if x(1)>1
    f_x(1)=x(1)-1;
else if x(1)<-1
        f_x(1)=x(1)+1;
    else  f_x(1)=0;
    end
dx(1)=x(2);
dx(2)=-2*m1/m2*a1*x(2)-me/m2*f_x(1)+Fm*(1-0.1*sin(we*t))+Fin*wg^2*sin(wg*t);
  %dx(2)=-2*me/m2*a1*x(2)-me/m2*(1+a2*sin(wg*t))*f_x(1)+Ft
  %Ft=Fm*(1-0.1*sin(we*t))+Fin*wg^2*sin(wg*t)
  dx=[dx(1);dx(2)];
运行如下
clear all;
clc;
H=[0 100];x0=[0,0];
[t,x]=ode45('f2',H,x0);
plot(x(:,1),x(:,2),'m-');
错误如下
??? Error: File: f2.m Line: 2 Column: 11
Unexpected MATLAB expression.

Error in ==> funfun\private\odearguments at 80
    if (nargin(ode) == 2)           

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> Untitled6 at 4
[t,x]=ode45('f2',H,x0);
着急 自己弄两天了 还是弄不出来 请高手帮忙
回复
分享到:

使用道具 举报

发表于 2012-10-5 23:33 | 显示全部楼层
正好最近也在解微分方程,试着弄了一下。其实没多大问题,主要是matlab语句的写法上。

1、function 定义时,不需要分号;
2、参数赋值时,要么之间加分号,要么另起一行;
3、条件语句的写法问题,好像是:if...elseif...else...end,elseif连写
最后,加上运行的结果:
f2.png

点评

建议补齐代码分享  发表于 2012-10-6 19:40

评分

1

查看全部评分

发表于 2012-10-6 12:14 | 显示全部楼层
1stOpt求解:

  1. Constant m1=1.8491,m2=1.3950,me=0.7951,
  2.          a1=0.03,wg=0.46,we=0.02,Fm=0.04,Fin=0.2;
  3. ConstStr f_x=if(x1>1,x1-1,if(x1<-1, x1+1,0));
  4. Variable t=[0,100],x1=0,x2=0;
  5. Plot x1[x],x2;
  6. ODEFunction x1'=x2;
  7.             x2'=-2*m1/m2*a1*x2-me/m2*f_x+Fm*(1-0.1*sin(we*t))+Fin*wg^2*sin(wg*t);
复制代码
jv1.jpg

评分

1

查看全部评分

发表于 2012-10-12 22:30 | 显示全部楼层
%%%下面是微分方程文件
function dx=f2(t,x)
m1=1.8491;
m2=1.3950;
me=0.7951;
a1=0.03;
wg=0.46;
we=0.02;
Fm=0.04;
Fin=0.2;
if x(1)>1
    f_x(1)=x(1)-1;
elseif x(1)<-1
    f_x(1)=x(1)+1;
else
    f_x(1)=0;
end
dx(1)=x(2);
dx(2)=-2*m1/m2*a1*x(2)-me/m2*f_x(1)+Fm*(1-0.1*sin(we*t))+Fin*wg^2*sin(wg*t);
  %dx(2)=-2*me/m2*a1*x(2)-me/m2*(1+a2*sin(wg*t))*f_x(1)+Ft
  %Ft=Fm*(1-0.1*sin(we*t))+Fin*wg^2*sin(wg*t)
  dx=[dx(1);dx(2)];

%%下面是主函数
clear all;
clc;
H=0.1:0.1:100;
x0=[0,0];
[t,x]=ode45('f2',H,x0);
plot(x(:,1),x(:,2),'m-');


顺便问一个技术性问题:论坛里的程序怎么贴上去的?可以显示行数,又可以直接点复制。


 楼主| 发表于 2012-10-13 09:47 | 显示全部楼层
发表于 2012-10-18 11:38 | 显示全部楼层
jixiantianzi 发表于 2012-10-12 22:30
%%%下面是微分方程文件
function dx=f2(t,x)
m1=1.8491;

编辑区选择插入代码文字就行,如图红圈标识位置
code2.png
code1.png

评分

1

查看全部评分

发表于 2012-10-18 20:52 | 显示全部楼层
happy 发表于 2012-10-18 11:38
编辑区选择插入代码文字就行,如图红圈标识位置

学到了,哈哈!
试验一把

  1. %%%下面是微分方程文件
  2. function dx=f2(t,x)
  3. m1=1.8491;
  4. m2=1.3950;
  5. me=0.7951;
  6. a1=0.03;
  7. wg=0.46;
  8. we=0.02;
  9. Fm=0.04;
  10. Fin=0.2;
  11. if x(1)>1
  12.     f_x(1)=x(1)-1;
  13. elseif x(1)<-1
  14.     f_x(1)=x(1)+1;
  15. else
  16.     f_x(1)=0;
  17. end
  18. dx(1)=x(2);
  19. dx(2)=-2*m1/m2*a1*x(2)-me/m2*f_x(1)+Fm*(1-0.1*sin(we*t))+Fin*wg^2*sin(wg*t);
  20.   %dx(2)=-2*me/m2*a1*x(2)-me/m2*(1+a2*sin(wg*t))*f_x(1)+Ft
  21.   %Ft=Fm*(1-0.1*sin(we*t))+Fin*wg^2*sin(wg*t)
  22.   dx=[dx(1);dx(2)];

  23. %%下面是主函数
  24. clear all;
  25. clc;
  26. H=0.1:0.1:100;
  27. x0=[0,0];
  28. [t,x]=ode45('f2',H,x0);
  29. plot(x(:,1),x(:,2),'m-');
复制代码


发表于 2013-6-8 09:23 | 显示全部楼层
哈哈 我也来试验一把 成功了
  1. a=1;
  2. b=2;
复制代码
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-14 14:01 , Processed in 0.084937 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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