为隐式 Runge-Kutta 方法编写函数（四阶）

Posted 2023-02-19

【中文标题】为隐式 Runge-Kutta 方法编写函数（四阶）

【英文标题】：Writing a function for the Implicit Runge-Kutta method (order four)

【发布时间】：2020-08-06 05:54:18

【问题描述】：

我正在尝试编写一个函数，该函数将使用 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;

进一步的尝试

我在 for 循环中包含了命令 fsolve。然而，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

【问题讨论】：

您是否尝试过从小处着手，实现隐式欧拉法或隐式梯形法？你对 Jacobi 矩阵了解多少？

@LutzLehmann 我没有尝试过更简单的隐式数值方法，但我会试一试。雅可比矩阵有什么帮助？

你有什么问题？看起来不错。我建议您使用内部循环进行求和，但如果 $\nu$=2 则不值得花时间。

@AlessandroTrigilio 我正在考虑使用循环进行求和，但明确地编写 $\xi_1$ 和 $\xi_2$ 似乎更容易。但是，这是一个非线性方程，那么如何计算 $y_n+1$？

我在对***.com/q/53910201/3088138 和***.com/a/61223515/3088138 的回答中做了类似的事情。这并不完全相同，但结构应该大致相同。要获得比在步骤开始时使用数据中的常数或线性函数更好的初始点，请使用前一段的插值函数来推断初始猜测。这应该会在执行时间上产生可测量的差异。

【参考方案1】：

我使用了四阶 Hammer 和 Hollingsworth (1955) 的 2 阶段方案，该方案在“R. May 和 J. Noye，“常微分方程的数值解：初值问题”中获得，North-holl . Math. Stud.，第 83 卷，第 1-94 页，1984 年 1 月"。

我的代码是：

function [tout,yout] = IRK4Solver(f,t,y0) 

t = t(:); % ensure that t is a column vector
N = length(t); 
h = t(2)-t(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

k1 = 0;
k2 = 0;

%calculate the loop
for n=2:N    
    k1 = f(t(n-1)+(0.5+sqrt(3)/6)*h,y(:,n-1)+1/4*k1+(1/4+sqrt(3)/6)*h*k2);
    k2 = f(t(n-1)+(0.5-sqrt(3)/6)*h,y(:,n-1)+1/4*k2+(1/4-sqrt(3)/6)*h*k1);
    y(:,n) = y(:,n-1)+(h/2)*(k1+k2);
end

tout = t;
yout = y;

我希望它是您正在寻找的。​​p>

【参考方案2】：

您需要在 for 循环中包含命令 fsolve