I Crank Nicolson method to solve a PDE

[hanabachi]

Hello,

I wrote a code to solve a non-linear PDE using Canrk nicolson method, but I'm still not able to get a correct final results. can anyone tell me what wrong with it?

 

[DrClaude]

There is a bug on line 42

If you don't post the code, there is not much anyone can do.

 

[hanabachi]

I though the code is visible in the attached files. Thank you for letting me know

 

[pasmith]

What's the maximum timestep you can use and still have the method be stable? Are you exceeding it?

 

[hanabachi]

There is a bug on line 42

If you don't post the code, there is not much anyone can do.

Here is the code

%% Crank-Nicolson Method
clear variables
close all
% 1. Space steps
xa = 0;
xb = 1;
dx = 1/40;
N = (xb-xa)/dx ;
x = xa:dx:xb;
%2.Time steps
ta = 0;
tb = 0.5;
dt = 1/3300;
M = (tb-ta)/dt ;
t = ta:dt:tb;
%3. Controling Parameters

%4. Define equations A, B , C and phi(x,t)
A = @(x,t) (50/3)*(x-0.5+4.95*t);
B = @(x,t) (250/3)*(x-0.5+0.75*t);
C = @(x,t) (500/3)*(x-0.375);
phi = @(x,t)(0.1*exp(-A(x,t)) +0.5*exp(-B(x,t))+exp(-C(x,t)))./...
(exp(-A(x,t)) +exp(-B(x,t))+exp(-C(x,t)));
% 5 Initial and boundary conditions
f = @(x) phi(x,t(1)); % Initial condition
g1 = @(t)phi(x(1),t); % Left boundary condition
g2 = @(t)phi(x(N+1),t); % Right boundary condition

r = dt / (2*dx^2);
a = -0.003*r;
b = 1 + 2*0.003*r;
c = -r*0.003;
d =(1+0.003*2*r);
e =1+0.003*2*r;
h =0.003*r;
% 6 Implementation of the explicit method
U = zeros(N+1,M+1);
U(2:N,1) = f(x(2:N)); % Put in the initial condition
U(1,:) = g1(t); % The boundary conditions, g1 and g2 at x = 0 and x = 1
U(N+1,:) = g2(t);

for i = 2:N   
     RHS= (1+2*0.003*r)*U(i,N)+0.003*r*U(i+1,N)+0.003*r*U(i-1,N)-r*dx*U(i,N)*(U(i,N)-U(i-1,N));
end
for j=2:M % Time Loop
    for i= 2:N % Space Loop
       RHS= (1+2*0.003*r)*U(i,N+1)+0.003*r*U(i+1,N+1)+0.003*r*U(i-1,N+1)-r*dx*U(i,N+1)*(U(i,N+1)-U(i-1,N+1));
    end
end
% Make some plots
T= 0:.1:M;
V = [];
for i= 1: length(T)
    P = find(t==T(i));
    V = [V P];
