龙格库塔法求解微分方程的例子，matlab

下面是一个用matlab求解微分方程的例子，使用的是龙格库塔法：

首先，我们定义微分方程和初始条件：

% 定义微分方程和初始条件
dydt = @(t,y) -2*t*y;
tspan = [0 2];
y0 = 1;

然后，我们使用ode45函数求解微分方程，作为比较的标准：

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

接下来，我们使用龙格库塔法求解微分方程：

% 使用龙格库塔法求解微分方程
h = 0.1; % 步长
t = tspan(1):h:tspan(2);
y = zeros(size(t));
y(1) = y0;

for i = 1:length(t)-1
    k1 = dydt(t(i), y(i));
    k2 = dydt(t(i)+h/2, y(i)+h*k1/2);
    k3 = dydt(t(i)+h/2, y(i)+h*k2/2);
    k4 = dydt(t(i)+h, y(i)+h*k3);
    y(i+1) = y(i) + h*(k1+2*k2+2*k3+k4)/6;
end

最后，我们画出两种方法求解的结果进行比较：

% 画出两种方法求解的结果进行比较
figure
plot(t, y, '-o', t, y, t, exp(-t.^2), '--')
legend('龙格库塔法', 'ode45', '真实解')
xlabel('t')
ylabel('y')

完整的代码如下：

% 定义微分方程和初始条件
dydt = @(t,y) -2*t*y;
tspan = [0 2];
y0 = 1;

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

% 使用龙格库塔法求解微分方程
h = 0.1; % 步长
t = tspan(1):h:tspan(2);
y = zeros(size(t));
y(1) = y0;

for i = 1:length(t)-1
    k1 = dydt(t(i), y(i));
    k2 = dydt(t(i)+h/2, y(i)+h*k1/2);
    k3 = dydt(t(i)+h/2, y(i)+h*k2/2);
    k4 = dydt(t(i)+h, y(i)+h*k3);
    y(i+1) = y(i) + h*(k1+2*k2+2*k3+k4)/6;
end

% 画出两种方法求解的结果进行比较
figure
plot(t, y, '-o', t, y, t, exp(-t.^2), '--')
legend('龙格库塔法', 'ode45', '真实解')
xlabel('t')
ylabel('y')