solve Nonlinear PDE and compare the analytical and numerical solutions
Show older comments
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')
Accepted Answer
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 Comments
Hana Bachi
on 19 Mar 2022
%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
on 20 Mar 2022
Edited: Hana Bachi
on 20 Mar 2022
This is how the solution should look like (as below), but my code shows something different
Torsten
on 20 Mar 2022
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
on 20 Mar 2022
Torsten
on 20 Mar 2022
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
on 20 Mar 2022
Torsten
on 20 Mar 2022
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
on 20 Mar 2022
Edited: Hana Bachi
on 20 Mar 2022
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.
More Answers (1)
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
10 Comments
Hana Bachi
on 20 Mar 2022
Torsten
on 20 Mar 2022
You don't have a licence for the optimization toolbox ?
Here is the method on how to check it out:
Hana Bachi
on 20 Mar 2022
Torsten
on 20 Mar 2022
Not
license checkout Control_Toolbox
but
license checkout Optimization_Toolbox
Hana Bachi
on 20 Mar 2022
Torsten
on 20 Mar 2022
What do you get as output when you enter
ver
?
Hana Bachi
on 20 Mar 2022
Hana Bachi
on 20 Mar 2022
Hana Bachi
on 20 Mar 2022
Torsten
on 20 Mar 2022
You should check whether the optimization toolbox is installed because it seems the optimization toolbox is not:
Categories
Find more on Eigenvalue Problems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