end
figure
subplot(131)
for j = 1:length(V)
hold on
plot(x,U(:,V(j)),'*-','linewidth',2.5,'DisplayName',sprintf('t = %1.4f',t(V(j))))
end
legend('-DynamicLegend','location','bestoutside');
a = ylabel('U-Values');
set(a,'Fontsize',14);
a = xlabel('X-Values');
set(a,'Fontsize',14);
a=title('Crank-Nicolson Solution');
set(a,'Fontsize',16);
grid;
% disp(u(:,V)'); % Each row corresponds to a particular value of t and
% Each column corresponds to a particular value of x
% Implement the exact solution and compare it to the exact solution
Exact =@(x,t) phi(x,t);
subplot(132)
for j = 1:length(V)
hold on
plot(x,Exact(x,t(V(j))),'*-','linewidth',2.5,'DisplayName',sprintf('t = %1.4f',t(V(j))))
end
legend('-DynamicLegend','location','bestoutside');
a = ylabel('U-Values');
set(a,'Fontsize',14);
a = xlabel('X-Values');
set(a,'Fontsize',14);
a=title(' Analytical Solution');
set(a,'Fontsize',16);
grid;
[X,T] = meshgrid(x,t);
% Exact2 = (0.1*exp(-A(X,T)) +0.5*exp(-B(X,T))+exp(-C(X,T)))./...
% (exp(-A(X,T)) +exp(-B(X,T))+exp(-C(X,T)));
Error =abs(Exact(X,T)'-U);
subplot(133)
for j = 1:length(V)
hold on
plot(x,Error(:,(V(j))),'*-','linewidth',2.5,'DisplayName',sprintf('t = %1.4f',t(V(j))))
end
legend('-DynamicLegend','location','bestoutside');
a = ylabel('Error');
set(a,'Fontsize',14);
a = xlabel('X-Values');
set(a,'Fontsize',14);
a=title(' Absolute Error ');
set(a,'Fontsize',16);
grid;

<mentor edit code tags>

 

[hanabachi]

What's the maximum timestep you can use and still have the method be stable? Are you exceeding it?

I used this dt = 1/3300

 

[bigfooted]

Well, to me it looks like you're not actually solving anything, you've just created a 2D matrix U and filled it.
You need to solve a matrix-vector system $$Ax=b$$ for your entire spatial domain at every time step. What is your matrix A and your vector b? Do you have lecture notes explaining how the Crank-Nicolson method works?
Have a look at the very simple 1D diffusion problem at wikipedia:
https://en.wikipedia.org/wiki/Crank–Nicolson_method
In this example you try to solve the diffusion problem
$$\frac{\partial u}{\partial t} = a \frac{\partial^2 u}{\partial x^2} $$, and the discretization scheme leads to:
$$-r u_{i+1}^{n+1} + (1+2r)u_i^{n+1}-ru_{i-1}^{n+1}=ru_{i+1}^{n}+(1-2r)u_i^{n}+ru_{i-1}^n, $$
with
$$r=\frac{a\Delta t}{2(\Delta x) ^2}$$.
The temporal and diffusion terms adds to the matrix A as well as the vector b as you move everything known at timestep t to the right and you keep the unknown terms at timestep t+1 in the matrix at the left. Can you construct the matrix A and vector b from this example?

 

[hanabachi]

actually the equation is

dU/dt=0U*du/dx+0.003*dU2/dx2

 

[member 428835]

actually the equation is

dU/dt=0U*du/dx+0.003*dU2/dx2

This is linear. I assume 0U is a typo though?

 

[bigfooted]

actually the equation is

dU/dt=0U*du/dx+0.003*dU2/dx2

Actually the partial differential equation is:
$$ \frac{\partial u}{\partial t} = -u\frac{\partial u}{\partial x} + 0.003 \frac{\partial^2 u }{\partial x^2}, x \in (0,1), t>0$$

But the question is to solve it using the Crank-Nicolson method. Do you have the difference scheme for this PDE and can you reproduce it here? We can continue from there.

 

[hanabachi]

This is linear. I assume 0U is a typo though?

0 should be negative signs (-), so the equation is the multiplication of U by its first derivative dU/dx. So the equation is not linear

 

[hanabachi]

Actually the partial differential equation is:
$$ \frac{\partial u}{\partial t} = -u\frac{\partial u}{\partial x} + 0.003 \frac{\partial^2 u }{\partial x^2}, x \in (0,1), t>0$$

But the question is to solve it using the Crank-Nicolson method. Do you have the difference scheme for this PDE and can you reproduce it here? We can continue from there.

 

[hanabachi]

Yes, this is what I'm using to solve the pde. please check the attached file

 

[bigfooted]

Can you please type it here in latex? My eyes are still hurting from looking at this document :-)

 

[fresh_42]

Yes, this is what I'm using to solve the pde. please check the attached file

Please be aware that a lot of people simply refuse to download data trash, i.e. documents that have to be deleted after use. I, for example. If it is too troublesome for you to type it out here, then it is too troublesome for me to download whatever you want to force me to.

Here is explained how you can type formulas on PF: https://www.physicsforums.com/help/latexhelp/

 

[hanabachi]

Can you please type it here in latex? My eyes are still hurting from looking at this document :-)

yes, here it is

