1

我正在尝试使用 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 = df\f(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
4

1 回答 1

0

我对如何在 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)。所以向量的前半部分是k(1)，后半部分是k(2)。使用这个新的组合向量，您可以更改F以便它k分别对两者进行操作。这导致

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

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

抱歉，我无法为 matlab 的实现提供更多帮助，但希望我能帮助解释如何实现牛顿法。

于 2016-03-08T16:25:42.497 回答