龙格库塔法求解三阶微分方程Matlab代码例子

以下是一个使用龙格库塔法求解三阶微分方程的Matlab代码示例：

function [T,Y] = thirdOrderODE_RK3(f,tspan,y0,h)
% f为方程右侧函数
% tspan为时间区间
% y0为初值
% h为步长
t0 = tspan(1);
tf = tspan(2);
T = t0:h:tf; % 时间序列
n = length(T);
Y = zeros(n,length(y0)); % 初始化Y
Y(1,:) = y0; % 设定初值
for i = 1:(n-1)
    k1 = h*feval(f,T(i),Y(i,:)).'; % 计算k1
    k2 = h*feval(f,T(i)+h/2,Y(i,:)+k1/2).'; % 计算k2
    k3 = h*feval(f,T(i)+h,Y(i,:)-k1+2*k2).'; % 计算k3
    Y(i+1,:) = Y(i,:) + (k1+4*k2+k3)/6; % 计算下一个点的值
end
end

例如，如果我们要求解如下的三阶微分方程：

y''' - 3y'' + 3y' - y = 0
y(0) = 1, y'(0) = 0, y''(0) = 0

可以使用以下代码调用上面的函数：

f = @(t,y) [y(2); y(3); y(1)+3*y(2)-3*y(3)]; % 定义右侧函数
tspan = [0 10]; % 时间区间
y0 = [1; 0; 0]; % 初值
h = 0.01; % 步长
[T,Y] = thirdOrderODE_RK3(f,tspan,y0,h); % 调用函数求解
plot(T,Y(:,1)); % 绘制y随时间的变化图像

这将绘制出y随时间的变化图像。