Matlab 使用四阶龙格库塔求解二阶隐式微分方程_ode45

---

一、ode45？

引用以下百度百科的内容：ode45，常微分方程的数值求解。MATLAB提供了求常微分方程数值解的函数。当难以求得微分方程的解析解时，可以求其数值解，Matlab中求微分方程数值解的函数有七个：ode45，ode23，ode113，ode15s，ode23s，ode23t，ode23tb 。
ode是Matlab专门用于解微分方程的功能函数。该求解器有变步长（variable-step）和定步长（fixed-step）两种类型。不同类型有着不同的求解器，其中ode45求解器属于变步长的一种，采用Runge-Kutta算法；其他采用相同算法的变步长求解器还有ode23。
ode45表示采用四阶-五阶Runge-Kutta算法，它用4阶方法提供候选解，5阶方法控制误差，是一种自适应步长（变步长）的常微分方程数值解法，其整体截断误差为(Δx)^5。解决的是Nonstiff(非刚性)常微分方程。
ode45是解决数值解问题的首选方法，若长时间没结果，应该就是刚性的，可换用ode15s试试。
这就是重点啦，ode45是使用了龙格库塔方法四阶五阶方法的一种求解常微分方程的方法，也是目前最常用的一种数值解计算方案。

二、案例解析

1.求解一阶常微分方程

考虑常微分方程为：
220sin($\omega (t)$) = Rs * y + Ls * dy + j * $\omega$ * y
除y和dy之外所有变量都已知，将方程进行变形，得到dy的表达式：
dy = ( 220sin($\omega (t)$) - Rs * y - j * $\omega$ * y ) / Ls
然后用odefun定义一个dy的表达式，格式如下：

odefun=@(t,y) (220*sin(wt) - Rs * y - j * w * y) / Ls; %定义dy表达式
tspan=[1 4]; %求解区间
y0=-2; %初值
[t,y]=ode45(odefun,tspan,y0);
plot(t,y) %作图
xlabel('t')
ylabel('dy')

2.求解二阶及高阶常微分方程

需要求解的高阶常微分方程：y’’ = -ty + e^ty’ + 3sin2t
求解的关键是将高阶转为一阶,odefun的书写.
F(y,y’,y’’…y(n-1),t)=0用变量替换,y1=y,y2=y’…注意odefun方程定义为行向量
dxdy=[y(1),y(2)…]

function Testode45
tspan=[3.9 4.0]; %求解区间
y0=[2 8]; %初值
[t,x]=ode45(@odefun,tspan,y0);
plot(t,x(:,1),'-o',t,x(:,2),'-*')
legend('y1','y2')
title('y'' ''=-t*y + e^t*y'' +3sin2t')
xlabel('t')
ylabel('y')
function y=odefun(t,x)
y=zeros(2,1); % 列向量
y(1)=x(2);
y(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t);
end
end

---