как реализовать newton-raphson для вычисления коэффициентов k (i) неявной руны kutta?

Question

Я пытаюсь реализовать RK неявное уравнение для конвекции-диффузии (1D) с коэффициентами мяча fdm_2nd и gauss: 'u_t = -uu_x + nu .u_xx'.

Моя цель - сравнить схему экспликации против схемы implcit. Явный rk, который хорошо работает с небольшим количеством вязкости. Кривая явной схемы показывает нам очень приятную ударную волну.

Мне нужна ваша помощь для правильной реализации решателя коэффициента k (i). Я не вижу, как реализовать метод newton для всех k (i). Мне нужно реализовать его для всех шагов времени? или как раз вовремя? Якобиан, возможно, ошибается, но я не вижу, где. Или, может быть, я использую джакобиан в неправильном направлении ...

Actualy, мой код работает, но я думаю, что это было где-то неправильно ... также неявная кривая не движется от начальных значений.

здесь моя функция:

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 

Answer 1

Я немного ржавый на точно, как это реализовать в 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)) 

для данный этап. Теперь идет сложная часть, вам нужно сделать 1-D вектор, построенный путем укладки 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 

Хорошо, теперь мы можем реализовать метод Ньютона. Вы будете делать это для каждого временного шага и до тех пор, пока вы не сходитесь по правильному ответу и не используете его во всех пространственных точках одновременно. Что вы делаете, так это вы догадались K и вычислили якобиан G(K,U) = K-FF(FF(u+dt*a*K). G(K,U) следует оценивать только в нуле, когда K - это правильное решение. Таким образом, другими словами, вы используете метод Ньютона по адресу K, и при поиске конвергенции вы должны видеть, что он сходится во всех точках. Я бы использовал метод Ньютона до max(abs(G(K,U)))< SolverTolerance.

Извините, я не могу больше помочь в реализации matlab, но, надеюсь, мне помогли объяснить, как реализовать метод newton.