大家帮我看看这个ODE45解微分方程程序哪里出错

科比vs流川枫

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);
着急 自己弄两天了 还是弄不出来 请高手帮忙

jixiantianzi

正好最近也在解微分方程，试着弄了一下。其实没多大问题，主要是matlab语句的写法上。

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

dingd

1stOpt求解：

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

复制代码

jixiantianzi

%%%下面是微分方程文件
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-');


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

科比vs流川枫

jixiantianzi 发表于 2012-10-12 22:30

                               
登录/注册后可看大图

%%%下面是微分方程文件
function dx=f2(t,x)
m1=1.8491;

我也不知道 同求 的确我是因为 else 和 if 没有连在一起写 而程序出错

happy

jixiantianzi 发表于 2012-10-12 22:30

                               
登录/注册后可看大图

%%%下面是微分方程文件
function dx=f2(t,x)
m1=1.8491;

编辑区选择插入代码文字就行，如图红圈标识位置

jixiantianzi

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. %%下面是主函数
    
25. clear all;
    
26. clc;
    
27. H=0.1:0.1:100;
    
28. x0=[0,0];
    
29. [t,x]=ode45('f2',H,x0);
    
30. plot(x(:,1),x(:,2),'m-');

复制代码

tiandong0213

哈哈 我也来试验一把 成功了

1. a=1;
   
2. b=2;
   

复制代码