solve Nonlinear PDE and compare the analytical and numerical solutions

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')

Answer 1

[Torsten]

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.

[Hana Bachi]

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

[Torsten]

%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]

This is how the solution should look like (as below), but my code shows something different

[Torsten]

Further, since there is an u*du/dx term in your pde, the update from t_j to t_j+1 in the loop can't be just solving a linear equation. The system you have to solve in each time step must be a system of nonlinear equations for which you have to use "fsolve". I don't know how you arrived at the discretization in your code - it's definitely wrong.

[Hana Bachi]

I used the method of burger equation solution, which nonlinear PDE, using crank nicolson method

[Torsten]

There must definitely appear quadratic terms in u in your discretization which result from the term u*du/dx. But you have no quadratic terms in the system you solve - you solve a simple linear equation in each time step.

For which equation did you find the above scheme ? Maybe

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

?

[Hana Bachi]

it was a solution for: du/dt +u* du/dx = d^2u/dx^2

and also from the following Mathwork website:

https://www.mathworks.com/matlabcentral/answers/1589549-solve-1d-burgers-equation-via-semi-implicit-adams-bashforth-crank-nicolson-scheme?s_tid=srchtitle

[Torsten]

What the author has programmed here is an approximation for the solution of the simple diffusion equation

du/dt = d^2u/dx^2

where he also made the mistake that he didn't incorporate the boundary conditions at both sides as you can see from the discontinuities in the curves at x=-6 and x=6.

[Hana Bachi]

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]

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]

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

[Hana Bachi]

The fsolver has an error to be run, I'm getting this error message

'fsolve' requires Optimization Toolbox.

Error in TorstenCp4 (line 26)

unew = fsolve(f,uold);

[Torsten]

You don't have a licence for the optimization toolbox ?

Here is the method on how to check it out:

https://de.mathworks.com/matlabcentral/answers/97495-how-can-i-test-the-license-for-a-toolbox-in-matlab

[Hana Bachi]

I think that I have the license as it shows below

>> license checkout Control_Toolbox

ans =

1

>> S = license('inuse','MATLAB')

S =

struct with fields:

feature: 'matlab'

user: 'oasis'

[Torsten]

Not

license checkout Control_Toolbox

but

license checkout Optimization_Toolbox

[Hana Bachi]

The answer is positive as well

>> license checkout Optimization_Toolbox

ans =

1

[Torsten]

What do you get as output when you enter

ver

?

[Hana Bachi]

I just run it on matlab online and it's working

[Hana Bachi]

form the instaled Matlab, I'm getting this error:

'fsolve' requires Optimization Toolbox.

Error in TorstenCp4 (line 26)

unew = fsolve(f,uold);

>>

[Hana Bachi]

>> ver

-----------------------------------------------------------------------------------------------------

MATLAB Version: 9.11.0.1837725 (R2021b) Update 2

MATLAB License Number: 875352

Operating System: Microsoft Windows 7 Professionnel Version 6.1 (Build 7600)

Java Version: Java 1.8.0_202-b08 with Oracle Corporation Java HotSpot(TM) 64-Bit Server VM mixed mode

-----------------------------------------------------------------------------------------------------

[Torsten]

You should check whether the optimization toolbox is installed because it seems the optimization toolbox is not:

https://de.mathworks.com/matlabcentral/answers/93913-how-can-i-determine-what-add-ons-i-have-installed