Matlab ODE Introduction

ODE是matlab中用于求解微分方程的一个功能函数，其中 ode45是一种求解器solver，采用四阶五阶龙格库塔法 Runge_Kutta。

ode45函数需要传入三个参数：ode函数，时间范围，初始条件。返回两个值：时间点和对应时刻的函数的值，也就是原函数的数值解。

[t, y] = ode45(ode_function, tspan, y0)

How to use

一阶

最简单的一阶线性常微分方程形式为：

$$ \frac{dy}{dx} + P(x)y = Q(x) $$

这种情况我们只需要把$\frac {dy} {dx}$移项到一侧，其余内容在另一侧即可。

例如对于 $y\frac{dy}{dx} + (x^2 + 1)y = x$，我们需要先通过整理得到$\frac{dy}{dx} = \frac{x - (x^2 + 1)y}{y}$，然后再在matlab中定义函数：

function dydx = myODE(x, y)
    dydx = (x - (x^2 + 1)*y) / y;
end
% 补充完整：
% 定义时间范围
tspan = [0 10];

% 定义初始条件
y0 = 1;

% 使用ode45求解微分方程
[t, y] = ode45(@myODE, tspan, y0);

%绘图：
plot(t, y);
xlabel('x');
ylabel('y');
title('Solution of the ODE');

初始值可以是 tspan中起始点对应的 y的值。

注意：ode方程一定要有两个传入参数，即使分离微分项之后，右侧没有 x或 y，传入参数也必须是两个。这是 ode45，ode15s等求解器的要求。此外，ode45的形参初始条件，在 ode方程中变成 y，可以在 ode方程中调用。

此外：别忘了把函数的定义移动到文件的最后，这是脚本文件的格式要求。

二阶

对于二阶，我们需要引入一个新的变量进行表示。同时，需要改变ode函数的返回值，让其变成一个向量：

function dydx = myODE(x, y)
    dydx = zeros(2, 1) % 生成一个2行1列的矩阵（向量）
end

两个元素分别代表了 y对 x的一阶导数和二阶导数。

同时，我们也需要有两个初始值：

y0 = [1; 0] % y0是一个2行1列的矩阵（向量）

两个元素分别代表了 y本身和 y对 x的一阶导数。

那么对于这个方程$\frac{d^2y}{dx^2} + 2\frac{dy}{dx} + 2y = 0$，我们先换元，并且移项整理，得到：

$\frac{dy}{dx} = z$

$\frac{dz}{dx} = -2z - 2y$

然后可以写出需要求解的ode方程形式：

function dydx = myODE(x, y)
    dydx = zeros(2, 1);
    dydx(1) = y(2);
    dydx(2) = -2 * y(2) - 2 * y(1);
end

再调用 ode45函数，可以得到一个完美的 y值。

此时的 y是一个矩阵，有两列，每一列有很多行。第一列表示 y的值在不同的 x情况下；第二列表示 y对 x的一阶导数值在不同的 x情况下。在绘制图像时，需要写清楚用 y矩阵的哪一列。

plot(x, y(:, 1)); % 使用y的第一列。matlab索引从1开始