四阶龙格库塔matlab
时间: 2023-06-27 21:06:04 浏览: 706
以下是四阶龙格-库塔法(RK4)的MATLAB代码示例:
```matlab
function [t, y] = rk4(f, tspan, y0, h)
% RK4 Classical 4th order Runge-Kutta method for
% solving the initial value problem y' = f(t,y), y(tspan(1)) = y0.
% USAGE: [t,y] = rk4(f, tspan, y0, h)
% INPUT:
% f - function handle of f(t,y)
% tspan - [t0 tf]
% y0 - initial condition; y(t0) = y0
% h - step size
% OUTPUT:
% t - column vector of times at which y was computed
% y - column vector of approximate solution corresponding to t
%
% Example:
% f = @(t,y) -y/t + 4*t^3;
% tspan = [1 2]; y0 = 2; h = 0.1;
% [t,y] = rk4(f, tspan, y0, h);
% plot(t,y,'o-',t,2*t.^2+2./t,'r')
% legend('RK4','Exact',0)
t0 = tspan(1); tf = tspan(2);
if h <= 0, error('Step size must be positive'); end
if tf < t0, error('The final time must be greater than the initial time'); end
t = (t0:h:tf)';
n = length(t);
y = zeros(n, 1);
y(1) = y0;
for i = 1:n-1
k1 = h*f(t(i), y(i));
k2 = h*f(t(i) + h/2, y(i) + k1/2);
k3 = h*f(t(i) + h/2, y(i) + k2/2);
k4 = h*f(t(i) + h, y(i) + k3);
y(i+1) = y(i) + (k1 + 2*k2 + 2*k3 + k4)/6;
end
end
```
使用示例:
```matlab
% 定义一个函数
f = @(t, y) -y/t + 4*t^3;
% 设定初始条件和时间范围
tspan = [1 2];
y0 = 2;
h = 0.1;
% 调用RK4函数求解
[t, y] = rk4(f, tspan, y0, h);
% 绘制结果
plot(t, y, 'o-', t, 2*t.^2+2./t, 'r');
legend('RK4', 'Exact', 0);
```
在上述示例中,我们定义了一个函数 $f(t,y)=-\frac{y}{t}+4t^3$,并通过 `rk4` 函数求解了其在 $t\in[1,2]$ 时刻的解 $y(t)$。我们同时也绘制了精确解 $y(t)=2t^2+\frac{2}{t}$ 以方便比较。