耦合的微分方程组不能求解出来

zhubaowei7

大家看看这个微分方程组 我在解的过程中遇到方程耦合的问题 解不出来了 大家帮帮忙

[ 本帖最后由 ChaChing 于 2010-6-28 18:47 编辑 ]

octopussheng

像你这种高自由度方程，可以用ode45/ode23/ode23tb试试

写成一阶方程形式，直接提交给matlab求解，论坛有很多这方面的例子，请使用搜索功能。

这种方程还是建议用fortran或者C来求解，毕竟matlab对高自由度的求解效率不是很高！

zhubaowei7

我就是用ode45做的  可是在选择状态变量时就不会了
x1=x,x2=x',x3=y,x4=y',x5=fan ,x6=fan'
                       x1'=x2
                       x2'=......x6'sin(fan)+fx/m
                       x3'=x4
                       x4'=......x6'cos(fan)+fy/m-g
                       x5'=x6
                       x6'=me/j+me[x2'sin(fan)-(x4'cos(fan))]+...
这样对吗  总感觉不太对啊
高手给看看  谢谢
其中  2  4  6式还是含有x6'  x2'   x4'  这该怎么办呢

[ 本帖最后由 ChaChing 于 2010-6-28 19:03 编辑 ]

octopussheng

强烈建议你把本版块中微分方程的部分搜索一下，这里再给你点提示

function dy=fff(t,y)

参数定义省略

dy=zeros(6,1)
dy1=y(2)
dy2=......
dy3=y(4)
dy4=......
dy5=y(6)
dy6=.....
dy=[dy1;dy2;dy3;dy4;dy5;dy6];

搜索微分方程，你会看到很多这方面的例子

无水1324

方程是不是可以做一定的化简

octopussheng

他这个比较复杂，不好化简！直接写成上面的形式还比较简单直观点，不过速度会慢！

zhubaowei7

方程两边都含最高阶的情况 比如
x''+2x'+8x=5y''+8xy+90
y''+4y'+9y=8yx''+7xy+70
参数定义怎么定呢?用dy=zeros(6,1)
dy1=y(2)
dy2=......
dy3=y(4)
dy4=......
dy5=y(6)
dy6=.....
dy=[dy1;dy2;dy3;dy4;dy5;dy6];
怎么做呢  谢谢

octopussheng

dy1=y(2)
dy2=-2*y(2)-8*y(1)+5*dy4+8*y(1)*y(3)+90
dy3=y(4)
dy4=-4*y(4)-9*y(3)+8*y(3)*dy2+7*y(1)*y(3)+70

求解下试试

zhubaowei7

function dudt=changode(t,u)
dudt=zeros(6,1); dudt(1)=u(4); dudt(2)=u(5); dudt(3)=u(6);
dudt(4)=-2*u(1)-u(2)-u(4)+3*sin(t);
dudt(5)=0.5*u(1)-u(2)+0.5*u(4)-0.5*u(6)+sin(t);
dudt(6)=(1/3)*(-4*u(3)-u(5))+2*sin(t);
dudt=dudt; %调用ode45求解该方程组，并画出结果图：
function changpingpin
u0=[-1;-3;-2;1;2;3]; tspan=[0 20]; [t,u]=ode45(@changode,tspan,u0);
plot(t,u(:,1),t,u(:,2),t,u(:,3)); legend('u(1)','u(2)','u(3)'); xlabel('t'); ylabel('u');

错误如下
Error in ==> changode at 3
dudt(1)=u(4);
怎么回事

[ 本帖最后由 ChaChing 于 2010-6-28 18:40 编辑 ]

花如月

怎么没有错误原因呢？楼主没帖全么？

eight

先看看各个置顶贴，请认真阅读

[ 本帖最后由 ChaChing 于 2010-6-28 18:38 编辑 ]

zhubaowei7

没看明白 什么意思

eight

原帖由 zhubaowei7 于 2007-8-3 13:45 发表

                               
登录/注册后可看大图

没看明白 什么意思

:@L 看这个：进版必读，建议先熟悉论坛界面及各种操作

xjzuo

这种简单错误自己看看基础书籍就能发现的, 希望以后能自己多思考:

%%% 将以下代码存为changpingpin.m,再运行即可
%%%%%%%%%%%%%%%%%%%%%%%%
function changpingpin
u0=[-1;-3;-2;1;2;3]; tspan=[0 20]; [t,u]=ode45(@changode,tspan,u0);
plot(t,u(:,1),t,u(:,2),t,u(:,3)); legend('u(1)','u(2)','u(3)'); xlabel('t'); ylabel('u');

function dudt=changode(t,u)
dudt=zeros(6,1); dudt(1)=u(4); dudt(2)=u(5); dudt(3)=u(6);
dudt(4)=-2*u(1)-u(2)-u(4)+3*sin(t);
dudt(5)=0.5*u(1)-u(2)+0.5*u(4)-0.5*u(6)+sin(t);
dudt(6)=(1/3)*(-4*u(3)-u(5))+2*sin(t);
%%%%%%%%%%%%%%%%%%%%%%%%

[ 本帖最后由 ChaChing 于 2010-6-28 18:43 编辑 ]

papayadong

先把最高阶的量统一放到方程左边，其余的放在右边，再编程就可以了！