\begin{array}{l}
\frac{U_{i}^{n+1/2} -U_{i}^{n}}{\Delta t/2} =-U_{i}^{n}\frac{U_{i}^{n} -U_{i-1}^{n}}{\Delta x} \ +0.003\ \frac{U_{i+1}^{n} -2*U_{i}^{n} +U_{i-1}^{n}}{\Delta x^{2}}\\
\\
\frac{U_{i}^{n+1} -U_{i}^{n+1/2}}{\Delta t/2} =-U_{i}^{n+1}\frac{U_{i}^{n+1} -U_{i-1}^{n+1}}{\Delta x} \ +0.003\ \frac{U_{i+1}^{n+1} -2*U_{i}^{n+1} +U_{i-1}^{n+1}}{\Delta x^{2}}\\
\\
adding\ ( 1) \ and\ ( 2)\\
r=\frac{\Delta t}{2\Delta x^{2}}\\
U_{i}^{n+1} +r*\Delta x*U_{i}^{n+1}\left( U_{i}^{n+1} -U_{i-1}^{n+1}\right) -0.003*r\left( U_{i+1}^{n+1} -2U_{i}^{n+1} +U_{i-1}^{n+1}\right) =\\
\ U_{i}^{n} -r*\Delta x*U_{i}^{n}\left( U_{i}^{n} -U_{i-1}^{n}\right) +0.003*r\left( U_{i+1}^{n} -2U_{i}^{n} +U_{i-1}^{n}\right)
\end{array}

 

[bigfooted]

OK, nice. So the left-hand side creates the matrix A and the right hand side created the vector b. Can you show now how to construct a row in matrix A and vector b? So what does the first and second row look like for instance?
row 1: $$[A_{11} A_{21} ... A_{N1}] \cdot U_1 = b_1$$

 

[hanabachi]

Actually the partial differential equation is:
$$ \frac{\partial u}{\partial t} = -u\frac{\partial u}{\partial x} + 0.003 \frac{\partial^2 u }{\partial x^2}, x \in (0,1), t>0$$

But the question is to solve it using the Crank-Nicolson method. Do you have the difference scheme for this PDE and can you reproduce it here? We can continue from there.

OK, nice. So the left-hand side creates the matrix A and the right hand side created the vector b. Can you show now how to construct a row in matrix A and vector b? So what does the first and second row look like for instance?
row 1: $$[A_{11} A_{21} ... A_{N1}] \cdot U_1 = b_1$$

Yes, that what it should looks like

 

[bigfooted]

But what are the exact values in your case for the A's and the b's?

 

[hanabachi]

b is RHS and A should be the factors a,b and c that I forgot to replace them in this experssion
RHS= (1+2*0.003*r)*U(i,N)+0.003*r*U(i+1,N)+0.003*r*U(i-1,N)-r*dx*U(i,N)*(U(i,N)-U(i-1,N));

 

[bigfooted]

There is a minus sign wrong in your RHS.
But the big problem is the definition of the matrix A. What are the coefficients that belong to U_i, U_{i+1} and U_{i-1}?

Specifically, what do you do with the nonlinear term
-r*\Delta x*U_{i}^{n+1}\left( U_{i}^{n+1} -U_{i-1}^{n+1}\right)

 

[hanabachi]

the nonlinear terms comes from the nonlinear PDE, and that what I got from Crank-Nicolson method

 

[pasmith]

The point is that for interior points the Crank-Nicholson method reduces to <br /> U_i^{n+1} - \tfrac12\Delta t f_i(U^{n+1}) = U_i^{n} + \tfrac12\Delta t f_i(U^{n}) where f_i = \frac{\partial U_i}{\partial t}. If f is non-linear, then you don't get a matrix equation; you get a non-linear algebraic system <br /> F_i(U^{n+1}) = U_i^{n} + \tfrac12\Delta t f_i(U^{n}) which has to be solved by, for example, Newton's method.

 

[bigfooted]

Exactly. So at this point you have to decide what to do with the nonlinear terms. Most common is to linearize it and use a fixed point iterator (simplest) or Newton's method (better). Then, instead of solving Au=b for u, you start with an estimate for u, and substitute it in the function F(u)=A(u)\cdot u -b(u) (A(u) is now depending on u) and you try to find u that minimizes F.