MATLAB 龙格库塔法求解 m(t) = M0 * e^(-at)

使用 MATLAB 龙格库塔法求解 m(t) = M0 * e^(-at)

本教程将介绍如何使用 MATLAB 的龙格库塔法求解微分方程 m(t) = M0 * e^(-at)。

1. 定义函数 m(t) 及其导数 dm/dt

function dm = m_derivative(t, m, a, M0)
    dm = -a * m;
end

function m = m_func(t, a, M0)
    m = M0 * exp(-a * t);
end

其中，m_derivative 函数表示 m(t) 的导数 dm/dt，m_func 函数表示 m(t) 的值。

2. 设置初始条件和参数

t0 = 0; % 初始时间
tf = 10; % 结束时间
h = 0.1; % 步长
a = 0.1; % a 参数
M0 = 1; % 初始值

3. 使用龙格库塔法求解

t = t0:h:tf; % 时间步长
m = zeros(size(t)); % 存储 m(t) 的值
m(1) = M0; % 初始值
for i = 1:length(t)-1
    k1 = h * m_derivative(t(i), m(i), a, M0);
    k2 = h * m_derivative(t(i)+h/2, m(i)+k1/2, a, M0);
    k3 = h * m_derivative(t(i)+h/2, m(i)+k2/2, a, M0);
    k4 = h * m_derivative(t(i)+h, m(i)+k3, a, M0);
    m(i+1) = m(i) + 1/6 * (k1 + 2*k2 + 2*k3 + k4);
end

4. 绘制图像

plot(t, m);
xlabel('t');
ylabel('m(t)');
title('龙格库塔法求解 m(t)');

完整代码

function dm = m_derivative(t, m, a, M0)
    dm = -a * m;
end

function m = m_func(t, a, M0)
    m = M0 * exp(-a * t);
end

t0 = 0; % 初始时间
tf = 10; % 结束时间
h = 0.1; % 步长
a = 0.1; % a 参数
M0 = 1; % 初始值

t = t0:h:tf; % 时间步长
m = zeros(size(t)); % 存储 m(t) 的值
m(1) = M0; % 初始值
for i = 1:length(t)-1
    k1 = h * m_derivative(t(i), m(i), a, M0);
    k2 = h * m_derivative(t(i)+h/2, m(i)+k1/2, a, M0);
    k3 = h * m_derivative(t(i)+h/2, m(i)+k2/2, a, M0);
    k4 = h * m_derivative(t(i)+h, m(i)+k3, a, M0);
    m(i+1) = m(i) + 1/6 * (k1 + 2*k2 + 2*k3 + k4);
end

plot(t, m);
xlabel('t');
ylabel('m(t)');
title('龙格库塔法求解 m(t)');

通过以上步骤，你就可以使用 MATLAB 的龙格库塔法求解 m(t) = M0 * e^(-at) 微分方程并绘制其图像。