# Crank Nicolson method to solve a PDE

• I
• hanabachi
In summary, the conversation is about a code written to solve a non-linear PDE using the Crank-Nicolson method. The code has a bug on line 42 and the person is not able to get the correct final results. They are discussing the maximum timestep that can be used and whether they are exceeding it. The code is then posted and it is pointed out that it is not actually solving anything, but just creating a matrix. The person is advised to look at a simple 1D diffusion problem on Wikipedia to understand how the Crank-Nicolson method works. The conversation ends with a request to construct the matrix A and vector b from the given example.
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?

#### Attachments

• CP4.pdf
99.6 KB · Views: 195
There is a bug on line 42

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

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

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

DrClaude said:
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
Matlab:
%% 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>

Last edited by a moderator:
pasmith said:
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

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?

actually the equation is

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

hanabachi said:
actually the equation is

dU/dt=0U*du/dx+0.003*dU2/dx2
This is linear. I assume 0U is a typo though?

hanabachi said:
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.

joshmccraney said:
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

bigfooted said:
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.

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

#### Attachments

• DocScanner Feb 18, 2022 7-56 PM.pdf
224.2 KB · Views: 160
Can you please type it here in latex? My eyes are still hurting from looking at this document :-)

hunt_mat, fresh_42 and (deleted member)
hanabachi said:
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/

jim mcnamara and member 428835
bigfooted said:
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}

jim mcnamara and fresh_42
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$$

bigfooted said:
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.

bigfooted said:
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

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

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

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

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

The point is that for interior points the Crank-Nicholson method reduces to $$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 $$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.

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.

## What is the Crank Nicolson method?

The Crank Nicolson method is a numerical method used to solve partial differential equations (PDEs). It is a combination of the explicit Euler method and the implicit Euler method, and is known for its stability and accuracy.

## How does the Crank Nicolson method work?

The Crank Nicolson method uses a finite difference approximation to discretize the PDE and then solves the resulting system of equations using a tridiagonal matrix algorithm. This method takes into account both the current and future time steps, resulting in a more accurate solution.

## What types of PDEs can the Crank Nicolson method solve?

The Crank Nicolson method can be used to solve a wide range of PDEs, including parabolic, hyperbolic, and elliptic equations. It is particularly useful for solving diffusion and heat transfer problems.

## Is the Crank Nicolson method always more accurate than other numerical methods?

No, the accuracy of the Crank Nicolson method depends on the specific PDE being solved and the chosen time and space discretization. In some cases, other methods may be more accurate. It is important to carefully consider the problem at hand and choose the most appropriate method.

## What are the advantages of using the Crank Nicolson method?

The Crank Nicolson method has several advantages, including its stability, accuracy, and ability to handle a wide range of PDEs. It is also relatively easy to implement and can handle non-uniform grids. Additionally, it is a second-order method, meaning that it converges to the exact solution at a faster rate than first-order methods.

• Differential Equations
Replies
4
Views
1K
• Differential Equations
Replies
6
Views
2K
• Differential Equations
Replies
2
Views
2K
• Differential Equations
Replies
5
Views
2K
• Differential Equations
Replies
2
Views
1K
• Differential Equations
Replies
8
Views
4K
• Differential Equations
Replies
12
Views
8K
• Differential Equations
Replies
1
Views
2K
• Differential Equations
Replies
1
Views
1K
• Differential Equations
Replies
1
Views
3K