求用MATLAB求解二阶微分方程组的算法

[rush]

题在附件中，希望用ode来解并画出图~我是一个初学者，但现在急需用到，请帮帮忙！如果不是太清楚题目可以联系我的qq，512434053，谢谢

[yfl]

>> syms x5 xu t
>> m5=1290;mu=20;k5=10^5;co=1.5*10^5;kt=2*10^5;xr=0.5*sin(8*t);
>> [x5,xu]=dsolve('m5*D2x5+k5*(x5-xu)+co*(Dx5-Dxu)=0','mu*D2xu-k5*(x5-xu)-co*(Dx5-Dxu)+kt*(xu-xr)=0')

x5 =

xr+C1*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 1)*t)+C2*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 2)*t)+C3*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 3)*t)+C4*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 4)*t)



xu =

-(m5*co^2*C1*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 1)^2*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 1)*t)+m5*co^2*C2*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 2)^2*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 2)*t)+m5*co^2*C3*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 3)^2*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 3)*t)+m5*co^2*C4*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 4)^2*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 4)*t)-co^2*kt*xr+co*mu*m5*C1*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 1)^3*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 1)*t)+co*mu*m5*C2*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 2)^3*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 2)*t)+co*mu*m5*C3*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 3)^3*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 3)*t)+co*mu*m5*C4*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 4)^3*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 4)*t)+mu*co^2*C1*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 1)^2*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 1)*t)+mu*co^2*C2*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 2)^2*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 2)*t)+mu*co^2*C3*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 3)^2*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 3)*t)+mu*co^2*C4*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 4)^2*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 4)*t)-m5*mu*k5*C1*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 1)^2*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 1)*t)-m5*mu*k5*C2*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 2)^2*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 2)*t)-m5*mu*k5*C3*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 3)^2*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 3)*t)-m5*mu*k5*C4*RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 4)^2*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 4)*t)-k5^2*mu*xr-k5^2*mu*C1*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 1)*t)-k5^2*mu*C2*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 2)*t)-k5^2*mu*C3*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 3)*t)-k5^2*mu*C4*exp(RootOf(_Z^4*mu*m5+(co*mu+co*m5)*_Z^3+(m5*k5+m5*kt+mu*k5)*_Z^2+co*kt*_Z+k5*kt,index = 4)*t))/(co^2*kt+k5^2*mu)
代入数值后的解：
x5 =

1/2*sin(8*t)+C1*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 1)*t)+C2*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 2)*t)+C3*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 3)*t)+C4*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 4)*t)



xu =

-1473621/225010000*C1*RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 1)^2*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 1)*t)-1473621/225010000*C2*RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 2)^2*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 2)*t)-1473621/225010000*C3*RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 3)^2*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 3)*t)-1473621/225010000*C4*RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 4)^2*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 4)*t)+1/2*sin(8*t)-387/450020000*C1*RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 1)^3*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 1)*t)-387/450020000*C2*RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 2)^3*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 2)*t)-387/450020000*C3*RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 3)^3*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 3)*t)-387/450020000*C4*RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 4)^3*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 4)*t)+1/22501*C1*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 1)*t)+1/22501*C2*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 2)*t)+1/22501*C3*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 3)*t)+1/22501*C4*exp(RootOf(129*_Z^4+982500*_Z^3+1945000*_Z^2+150000000*_Z+100000000,index = 4)*t)
其中有些符号我也不认识，因为我才大二，不知对不对，你看看，如果是对的，那么请联系我，我再帮你画图，如不对，也不要见笑，因为我也是在学习吗，呵呵。我的QQ：602561167

[yfl]

补充化解：
>> pretty(x5)

1/2 sin(8 t) + C1 exp(RootOf(%1, index = 1) t)

+ C2 exp(RootOf(%1, index = 2) t) + C3 exp(RootOf(%1, index = 3) t)

+ C4 exp(RootOf(%1, index = 4) t)

4 3 2
%1 := 129 _Z + 982500 _Z + 1945000 _Z + 150000000 _Z + 100000000
>> pretty(xu)

1473621 2
- --------- C1 RootOf(%1, index = 1) exp(RootOf(%1, index = 1) t)
225010000

1473621 2
- --------- C2 RootOf(%1, index = 2) exp(RootOf(%1, index = 2) t)
225010000

1473621 2
- --------- C3 RootOf(%1, index = 3) exp(RootOf(%1, index = 3) t)
225010000

1473621 2
- --------- C4 RootOf(%1, index = 4) exp(RootOf(%1, index = 4) t)
225010000

+ 1/2 sin(8 t)

387 3
- --------- C1 RootOf(%1, index = 1) exp(RootOf(%1, index = 1) t)
450020000

387 3
- --------- C2 RootOf(%1, index = 2) exp(RootOf(%1, index = 2) t)
450020000

387 3
- --------- C3 RootOf(%1, index = 3) exp(RootOf(%1, index = 3) t)
450020000

387 3
- --------- C4 RootOf(%1, index = 4) exp(RootOf(%1, index = 4) t)
450020000

+ 1/22501 C1 exp(RootOf(%1, index = 1) t)

+ 1/22501 C2 exp(RootOf(%1, index = 2) t)

+ 1/22501 C3 exp(RootOf(%1, index = 3) t)

+ 1/22501 C4 exp(RootOf(%1, index = 4) t)

4 3 2
%1 := 129 _Z + 982500 _Z + 1945000 _Z + 150000000 _Z + 100000000

[rush]

谢谢你！但你这样求算的结果是不对的！这样的方程组最好用ode来解！我也是初学不怎么会！

[never]

小女拜求：用四阶龙格库塔解一阶或者二阶微分方程的程序例子，希望能得到相关资料 QQ：50713362 随时在线

[laosam280]

你怎么不给初值条件？
那我就按照系统默认的初值条件给你求了啊。
对于这样的方程组，首先你要转化成一阶方程组，才能用Runge-Kutta法进行迭代求解，不能直接用dsolve命令，那是求解一阶线性方程和方程组的。
我的求解程序如下，希望可以给到你一点帮助。
function dx=myfun01(t,x)
m5=1290; mu=20; c0=150000; k5=100000; kt=200000; xr=0.5*sin(8*t);
dx=[x(2);
-k5/m5*(x(1)-x(3))-c0/m5*(x(2)-x(4));
x(4);
k5/mu*(x(1)-x(3))+c0/mu*(x(2)-x(4))+kt*(x(3)-xr)];

%以上是方程组，保存到当前工作目录下，以下是求解代码。
[t,x]=ode45('myfun01',[时间范围],[初值]);
plot（看你自己要画什么类型的图，所有的数据都在向量x里）

[yfl]

楼上的你好，那个dx应该是个系数矩阵啊，你为啥直接那么写了？小弟正在学习他的解法，希望你能给个详细的解释，或者给出上面方程组的程序代码，谢了！！！