如何实现牛顿-拉夫森来计算隐式龙格库塔的 k(i) 系数

Question

>我正在尝试实现具有fdm_2nd和高斯屠夫系数的 RK 隐式 2 阶对流扩散方程 （1D）：'u_t = -uu_x + nu .u_xx"。

我的目标是比较显式和隐式方案。显式 rk，适用于少量粘度。 显性谭的曲线向我们展示了一个非常好的冲击波。

我需要您的帮助才能正确实现 k（i） 系数的求解器。我不明白如何为所有 k（i） 实现牛顿方法。我需要为所有时空步骤实现它吗？还是及时？雅各布可能是错的，但我看不出在哪里。或者也许我在错误的方向上使用雅可比亚......

实际上，我的代码有效，但我认为某处是错误的......此外，隐式曲线不会从初始值移动。

这是我的功能：

function [t,u] = burgers(t0,U,N,dx) 
nu=0.01; %coefficient de viscosité
A=(diag(zeros(1,N))-diag(ones(1,N-1),1)+diag(ones(1,N-1),-1)) / (2*dx);
B=(-2*diag(ones(1,N))+diag(ones(1,N-1),1)+diag(ones(1,N-1),-1)) / (dx).^2;
t=t0;
u = - A * U.^2 + nu .* B * U;

雅可比：

function Jb = burJK(U,dx,i)
%Opérateurs
    a(1,1) = 1/4;
    a(1,2) = 1/4 - (3).^(1/2) / 6;
    a(2,1) = 1/4 + (3).^(1/2) / 6;
    a(2,2) = 1/4;
Jb(1,1) = a(1,1) .* (U(i+1,1) - U(i-1,1))/ (2*dx) - 1; 
Jb(1,2) = a(1,2) .* (U(i+1,1) - U(i-1,1))/ (2*dx);
Jb(2,1) = a(2,1) .* (U(i+1,2) - U(i-1,2))/ (2*dx);
Jb(2,2) = a(2,2) .* (U(i+1,2) - U(i-1,2))/ (2*dx) - 1;

这是我的牛顿代码：

iter = 1;
iter_max = 100;
k=zeros(2,N);
k(:,1)=[0.4;0.6];
[w_1,f1] =burgers(n + c(1) * dt,uu + dt * (a(1,:) * k(:,iter)),iter,dx);
[w_2,f2] =burgers(n + c(2) * dt,uu + dt * (a(2,:) * k(:,iter)),iter,dx);
f1 = -k(1,iter) + f1;
f2 = -k(1,iter) + f2;
f(:,1)=f1;
f(:,2)=f2;
df = burJK(f,dx,iter+1);
while iter<iter_max-1 % K_newton 
    delta = dff(iter,:)';
    k(:,iter+1) = k(:,iter) - delta;
    iter = iter+1;
    [w_1,f1] =burgers(n + c(1) * dt,uu + dt * (a(1,:) * k(:,iter+1)),N,dx);
    [w_2,f2] =burgers(n + c(2) * dt,uu + dt * (a(2,:) * k(:,iter+1)),N,dx);
    f1 = -k(1,iter+1) + f1;
    f2 = -k(1,iter+1) + f2;
    f(:,1)=f1;
    f(:,2)=f2;
    df = burJK(f,dx,iter);
    if iter>iter_max
        disp('#');
    else
        disp('ok');
    end
end

Answer 1

我对如何在 matlab 中实现这一点有点生疏，但我可以引导您完成一般步骤，希望这会有所帮助。 首先，我们可以考虑您正在求解的方程，以适合可以提出为的一般问题类别

du/dt = F(u), where F is a linear or nonlinear function 

对于 Runge Kutta 方案，您通常会将问题重新转换为这样的内容

k(i) = F(u+dt*a(i,i)*k(i)+ a(i,j)*k(j))

对于给定阶段。 现在是棘手的部分，您需要通过将k(1)堆叠到k(2)上来构建 1-D 矢量。所以向量的元素的前半部分是k(1)的，后半部分是k(2)的。使用这个新的组合向量，您可以更改F以便它分别在两个k上运行。这导致

K = FF(u+dt*a*K) where FF is F for the new double k vector, K

好了，现在我们可以实现牛顿方法了。 您将为每个时间步长执行此操作，直到您收敛到正确答案并同时在所有空间点上使用它。 你所做的是你猜测一个K并计算G(K,U) = K-FF(FF(u+dt*a*K)的雅可比。 只有当K处于正确的解决方案时，G(K,U)的值才应为零。换句话说，你牛顿的方法K，在寻找收敛时，你需要看到它在所有点上都收敛。 我会运行牛顿的方法直到max(abs(G(K,U)))< SolverTolerance.

抱歉，我在 matlab 实现方面无法提供更多帮助，但希望我帮助解释了如何实现牛顿方法。