声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1972|回复: 4

[编程技巧] matlab求解微分方程,高手请进来帮帮忙

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

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

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

x
降阶以后的系统微分方程如下:
dx=[x(2);fx/M+md*cos(tao+fai);x(4);fy/M+md*sin(tao+fai)-G]; % x=x(1) x'=x(2) y=x(3) y'=x(4)
以下为系统里需要加入的非线性力
fx=(sqrt((x(1)-2*x(4))^2+(x(1)+2*x(2))^2)/(1-x(1)^2-x(3)^2))*(3*x(1)*V-G1*sin(a)-2*S*cos(a))+(RC/L^2)*(8*pi/sqrt(x(1)^2+x(3)^2))*(1-1/sqrt(1-x(1)^2-x(3)^2))*cos(a);
fy=(sqrt((x(1)-2*x(4))^2+(x(1)+2*x(2))^2)/(1-x(1)^2-x(3)^2))*(3*x(3)*V-G1*cos(a)-2*S*sin(a))+(RC/L^2)*(8*pi/sqrt(x(1)^2+x(3)^2))*(1-1/sqrt(1-x(1)^2-x(3)^2))*sin(a);
非线性力中所需的参数
V=(2+(X(3)*cos(a)-x(1)*sin(a))*G1)/(1-x(1)^2-x(3)^2);
S=(x(1)*cos(a)+x(3)*sin(a))/(1-(x(1)*cosa+x(3)*sin(a))^2);
G1=(2/sqrt(1-x(1)^2-x(3)^2))*(pai/2+atan((x(3)*cos(a)-x(1)*sin(a))/(1-x(1)^2-x(3)^2)));
a=atan((x(3)+2*x(2))/(x(1)-2*x(4)))-(pi/2)*sin((x(3)+2*x(2))/(x(1)-2*x(4)))-(pai/2)*sin(x(3)+2*x(2));
直接用ODE求解总会报错,
??? At compilation, "fx" was determined to be a variable and this
variable is uninitialized.  "fx" is also a function name and previous versions of MATLAB would have called the
function.
However, MATLAB 7 forbids the use of the same name in the same
context as both a function and a variable.

Error in ==> ymwd at 3
dx=[x(2);fx/M+md*cos(tao+fai);x(4);fy/M+md*sin(tao+fai)-G]; % x=x(1) x'=x(2) y=x(3) y'=x(4)

Error in ==> odearguments at 111
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
这类方程该怎么用matlab解?请高手帮帮忙啊 ,好几天了,一直纠结这个上面
回复
分享到:

使用道具 举报

发表于 2011-5-19 22:43 | 显示全部楼层
本帖最后由 meiyongyuandeze 于 2011-5-19 22:46 编辑

回复 1 # lhy 的帖子

首先你给的程序在书写上就有很多的错误:
1. V=(2+(X(3)*cos(a)-x(1)*sin(a))*G1)/(1-x(1)^2-x(3)^2);中的X(3)应该是x(3);
2. a=atan((x(3)+2*x(2))/(x(1)-2*x(4)))-(pi/2)*sin((x(3)+2*x(2))/(x(1)-2*x(4)))-(pai/2)*sin(x(3)+2*x(2)); 中的pai是什么啊,感觉应该是pi吧?
3.dx=[x(2);fx/M+md*cos(tao+fai);x(4);fy/M+md*sin(tao+fai)-G];中的tao,fai G不知道是怎么定义。
最好能给把参数的定义给全,你还是自己仔细看下的程序书写吧!

评分

1

查看全部评分

 楼主| 发表于 2011-5-20 10:35 | 显示全部楼层
谢谢指出错误,后面参数的定义没有粘上来,个人认为不在参数是否定义完全上,想请教个计算该类微分方程的方法。这是完整的M文件
function dx=ymwd(t,x)
global M md;
dx=[x(2);fx/M+md*cos(tao+fai);x(4);fy/M+md*sin(tao+fai)-G]; % x=x(1) x'=x(2) y=x(3) y'=x(4)
fx=(sqrt((x(1)-2*x(4))^2+(x(1)+2*x(2))^2)/(1-x(1)^2-x(3)^2))*(3*x(1)*V-G1*sin(a)-2*S*cos(a))+(RC/L^2)*(8*pi/sqrt(x(1)^2+x(3)^2))*(1-1/sqrt(1-x(1)^2-x(3)^2))*cos(a);
fy=(sqrt((x(1)-2*x(4))^2+(x(1)+2*x(2))^2)/(1-x(1)^2-x(3)^2))*(3*x(3)*V-G1*cos(a)-2*S*sin(a))+(RC/L^2)*(8*pi/sqrt(x(1)^2+x(3)^2))*(1-1/sqrt(1-x(1)^2-x(3)^2))*sin(a);
% 简写G1=G,a表示罗马数字
V=(2+(x(3)*cos(a)-x(1)*sin(a))*G1)/(1-x(1)^2-x(3)^2);
S=(x(1)*cos(a)+x(3)*sin(a))/(1-(x(1)*cosa+x(3)*sin(a))^2);
G1=(2/sqrt(1-x(1)^2-x(3)^2))*(pi/2+atan((x(3)*cos(a)-x(1)*sin(a))/(1-x(1)^2-x(3)^2)));
a=atan((x(3)+2*x(2))/(x(1)-2*x(4)))-(pi/2)*sin((x(3)+2*x(2))/(x(1)-2*x(4)))-(pi/2)*sin(x(3)+2*x(2));
% 无量纲化转换
tao=omg*t;md=e/C;qs=n*omg*R*L*(R/C)^2*(L/2*R)^2;
M=(omg^2*m*C)/qs;
G=g/C*omg^2;
X=x*C;Y=y*C;
m=2;    %kg
C=0.1; %mm
n=0.018;%Pa.s
R=25;%mm
L=10;%mm
g=9.8;%N/kg^2
发表于 2011-5-20 16:14 | 显示全部楼层
"fx" is also a function name and previous versions of MATLAB would have called the
function.
However, MATLAB 7 forbids the use of the same name in the same
context as both a function and a variable.


看这句 ,你似乎已经定义过fx为某个函数名了,改个称呼吧

评分

1

查看全部评分

发表于 2011-5-21 08:21 | 显示全部楼层
你这类微分方程也没有什么很特别的地方啊,感觉你M文件中的某些语句的位置可能不对!
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-11 06:55 , Processed in 0.074134 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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