我有一个Matlab函数,用于对一阶ODE进行Runge-Kutta4k近似,我想将其调整为适用于二阶ODE。有人能帮我开始吗?这是我的第一个订单RK4K:
function [y,t,h] = rungekutta4kh(F,y0,a,b,n)
% Euler ODE solver
t = linspace(a,b,n);
h = t(2)-t(1);
y = zeros(n,1);
y(1) = y0;
for i=2:n
s1 = feval(F,t(i-1),y(i-1));
s2 = feval(F,t(i-1)+h/2,y(i-1)+ (h/2)*s1);
s3 = feval(F,t(i-1)+h/2,y(i-1)+ (h/2)*s2);
s4 = feval(F,t(i-1)+h,y(i-1)+ h*s3);
y(i) = y(i-1) + ...
(h/6)*( s1 + 2*s2 + 2*s3 + s4 );
end
将其转换为一阶系统,这意味着函数F
变为向量值,并且应该在y
列表的构造中包括y0
的形状。
y''=f(x,y,y')
转换为
y0'=y1
y1'=f(x,y0,y1)
以便(使用matlab索引)
F(x,y) = [ y(2), f(x,y(1),y(2)) ]