如何用Matlab编写龙格-库塔求解二阶微分方程

下面是Matlab代码实现龙格-库塔求解二阶微分方程：

function [t, y] = RK2_solver(f, y0, tspan, h)
% f: the right-hand side function in the differential equation y'' = f(t, y, y')
% y0: initial condition, y0(1) is the initial value of y, y0(2) is the initial value of y'
% tspan: [t0, tf], the time span of the simulation
% h: step size

% Initialize the output arrays
t = tspan(1):h:tspan(2);
y = zeros(2, length(t));
y(:, 1) = y0;

% Define the Runge-Kutta coefficients
a = [0 1/2];
b = [0 1];
c = [1/2 1/2];

% Iterate over time steps
for i = 1:length(t)-1
    % Calculate the intermediate variables
    k1 = h * [y(2, i); f(t(i), y(1, i), y(2, i))];
    k2 = h * [y(2, i)+b(2)*k1(2); f(t(i)+c(2)*h, y(1, i)+a(2)*k1(1), y(2, i)+b(2)*k1(2))];
    % Update the solution
    y(:, i+1) = y(:, i) + (a(1)*k1 + a(2)*k2);
end
end

其中，`f` 是二阶微分方程中的右侧函数，`y0` 是初始条件（包括 $y(0)$ 和 $y'(0)$），`tspan` 是时间范围，`h` 是时间步长。函数返回一个时间数组 `t` 和一个对应的状态数组 `y`，其中 `y(1,:)` 是 $y$ 的值，`y(2,:)` 是 $y'$ 的值。

例如，假设要求解下面的二阶微分方程：

$$y''(t) + 2y'(t) + 2y(t) = \sin(t)$$

并且初始条件为 $y(0) = 1, y'(0) = 0$，时间范围为 $t \in [0, 10]$，时间步长为 $h=0.1$，则可以用下面的代码求解：

% Define the right-hand side function
f = @(t, y, yp) sin(t) - 2*yp - 2*y;

% Define the initial conditions and time span
y0 = [1; 0];
tspan = [0, 10];
h = 0.1;

% Solve the differential equation
[t, y] = RK2_solver(f, y0, tspan, h);

% Plot the solution
plot(t, y(1,:))
xlabel('t')
ylabel('y')
title('Solution of y'''' + 2y'' + 2y = sin(t)')

运行后，可以得到下面的图像：


这个图像显示了在 $t\in[0,10]$ 范围内的 $y$ 的解。