Return to Answer

Added more information also to address a suggestion in the comment

Source Link

edited May 31, 2016 at 14:17

1.3k
6
8

I was curious about your problem and it was not difficult for me to test CN method for your two equations. Added - parts beginning with "added" were written afterwards to address the comment.

The point given by David is important. Here are three plots for initial condition (added - that do not fulfill the suggested condition $\alpha u_{xx} + \beta v_{xx}=0$), the first time step and the 5th one, computed by implicit Euler method, the blue curve is u, the other one is v:

The diffusion coefficients $\alpha=\beta=1$ are constant and equal 1. Added - The solution changes rapidly in the first time step for arbitrary small time step that is unusual. Afterwards standard numerical diffusion can be observed (the 5th time step).

The same computations with Crank-Nicolson are strange as you observed (the 1st and 5th time step):

I suppose someone can give you later comment on your PDEs directly, but from this simplified numerical analysis I would guess that you have no chance to obtain reasonable solution with Crank-Nicolson method (added - if your initial functions do not fulfill $(\alpha u_x)_x + (\beta v_x)_x=0$). I doubt if the numerical solution with implicit Euler method is reasonable, but if yes, it suggests that for initial condition that do not fulfill the conditionthis PDE described by David, you obtain for arbitrary small time a jump in the solution, so you have some "infinite" speed of evolution. Added - The situation is even worse for explicit Euler methods as it is noted in the comment.

This is only heuristic answer, but at least it confirms your doubts about Crank-Nicolson method for your problem. HereBelow is the code in Mathematica that produce the pictures:.

Added - if your initial functions fulfill the required stationary PDE (but in the approximative form using e.g. finite difference method as in the provided code) then you may use implicit ($\theta=1$ in the code) or explicit ($\theta=0$) Euler method or Crank-Nicolson method ($\theta=0.5$) for time discretization, they will work because these methods preserve this property also for evolved numerical functions if it is fulfilled for initial functions.

I was curious about your problem and it was not difficult for me to test CN method for your two equations.

The point given by David is important. Here are three plots for initial condition, the first time step and the 5th one, computed by implicit Euler method, the blue curve is u, the other one is v:

The diffusion coefficients are constant and equal 1.

The computations with Crank-Nicolson are strange as you observed (the 1st and 5th time step):

I suppose someone can give you later comment on your PDEs directly, but from this simplified numerical analysis I would guess that you have no chance to obtain reasonable solution with Crank-Nicolson method. I doubt if the numerical solution with implicit Euler method is reasonable, but if yes, it suggests that for initial condition that do not fulfill the condition described by David you obtain for arbitrary small time a jump in the solution, so you have some "infinite" speed of evolution.

This is only heuristic answer, but at least it confirms your doubts about Crank-Nicolson method for your problem. Here is the code in Mathematica that produce the pictures:

I was curious about your problem and it was not difficult for me to test CN method for your two equations. Added - parts beginning with "added" were written afterwards to address the comment.

The diffusion coefficients $\alpha=\beta=1$ are constant. Added - The solution changes rapidly in the first time step for arbitrary small time step that is unusual. Afterwards standard numerical diffusion can be observed (the 5th time step).

The same computations with Crank-Nicolson are strange as you observed (the 1st and 5th time step):

I suppose someone can give you later comment on your PDEs directly, but from this simplified numerical analysis I would guess that you have no chance to obtain reasonable solution with Crank-Nicolson method (added - if your initial functions do not fulfill $(\alpha u_x)_x + (\beta v_x)_x=0$). I doubt if the numerical solution with implicit Euler method is reasonable, but if yes, it suggests that for initial condition that do not fulfill this PDE described by David, you obtain for arbitrary small time a jump in the solution, so you have some "infinite" speed of evolution. Added - The situation is even worse for explicit Euler methods as it is noted in the comment.

This is only heuristic answer, but at least it confirms your doubts about Crank-Nicolson method for your problem. Below is the code in Mathematica that produce the pictures.

