Question

二阶微分方程组的MATLAB解法

请告诉一种用ode45来解二阶微分方程组的方法，最好把算法写出来

Answer 1

一个二阶微分方程：
y''+y'+y=sin(t)

初始条件为y(0)=5,y'(0)=6。

过程：
先降阶为一阶微分方程组
y'=z
z'=-z-y+sin(t)

编制如下函数m文件
function dy=weifen(t,x)
dy=zeros(2,1);
%y=x(1)
%z=x(2)
dy(1)=x(2);
dy(2)=sin(t)-x(2)-x(1);

*******************
然后用ode45解方程
[t,y]=ode45(@weifen,[0 20],[5 6])
plot(t,y)就可以画出y和y'的图像

Answer 2

下面的程序是9自由度，二阶微分方程的rk解法：第一部分是定义函数：
定义函数
function ydot=f(t,Y,P)
m1=700
m2=150;%
m3=150;%对
m4=1700;%
m5=800;%
m6=600;%
i1=4;%
i2=3;%
i3=3;%

r1=0.225;%
r2=0.28;%
r3=0.28;%
k1=7000000;
k3=600000;
k4=600000;
kp=2000000;
ee=80000000000;
aa=0.000133;
nline=2;
lline=20;
kr1=nline*ee*aa/lline;
kr2=kr1;
khead=800000;
kr1s=kr1*khead/(kr1+khead);
kr2s=kr2*khead/(kr2+khead);
km11=k1+kr1+kr2;
km22=kr1+kr1s+k3;
km33=kr2+kr2s+k4;
km44=k4;
km55=k3+kp;
km66=kp;
km77=(kr1+kr2)*r1*r1*1.01;
km88=(kr1+kr1s)*r2*r2;
km99=(kr2+kr2s)*r3*r3;
km17=(kr1-kr2)*r1;
km27=-kr1*r1;
km37=kr2*r1;
km18=kr1*r2;
km28=(-kr1+kr1s)*r2;
km78=kr1*r1*r2;
km19=-kr2*r3;
km39=(kr2-kr2s)*r3;
km79=kr2*r1*r3;
M=[m1 0 0 0 0 0 0 0 0;0 m2 0 0 0 0 0 0 0;0 0 m3 0 0 0 0 0 0;0 0 0 m4 0 0 0 0 0;0 0 0 0 m5 0 0 0 0;0 0 0 0 0 m6 0 0 0;0 0 0 0 0 0 i1 0 0;0 0 0 0 0 0 0 i2 0;0 0 0 0 0 0 0 0 i3];
K=[km11 -kr1 -kr2 0 0 0 km17 km18 km19;-kr1 km22 0 0 -k3 0 km27 km28 0;-kr2 0 km33 -k4 0 0 km37 0 km39;0 0 -k4 km44 0 0 0 0 0;0 -k3 0 0 km55 -kp 0 0 0;0 0 0 0 -kp kp 0 0 0;km17 km27 km37 0 0 0 km77 km78 km79;km18 km28 0 0 0 0 km78 km88 0;km19 0 km39 0 0 0 km79 0 km99];
C=zeros(9,9);
%C=0.01*[km11 -kr1 -kr2 0 0 0 km17 km18 km19;-kr1 km22 0 0 -k3 0 km27 km28 0;-kr2 0 km33 -k4 0 0 km37 0 km39;0 0 -k4 km44 0 0 0 0 0;0 -k3 0 0 km55 -kp 0 0 0;0 0 0 0 -kp kp 0 0 0;km17 km27 km37 0 0 0 km77 km78 km79;km18 km28 0 0 0 0 km78 km88 0;km19 0 km39 0 0 0 km79 0 km99];
I=eye(9,9);
A=[zeros(9,9),I;-inv(M)*K,-inv(M)*C];
ydot=A*Y+P;