solve Nonlinear PDE and compare the analytical and numerical solutions

Question

Hana Bachi 2022-3-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1675644-solve-nonlinear-pde-and-compare-the-analytical-and-numerical-solutions

编辑： Torsten 2022-3-21

I wrote the following code for this problem, but the code stuck at ligne 56. what wrong with is and how can I avoid that problem?

clear all; close all; clc;
% Construct spatial mesh
Nx = 1650;           % Number of spatial steps
xl = 0;            % Left x boundary
xr = 1;             % Right x boundary
dx = (xr-xl)/Nx;    % Spatial step
x = xl:dx:xr;       % x
% Construct temporal mesh
tf = 0.5;             % Final time
dt = 1/3300;        % Timestep
Nt = round(tf/dt);  % Number of timesteps
t = 0:dt:tf;        % t
% Parameters
umax = 15;                      % Found by a perturbation of t=10^-2
C = umax*dt/dx;                 % Convective Stability / Courant Number
disp("Courant Number: "+C)
V = dt/(dx*dx);                 % Viscous Stability / Diffusion Number
disp("Diffusion Number: "+V)
% Dephine phi(x)
phi=zeros(Nx+1,Nt+1);
tt=zeros(1,Nt+1);  
for j=1:Nt
    x(1,j+1)=x(1,j)+dx;
end
for j=1:Nt
    tt(1,j+1)=tt(1,j)+dt;
end
for i=1:Nx+1
    for j=1:Nt+1
        A(j,i)=(50/3)*(x(j)-0.5+4.95*tt(i));
        B(j,i)=(250/3)*(x(j)-0.5+0.75*tt(i));
        C(j,i)=(500/3)*(x(j)-0.375);
        phi(j,i)=(0.1*exp(-A(j,i))+0.5*exp(-B(j,i))+exp(-C(j,i)))/(exp(-A(j,i))+exp(-B(j,i))+exp(-C(j,i)));
    end
end
% Define Boundary & Initial Conditions
% Boundary conditions (Dirichlet)
u(:,1)=phi(:,1);
%left boundary condition
u(1,1)=phi(1,1);
%right boundary condition
u(Nx+1,1)=phi(Nx+1,1);
% Define the right (MMr) and left (MMl) matrices in the Crank-Nicolson method
aal(1:Nx-2) = -V;                   % Below the main diagonal in MMl
bbl(1:Nx-1) = 2+2*V;                % The main diagonal in MMl
ccl(1:Nx-2) = -V;                   % Above the main diagonal in MMl
MMl = diag(bbl,0)+diag(aal,-1)+diag(ccl,1);     % Building the MMl matrix
aar(1:Nx-2) = V;                    % Below the main diagonal in MMr
bbr(1:Nx-1) = 2-2*V;                % The main diagonal in MMr
ccr(1:Nx-2) = V;                    % Above the main diagonal in MMr
MMr = diag(bbr,0)+diag(aar,-1)+diag(ccr,1);     % Building the MMr matrix
% Implementation of the Crank-Nicholson method
for j = 1:Nt                        
    u(2:Nx,j+1) = inv(MMl)*MMr*u(2:Nx,j);   
end
figure()
clf
plot(x,phi(:,ts+1),'b--o')
hold on
%plot(x,phi(:,1001))
xlabel('x')
ylabel('U[x,t]')
title('Analytical Solution')
figure()
clf
plot(x,u(:,:),'g')
hold on
%plot(x,phi(:,1001))
xlabel('x')
ylim([0 1.2])
ylabel('U[x,t]')
title('Numerical Solution')

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2022-3-19

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1675644-solve-nonlinear-pde-and-compare-the-analytical-and-numerical-solutions#answer_921994

编辑：Torsten 2022-3-19

在 MATLAB Online 中打开

In the recursion

% Implementation of the Crank-Nicholson method
for j = 1:Nt                        
    u(2:Nx,j+1) = inv(MMl)*MMr*u(2:Nx,j);   
end

you never address u(1,:) and u(Nx+1,:) which contain the boundary values of your problem.

So this loop cannot be correct to get the solution of your problem.

10 个评论
显示 8更早的评论隐藏 8更早的评论

Torsten 2022-3-19

编辑：Torsten 2022-3-19

在 MATLAB Online 中打开

%clear all; close all; clc;
% Construct spatial mesh
Nx = 1650;           % Number of spatial steps
xl = 0;            % Left x boundary
xr = 1;             % Right x boundary
dx = (xr-xl)/Nx;    % Spatial step
x = xl:dx:xr;       % x
% Construct temporal mesh
tf = 0.5;             % Final time
dt = 1/3300;        % Timestep
Nt = round(tf/dt);  % Number of timesteps
t = 0:dt:tf;        % t
% Parameters
umax = 15;                      % Found by a perturbation of t=10^-2
C = umax*dt/dx;                 % Convective Stability / Courant Number
disp("Courant Number: "+C)
V = dt/(dx*dx);                 % Viscous Stability / Diffusion Number
disp("Diffusion Number: "+V)
% Dephine phi(x)
phi=zeros(Nx+1,Nt+1);
tt=zeros(1,Nt+1);  
for j=1:Nt
    x(1,j+1)=x(1,j)+dx;
end
for j=1:Nt
    tt(1,j+1)=tt(1,j)+dt;
end
for i=1:Nx+1
    for j=1:Nt+1
        A(j,i)=(50/3)*(x(j)-0.5+4.95*tt(i));
        B(j,i)=(250/3)*(x(j)-0.5+0.75*tt(i));
        C(j,i)=(500/3)*(x(j)-0.375);
        phi(j,i)=(0.1*exp(-A(j,i))+0.5*exp(-B(j,i))+exp(-C(j,i)))/(exp(-A(j,i))+exp(-B(j,i))+exp(-C(j,i)));
    end
