我正在尝试编写一个函数,该函数将使用 4 阶隐式 Runge-Kutta 方法 (IRK) 求解 ODES 系统,但我无法正确定义循环。这里我们定义 IRK
任何建议将不胜感激!
function [tout,yout] = IRK4Solver(f,t,y0)
t = t(:); % ensure that t is a column vector
N = length(t);
h = (t(end)-t(1))/(N-1); % calculate h by assuming t gridpoints are equidistant
d = length(y0); % defines the size of y0
y0 = y0(:); % ensures that y0 is a column vector
y = zeros(d,N);% allocate the output array y
y(:,1) = y0; % assign y0 to the first column of y
% The entries of the following tableau are provided in the lecture notes
c = [1/2-sqrt(3)/6;
1/2+sqrt(3)/6];
A = [1/4, 1/4-sqrt(3)/6;
1/4+sqrt(3)/6, 1/4];
b = [1/2;1/2];
%calculate the loop
for n=1:N
xi_1 = y(:,n)+h.*A(1,1).*f(t(n)+c(1).*h,xi_1)+h.*A(1,2)f(t(n)+c(2).*h,xi_2);
xi_2 = y(:,n)+h.*A(2,1).*f(t(n)+c(1).*h,xi_1)+h.*A(2,2)f(t(n)+c(2).*h,xi_2);
y(:,n+1) = y(:,n)+h.*b(1).*f(t(n)+c(1).*h,xi_1)+h.*b(2)f(t(n)+c(2).*h,xi_2);
end
tout = t;
yout = y;
进一步尝试
我已经包含了命令
fsolve
在我的 for 循环中。然而,porgram 仍然不会运行。for n=1:N
eq=@(xi) [xi(1:3)-(y(:,n)+h.*A(1,1).*f(t(n)+c(1).*h,xi(1:3))+h.*A(1,2)f(t(n)+c(2).*h,xi(1:3)));
xi(1:3)-(y(:,n)+h.*A(2,1).*f(t(n)+c(1).*h,xi(1:3))+h.*A(2,2)f(t(n)+c(2).*h,xi(1:3)))];
xistar=fsolve(eq,[1 1 1;1 1 1]);
y(:,n+1) = y(:,n)+h.*b(1).*f(t(n)+c(1).*h,xistar(1:3))+h.*b(2)f(t(n)+c(2).*h,xistar(1:3));
end
最佳答案
您需要包含命令 fsolve
在你的 for 循环中
关于matlab - 为隐式 Runge-Kutta 方法编写函数(四阶),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/61378206/