- #1

- 12

- 2

You are using an out of date browser. It may not display this or other websites correctly.

You should upgrade or use an alternative browser.

You should upgrade or use an alternative browser.

- #1

- 12

- 2

- #2

Mentor

- 8,241

- 5,123

There is a bug on line 42

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

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

- #3

- 12

- 2

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

- #4

Homework Helper

2022 Award

- 2,673

- 1,280

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

- #5

- 12

- 2

Here is the codeThere is a bug on line 42

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

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:

- #6

- 12

- 2

I used this dt = 1/3300What's the maximum timestep you can use and still have the method be stable? Are you exceeding it?

- #7

Gold Member

- 677

- 198

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?

- #8

- 12

- 2

actually the equation is

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

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

- #9

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

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

- #10

Gold Member

- 677

- 198

Actually the partial differential equation is:actually the equation is

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

$$ \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.

- #11

- 12

- 2

0 should be negative signs (-), so the equation is the multiplication of U by its first derivative dU/dx. So the equation is not linearThis is linear. I assume 0U is a typo though?

- #12

- 12

- 2

$$ \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.

- #13

- 12

- 2

- #14

Gold Member

- 677

- 198

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

- #15

- 18,267

- 20,231

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.Yes, this is what I'm using to solve the pde. please check the attached file

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

- #16

- 12

- 2

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

\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}

- #17

Gold Member

- 677

- 198

row 1: $$[A_{11} A_{21} ... A_{N1}] \cdot U_1 = b_1$$

- #18

- 12

- 2

$$ \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, that what it should looks like

row 1: $$[A_{11} A_{21} ... A_{N1}] \cdot U_1 = b_1$$

- #19

Gold Member

- 677

- 198

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

- #20

- 12

- 2

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

- #21

Gold Member

- 677

- 198

But the big problem is the definition of the matrix A. What are the coefficients that belong to [itex]U_i[/itex], [itex]U_{i+1}[/itex] and [itex]U_{i-1}[/itex]?

Specifically, what do you do with the nonlinear term

[itex]-r*\Delta x*U_{i}^{n+1}\left( U_{i}^{n+1} -U_{i-1}^{n+1}\right)[/itex]

- #22

- 12

- 2

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

- #23

Homework Helper

2022 Award

- 2,673

- 1,280

U_i^{n+1} - \tfrac12\Delta t f_i(U^{n+1}) = U_i^{n} + \tfrac12\Delta t f_i(U^{n})[/tex] where [itex]f_i = \frac{\partial U_i}{\partial t}[/itex]. If [itex]f[/itex] is non-linear, then you don't get a matrix equation; you get a non-linear algebraic system [tex]

F_i(U^{n+1}) = U_i^{n} + \tfrac12\Delta t f_i(U^{n})[/tex] which has to be solved by, for example, Newton's method.

- #24

Gold Member

- 677

- 198

Share:

- Replies
- 2

- Views
- 964

- Replies
- 7

- Views
- 335

- Replies
- 6

- Views
- 2K

- Replies
- 4

- Views
- 370

- Replies
- 5

- Views
- 1K

- Replies
- 3

- Views
- 408

- Replies
- 13

- Views
- 2K

- Replies
- 8

- Views
- 3K

- Replies
- 2

- Views
- 321

- Replies
- 7

- Views
- 437