Прямая и обратная интеграция PDE - PullRequest
Новая
       67

Прямая и обратная интеграция PDE

0 голосов
/ 15 марта 2020

У меня есть следующий PDE с конечным условием и граничным условием, как показано ниже

PDE Equaton

Я пытаюсь интегрировать этот pde назад во времени, однако решение кажется взорванным. Я дискретизировал PDE, используя схему конечных разностей, и попробовал следующий код MATLAB. Однако решение кажется взорванным. Может ли кто-нибудь указать на ошибку в коде, которую я пытаюсь реализовать, или на любую другую причину, по которой решение кажется взорванным. Я попытался интегрировать его во времени, что, кажется, работает нормально. Однако мне необходимо интегрировать его назад во времени. 

%clear all 
clc

%% Defining the variables
v_max = 4;

%%%--- Space Domain----%%%
x_min = 0;     
x_max = 1;      
steps_space=10; %% No of discretization points in space
dx =(x_max-x_min)/steps_space;
x =0:dx:x_max;
x = x(1:end-1);

%%----Time Domain ----%%%
t=0;
t_max=2;
dt = dx/v_max; % To ensure CFL condition
steps_time = int64(t_max/dt); 


%%---- Initial Influx (Input)---%%        
lambda = 2*ones(steps_time,1);

%%%--- Initializing the Density ----%%%

row_final = 1*ones(steps_space,1);%final_con;
row1 = zeros(steps_time,steps_space);
row1(end,:) =[row_final];

%%-- Initial Velocity--%%%
v_final = v_max/(1+trapz(dx,row_final));
v = zeros(steps_time,1);
v(end) = v_final;    

for n= (steps_time-1) :-1:1

    %%% Using the First Order Upwind  

for(i=1:1:steps_space)


    if (i==1)
    %%% Boundary conditoin
    row1(n,i) = row1(n+1,i)+(dt/dx)*((v(n+1)*row1(n+1,i))-lambda(n+1));

    else
    %%% Forward PDE
    row1(n,i) = row1(n+1,i)+(v(n+1)*(dt/dx))*(row1(n+1,i)-row1(n+1,i-1));

    end

end
%%--- Updating the velocity ----%%%
v(n,:) = v_max/(1+trapz(dx,row1(n,:)));

%%--- Updating the outflux ---%%%
outflux1(n) =v(n)*row1(n,end);    


end

1 Ответ

0 голосов
/ 17 марта 2020

Попробуйте эту схему различий

enter image description here

for n = (steps_time-1) :-1:1
    %%% Using the First Order Upwind
    for i = 1:steps_space-1
        row1(n,i) = (row1(n+1,i)*dx+row1(n,i+1)*dt*v(n+1))/(dx+v(n+1)*dt);
    end
    surf(row1)
    pause(0.5)
    %%--- Updating the velocity ----%%%
    v(n,:) = v_max/(1+trapz(x,row1(n,:)));

    %%--- Updating the outflux ---%%%
    outflux1(n) =v(n)*row1(n,end);
end

Результат

enter image description here

...