end
% Define Boundary & Initial Conditions
% Boundary conditions (Dirichlet)
u(:,1)=phi(:,1);
%left boundary condition
u(1,1)=phi(1,1);
%right boundary condition
u(Nx+1,1)=phi(Nx+1,1);
% Define the right (MMr) and left (MMl) matrices in the Crank-Nicolson method
aal(1:Nx-2) = -V;                   % Below the main diagonal in MMl
bbl(1:Nx-1) = 2+2*V;                % The main diagonal in MMl
ccl(1:Nx-2) = -V;                   % Above the main diagonal in MMl
MMl = diag(bbl,0)+diag(aal,-1)+diag(ccl,1);     % Building the MMl matrix
aar(1:Nx-2) = V;                    % Below the main diagonal in MMr
bbr(1:Nx-1) = 2-2*V;                % The main diagonal in MMr
ccr(1:Nx-2) = V;                    % Above the main diagonal in MMr
MMr = diag(bbr,0)+diag(aar,-1)+diag(ccr,1);     % Building the MMr matrix
% Implementation of the Crank-Nicholson method
[L,U] = lu(MMl);
for j = 1:Nt
    u(2:Nx,j+1) = U \ (L \ (MMr*u(2:Nx,j))) ;
end  
%u(2:Nx,2:Nt+1) = MMl\(MMr*u(2:Nx,1:Nt))
%for j = 1:Nt                        
%    u(2:Nx,j+1) = MM1\(MMr*u(2:Nx,j);  %inv(MMl)*MMr*u(2:Nx,j);   
%end
figure()
clf
%plot(x,phi(:,ts+1),'b--o')
plot(x.',phi(:,800),'b--o')
hold on
%plot(x,phi(:,1001))
xlabel('x')
ylabel('U[x,t]')
title('Analytical Solution')
figure()
clf
plot(x,u(800,:),'g')
hold on
%plot(x,phi(:,1001))
xlabel('x')
ylim([0 1.2])
ylabel('U[x,t]')
title('Numerical Solution')

As predicted, this works but gives the wrong solution.

so how can I implement those boundary conditions while I can keep the loop running?

Without diving deep into the CN discretization, I can't tell you exactly how you have to modify your code.

You must write down the discretized system again and compare with your implementation.

My guess is that the line

u(2:Nx,j+1) = U \ (L \ (MMr*u(2:Nx,j))) ;

must read somehow

u(2:Nx,j+1) = U \ (L \ (MMr*u(2:Nx,j) + b)) ;

where b in its first component depends on phi(j,1) and phi(j+1,1), in its last component depends on phi(j,Nx+1) and phi(j+1,Nx+1) and all other components are 0.

Hana Bachi 2022-3-20

编辑：Hana Bachi 2022-3-20

how can that be fixed without rewriting the code from scratch

my eauqtion has a factor of 0.003

du/dt +u* du/dx = 0.003*d^2u/dx^2

Torsten 2022-3-20

编辑：Torsten 2022-3-20

在 MATLAB Online 中打开

It can't be fixed.

At least the part of the code where you implement Crank-Nicholson has to be completely removed, i.e. the part starting with

% Define Boundary & Initial Conditions
% Boundary conditions (Dirichlet)

has to be written anew.

请先登录，再进行评论。

Answer 2

Torsten 2022-3-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1675644-solve-nonlinear-pde-and-compare-the-analytical-and-numerical-solutions#answer_922414

编辑：Torsten 2022-3-21

在 MATLAB Online 中打开

D = 3e-3;
xstart = 0.0;
xend = 1.0;
dx = 1/400;
tstart = 0.0;
tend = 0.5;
dt = 0.01;
  
X = (xstart:dx:xend).';
T = (tstart:dt:tend);
nx = numel(X);
nt = numel(T);
  
U_ana = u_ana(T,X);
  
U = zeros(nx,nt);
U(:,1) = U_ana(:,1);
  
told = T(1);
uold = U(:,1);
  
for j = 2:nt
  tnew = told + dt;
  f = @(u)fun(u,uold,tnew,told,dt,X,nx,dx,D);
  unew = fsolve(f,uold);
  U(:,j) = unew;    
  told = tnew;
  uold = unew;
  j
end
  
plot(X,U(:,10))
hold on
plot(X,U_ana(:,10))  
function res = fun(u,uold,t,told,dt,X,nx,dx,D) 
  res         = zeros(nx,1); 
  res(1)      = u(1) - 0.5*(u_ana(told,X(1)) + u_ana(t,X(1)));
  res(2:nx-1) = (u(2:nx-1)-uold(2:nx-1))/dt - 0.5*( ...
                 (-uold(2:nx-1).*(uold(3:nx)-uold(1:nx-2))/(2*dx) + ...
                 D*(uold(3:nx)-2*uold(2:nx-1)+uold(1:nx-2))/dx^2) + ...
                 (-u(2:nx-1).*(u(3:nx)-u(1:nx-2))/(2*dx) + ...
                 D*(u(3:nx)-2*u(2:nx-1)+u(1:nx-2))/dx^2));
  res(nx)     = u(nx) - 0.5*(u_ana(told,X(nx)) + u_ana(t,X(nx)));
end   
function out = u_ana(t,x)
  A = 50/3*(x-0.5+4.95*t);
  B = 250/3*(x-0.5+0.75*t);
  C = 500/3*(x-0.375);
  out = (0.1*exp(-A)+0.5*exp(-B)+exp(-C))./(exp(-A)+exp(-B)+exp(-C));
end