四阶龙格库塔matlab计算例题,四阶龙格库塔法matlab实现

以下是一个使用四阶龙格-库塔方法（RK4）求解常微分方程的matlab示例：

假设我们要求解以下常微分方程：

y' = -2ty^2，y(0) = 1

首先，将该微分方程转化为matlab函数：

function dydt = eqn(t,y)
dydt = -2*t*y^2;
end

然后，我们可以使用matlab内置的ode45函数生成一个解析解：

[t,y] = ode45(@eqn,[0 2],1);
plot(t,y)

接下来，我们使用RK4方法生成数值解。首先，我们定义时间步长和求解区间：

h = 0.1; % 时间步长
tspan = [0 2]; % 求解区间

然后，我们使用一个for循环来迭代计算数值解：

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

在每次迭代中，我们计算四个斜率（k1，k2，k3和k4），然后使用加权平均值（1/6 * k1 + 2/6 * k2 + 2/6 * k3 + 1/6 * k4）来更新y值。

最后，我们可以将数值解与解析解进行比较：

plot(t,y,'r',t,y_exact,'b')
legend('RK4','Exact')

完整的代码如下：

function dydt = eqn(t,y)
dydt = -2*t*y^2;
end

h = 0.1; % 时间步长
tspan = [0 2]; % 求解区间

% 解析解
[t_exact,y_exact] = ode45(@eqn,[0 2],1);

% RK4方法
t = tspan(1):h:tspan(2);
y = zeros(size(t));
y(1) = 1;
for i = 1:length(t)-1
    k1 = h * eqn(t(i),y(i));
    k2 = h * eqn(t(i)+h/2, y(i)+k1/2);
    k3 = h * eqn(t(i)+h/2, y(i)+k2/2);
    k4 = h * eqn(t(i)+h, y(i)+k3);
    y(i+1) = y(i) + 1/6*(k1+2*k2+2*k3+k4);
end

% 绘图
plot(t,y,'r',t_exact,y_exact,'b')
legend('RK4','Exact')