Source Link

created May 30, 2016 at 21:45

Peter Frolkovič

1.3k
6
8

I was curious about your problem and it was not difficult for me to test CN method for your two equations.

The point given by David is important. Here are three plots for initial condition, the first time step and the 5th one, computed by implicit Euler method, the blue curve is u, the other one is v:

The diffusion coefficients are constant and equal 1.

The computations with Crank-Nicolson are strange as you observed (the 1st and 5th time step):

I suppose someone can give you later comment on your PDEs directly, but from this simplified numerical analysis I would guess that you have no chance to obtain reasonable solution with Crank-Nicolson method. I doubt if the numerical solution with implicit Euler method is reasonable, but if yes, it suggests that for initial condition that do not fulfill the condition described by David you obtain for arbitrary small time a jump in the solution, so you have some "infinite" speed of evolution.

This is only heuristic answer, but at least it confirms your doubts about Crank-Nicolson method for your problem. Here is the code in Mathematica that produce the pictures:

Model data
t0 = 0.; T = 0.1; L = 10.; alpha = 1.; beta = 1.; 
u0[x_] := If[Abs[x - 3] <= 1., 1., 0.]; 
v0[x_] := If[Abs[x - 7] <= 1., -1., 0.]; 

Method's input
Nd = 5; M = 20; theta = 0.5; \[Tau] = T/Nd; h = L/(M + 1); 
Mat = SparseArray[{}, {2*M, 2*M}]; 
ps = Array[{}, {2*M}]; 

Method
For[i = 1, i <= M, i++, If[i != 1, Mat[[i,i - 1]] = (-alpha)*\[Tau]*theta]; Mat[[i,i]] = h^2 + 2*alpha*\[Tau]*theta; 
   If[i != M, Mat[[i,i + 1]] = (-alpha)*\[Tau]*theta]; Mat[[i,i + M]] = -h^2; ]
For[i = M + 1, i <= 2*M, i++, If[i != M + 1, Mat[[i,i - 1]] = (-beta)*\[Tau]*theta]; Mat[[i,i]] = h^2 + 2*beta*\[Tau]*theta; 
   If[i != 2*M, Mat[[i,i + 1]] = (-beta)*\[Tau]*theta]; Mat[[i,i - M]] = -h^2; ]

solver = LinearSolve[Mat]; 

For[i = 0, i <= M + 1, i++, u[i, 0] = u0[i*h]; v[i, 0] = v0[i*h]]

For[n = 1, n <= Nd, n++, 
  For[i = 1, i <= M, i++, ps[[i]] = h^2*(u[i, n - 1] - v[i, n - 1]) + (1 - theta)*\[Tau]*alpha*(u[i - 1, n - 1] - 2*u[i, n - 1] + 
         u[i + 1, n - 1]); ps[[M + i]] = h^2*(v[i, n - 1] - u[i, n - 1]) + (1 - theta)*\[Tau]*beta*(v[i - 1, n - 1] - 2*v[i, n - 1] + 
         v[i + 1, n - 1])]; r = solver[ps]; u[0, n] = u[M + 1, n] = v[0, n] = v[M + 1, n] = 0.; 
   For[i = 1, i <= M, i++, u[i, n] = r[[i]]; v[i, n] = r[[M + i]]]; ]

n = 0; ListLinePlot[{Table[{i*h, u[i, n]}, {i, 0, M + 1}], Table[{i*h, v[i, n]}, {i, 0, M + 1}]}, PlotRange -> All]
n = 1; ListLinePlot[{Table[{i*h, u[i, n]}, {i, 0, M + 1}], Table[{i*h, v[i, n]}, {i, 0, M + 1}]}, PlotRange -> All]
n = 5; ListLinePlot[{Table[{i*h, u[i, n]}, {i, 0, M + 1}], Table[{i*h, v[i, n]}, {i, 0, M + 1}]}, PlotRange -> All]