1

fdm_2nd とガウス肉屋係数を使用して、RK 暗黙的な 2 次から対流拡散方程式 (1D) を実装しようとしています: 'u_t = -uu_x + nu .u_xx' 。

私の目標は、明示的スキームと暗黙的スキームを比較することです。粘度の少ない数値でうまく機能する明示的なrk。明示的なスキームの曲線は、非常に素晴らしい衝撃波を示しています。

k(i) 係数のソルバーを正しく実装するには、あなたの助けが必要です。すべての k(i) に対して newton メソッドを実装する方法がわかりません。すべての時空間ステップに実装する必要がありますか? またはジャストインタイム？ヤコビアンは間違っているかもしれませんが、どこにあるのかわかりません。それとも、ジャコビアンを間違った方向に使用しているのかもしれません...

実際、私のコードは機能しますが、どこか間違っていると思います...また、暗黙の曲線は初期値から移動しません。

ここに私の機能：

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

ルンゲ クッタ スキームの場合、通常、次のように問題を作り直します。

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

特定の段階のために。ここでトリッキーな部分が来k(1)ますk(2)。したがって、ベクトルの要素の前半は で、k(1)後半はk(2)です。この新しい結合ベクトルを使用して、 を変更して、2 つの を別々Fに操作できるようにします。kこれにより、

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

OK、これでニュートン法を実装できます。これを時間ステップごとに行い、正しい答えに収束するまでそれをすべての空間ポイントで同時に使用します。あなたがすることは、 a を推測しK、 のヤコビアンを計算することです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 に答える