Matlab求解微分方程，入门实例教程

Matlab 求解常微分方程非常方便，不用你自己写龙格库塔，调用函数就行。本文从一阶到方程组，给完整例子。

一阶常微分方程，初值问题

求解形式：

$ \frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0 $

用 ode45，这是 Matlab 最常用的求解器，自适应步长四阶五阶龙格库塔，大部分情况用它就对了。

例子：求解 $\frac{dy}{dt} = -2y + sin(t)$，初始条件 $y(0) = 1$

步骤：

1. 定义微分方程函数

function dydt = myode(t, y)
    dydt = -2*y + sin(t);
end

也可以写匿名函数：

myode = @(t,y) -2*y + sin(t);

2. 设置求解区间和初始条件

tspan = [0 10]; % 从t=0算到t=10
y0 = 1; % 初始条件

3. 调用 ode45 求解

[t, y] = ode45(myode, tspan, y0);

4. 画图看结果

figure;
plot(t, y, 'LineWidth', 1.5);
xlabel('t');
ylabel('y');

就四步，出来结果了。

二阶微分方程怎么处理

高阶方程要转换成一阶方程组才能求解。比如：

$ y'' + 2y' + 2y = 0, \quad y(0)=0, y'(0)=1 $

做变量替换：

$ y_1 = y $
$ y_2 = y' $

得到方程组：

$ y_1' = y_2 $
$ y_2' = -2y_2 - 2y_1 $

写函数：

function dydt = second_order(t, y)
    dydt = zeros(2, 1);
    dydt(1) = y(2);
    dydt(2) = -2*y(2) - 2*y(1);
end

求解：

tspan = [0 10];
y0 = [0; 1]; % y1(0)=0, y2(0)=1
[t, y] = ode45(@second_order, tspan, y0);

figure;
plot(t, y(:,1), 'LineWidth', 1.5); % y(:,1)就是y
xlabel('t');
ylabel('y');

就这么转换，n阶方程转n个一阶方程，依此类推。

微分方程组

一阶方程组直接写，例子：洛伦兹吸引子

方程组：
$ x' = \sigma(y-x) $
$ y' = x(\rho - z) - y $
$ z' = xy - \beta z $

参数：σ=10, ρ=28, β=8/3

代码：

function dydt = lorenz(t, y, sigma, rho, beta)
    dydt = zeros(3, 1);
    dydt(1) = sigma * (y(2) - y(1));
    dydt(2) = y(1) * (rho - y(3)) - y(2);
    dydt(3) = y(1)*y(2) - beta*y(3);
end

调用，参数传进去：

sigma = 10;
rho = 28;
beta = 8/3;
tspan = [0 100];
y0 = [1; 1; 1];

[t, y] = ode45(@(t,y) lorenz(t,y,sigma,rho,beta), tspan, y0);

figure;
plot3(y(:,1), y(:,2), y(:,3));
xlabel('x');
ylabel('y');
zlabel('z');
title('洛伦兹吸引子');

跑出来就是那个著名蝴蝶图，混沌吸引子。

怎么选求解器

Matlab 好几个 ode 求解器，选哪个：

• ode45：大部分情况用这个，非刚性问题，自适应步长，精度不错，默认就选它
• ode23：精度要求低一点，比 ode45 快一点
• ode15s：刚性问题，迭代收敛慢的时候用这个
• ode23s：刚性问题，精度低一点快一点
• ode113：更高精度非刚性

不知道是不是刚性，先试 ode45，跑出来慢或者不收敛再换 ode15s。

边值问题怎么解

初值问题给起点，边值问题给两端边界条件，用 bvp4c 求解。

例子：$ y'' + |y| = 0 $，边界条件 $y(0)=0, y(4)=-2$

代码：

function ode_bvp_example
    solinit = bvpinit(linspace(0,4,10), [1 0]);
    sol = bvp4c(@ode, @bc, solinit);
    y = deval(sol, linspace(0,4));
    plot(linspace(0,4), y);

    function dydx = ode(x,y)
        dydx = [y(2); -abs(y(1))];
    end

    function res = bc(ya, yb)
        res = [ya(1); yb(1) + 2];
    end
end

就能解出来，边值问题可能多个解，初始猜测猜得好才能收敛到你要的解。

小结

求解微分方程步骤：

1. 高阶转一阶方程组
2. 写微分方程函数，输入t,y输出导数
3. 设置求解区间和初始条件
4. ode45 求解，画图

大部分非刚性问题 ode45 搞定，刚性换 ode15s，边值用 bvp4c，很简单。