Crank–Nicolson Stability for Forward and Backward Time Problems
Article Quant Q&A · Author: Si Mo
Summary
The document tests a Crank–Nicolson finite-difference scheme on a linear parabolic equation with a known exponential solution. With exact initial and spatial boundary values, reported mean absolute error falls as the spatial and time grids are refined while keeping their squared-grid ratio fixed. The author then explores two modifications: replacing exact spatial boundary values with third-order smoothing conditions, and solving from a terminal condition by stepping backward in time or transforming the time coordinate.
The reported smoothing results show errors decreasing only slowly, while both backward approaches produce very large errors or numerical overflow. The examples suggest that correct boundary equations, coefficient signs, and time direction matter, and that a scheme working in the forward initial-value setting does not guarantee a stable backward solve. The document does not establish a general stability result or isolate the precise coding errors; the script is partial, and the displayed results alone cannot distinguish implementation mistakes from ill-conditioning intrinsic to the problem.
Key ideas
- A known exact solution provides a way to assess a finite-difference implementation.
- The reported forward results improve as the grids are refined at a fixed time-to-space scaling ratio.
- Third-order smoothing boundary conditions produce much larger errors than exact boundary values in this example.
- Backward stepping and a time-reversed PDE both show explosive errors, so their discretization and boundary treatment need careful review.
- The numerical experiments do not prove that the backward problems are intrinsically unstable.
Tags
Full text
# crank nicolson - forward vs backward equations
# crank nicolson - forward vs backward equations
I have implemented a script in python to solve a differential equation in $(x,t)$ using the crank-nicolson method. To start with, I was testing a case in which I know a solution: $$u_{xx} + u_{x} - u - u_{t} = 0, \ \mbox{ with } \ X_{0} \leq x \leq X_{1}, \ T_{0} \leq t \leq T_{1}.$$ This is solved by $$u(x, t) = exp(x+t),$$ if we impose appropriate boundary conditions.
I was relying on this paper. Though, looking at equations (16/17), the role of $A_i,$ $C_i$ seem to be inverted, so I have inverted them. The article refers to the case in which you have a condition at $T_0$, i.e. $u(x,T_0) = f(x),$ so I used the boundary condition: $$u(x,T_0) = exp(x+T_0).$$ I also used the boundaries at the upper/lower $x-$ axis: $$u(X_0,t) = exp(X_0+t), \ \ u(X_1,t) = exp(X_1+t),$$ using first and last equations in (18).
Until here, everything works fine. Keeping fixed $k/h^2$ and increasing $k$, the errors between the solution and the true solution $(u = exp(x+t))$ converge to zero. Here is the output:
```
N.B. N_x = 10 * multiplier, N_t = 4 * multiplier
multiplier: 1 , N_x: 10 , N_t: 4 , k: 0.500 , h^2: 0.6400 , k/h^2: 0.781 --> mean abs err is: 0.9999
multiplier: 3 , N_x: 30 , N_t: 36 , k: 0.056 , h^2: 0.0711 , k/h^2: 0.781 --> mean abs err is: 0.1109
multiplier: 5 , N_x: 50 , N_t: 100 , k: 0.020 , h^2: 0.0256 , k/h^2: 0.781 --> mean abs err is: 0.0403
multiplier: 10 , N_x: 100 , N_t: 400 , k: 0.005 , h^2: 0.0064 , k/h^2: 0.781 --> mean abs err is: 0.0101
multiplier: 30 , N_x: 300 , N_t: 3600 , k: 0.001 , h^2: 0.0007 , k/h^2: 0.781 --> mean abs err is: 0.0011
multiplier: 100 , N_x: 1000 , N_t: 40000 , k: 0.000 , h^2: 0.0001 , k/h^2: 0.781 --> mean abs err is: 0.0001
```
I also wanted to try a couple of variations, for which I am getting problems:
- The first variation is to use the smoothing conditions at the upper/lower $x-$boundaries, as suggested by the paper. Basically, first and last equations become: $u_{0,n+1} - 3\cdot u_{1,n+1} + 3\cdot u_{2,n+1} -u_{3,n+1} = 0;$ $u_{I-3,n+1} - 3\cdot u_{I-2,n+1} + 3\cdot u_{I-1,n+1} -u_{I,n+1} = 0;$
- The second variation is to impose the boundary condition at $T_1$ rather than $T_0$ and solve backward. For this, I have tried two methods: The first method is to run the same PDEs, but backwards. Basically, in eq. (16), the recursively known vector $D_i$ becomes the LHS rather than the RHS and we solve for $u_{.,n}$ with $u_{.,n+1}$ known from previous step. I also changed the boundary conditions replacing the conditions at $t$ with $T_1-t$, i.e.: $b.1) \ u(x, T_1) = exp(x+T_1);$ Then, recursively at time-step $j$: $b.2) \ u(X_0, (T_1-j)*d\_t) = exp(X_0+(T_1-j)*d\_t);$ $b.3) \ u(X_1, (T_1-j)*d\_t) = exp(X_1+(T_1-j)*d\_t).$ The second method is to define $\tilde{u}(x, t) = u(x, T-t),$ which solves a slight variation of the original PDE. Basically, the sign of $\tilde{u}_t$ changes: $$\tilde{u}_{xx} + \tilde{u}_{x} - \tilde{u} + \tilde{u}_{t} = 0.$$ Taking this into account and going through the calculations eq. 15) and 16) - unless I have made wrong calculations - what changes is that the $4\cdot h^2$ term within $B_i$ in eq. 16) - both RHS and LHS - change sign. As it comes to the boundary conditions, they are the same as $b.1)-b.3)$ above except that we need to make the change $t \rightarrow T_1-t$ only in the RHS of the equations and not in the $\tilde{u}(x, .)$, i.e: $b.1) \ \tilde{u}(x, T_0) = exp(x+T_1);$ Then, recursively at time-step $j$: $b.2) \ \tilde{u}(X_0, j*d_t) = exp(X_0+(T_1-j)*d\_t);$ $b.3) \ \tilde{u}(X_1, j*d_t) = exp(X_1+(T_1-j)*d\_t).$ Then our solution is simply $u(x,t) = \tilde{u}(x, T_1-t)$
When it comes to 1. (smoothing conditions), the convergence is very slow, and I am not even sure it is converging in the first place:
```
multiplier: 1 , N_x: 10 , N_t: 4 , k: 0.500 , h^2: 0.6400 , k/h^2: 0.781 --> mean abs err is: 21.4295
multiplier: 3 , N_x: 30 , N_t: 36 , k: 0.056 , h^2: 0.0711 , k/h^2: 0.781 --> mean abs err is: 14.7544
multiplier: 5 , N_x: 50 , N_t: 100 , k: 0.020 , h^2: 0.0256 , k/h^2: 0.781 --> mean abs err is: 13.5737
multiplier: 10 , N_x: 100 , N_t: 400 , k: 0.005 , h^2: 0.0064 , k/h^2: 0.781 --> mean abs err is: 12.7183
multiplier: 30 , N_x: 300 , N_t: 3600 , k: 0.001 , h^2: 0.0007 , k/h^2: 0.781 --> mean abs err is: 12.1621
multiplier: 100 , N_x: 1000 , N_t: 40000 , k: 0.000 , h^2: 0.0001 , k/h^2: 0.781 --> mean abs err is: 11.9701
```
When it comes to 2. (boundary condition at $T_1$), the solution seems to explode:
```
forward or backward: backward , x_bound_cond: exact_solution , method_if_backward: backward_equations
multiplier: 1 , N_x: 10 , N_t: 4 , k: 0.500 , h^2: 0.6400 , k/h^2: 0.781 --> mean abs err is: 14957.2237
multiplier: 3 , N_x: 30 , N_t: 36 , k: 0.056 , h^2: 0.0711 , k/h^2: 0.781 --> mean abs err is: 7681.0148
multiplier: 5 , N_x: 50 , N_t: 100 , k: 0.020 , h^2: 0.0256 , k/h^2: 0.781 --> mean abs err is: 11680.6986
multiplier: 10 , N_x: 100 , N_t: 400 , k: 0.005 , h^2: 0.0064 , k/h^2: 0.781 --> mean abs err is: 22387.3289
multiplier: 30 , N_x: 300 , N_t: 3600 , k: 0.001 , h^2: 0.0007 , k/h^2: 0.781 --> mean abs err is: 66052.2187
-----------------------------------
forward or backward: backward , x_bound_cond: exact_solution , method_if_backward: y_tilde
multiplier: 1 , N_x: 10 , N_t: 4 , k: 0.500 , h^2: 0.6400 , k/h^2: 0.781 --> mean abs err is: 20365556.8251
multiplier: 3 , N_x: 30 , N_t: 36 , k: 0.056 , h^2: 0.0711 , k/h^2: 0.781 --> mean abs err is: 314523381202040153368016313980598438510416920683660041464578048.0000
multiplier: 5 , N_x: 50 , N_t: 100 , k: 0.020 , h^2: 0.0256 , k/h^2: 0.781 --> mean abs err is: 2775367546745822588503047993652382588512277657755480108699830290855659378104663930014105617681772368669677168984188472128162486164535007131103223952412569027060519374860855562370817824233581897145388919422976.0000
<ipython-input-6-e291cb7cd65f>:58: RuntimeWarning: overflow encountered in double_scalars
d[i] = d[i] - w * d[i-1]
<ipython-input-6-e291cb7cd65f>:58: RuntimeWarning: invalid value encountered in double_scalars
d[i] = d[i] - w * d[i-1]
multiplier: 10 , N_x: 100 , N_t: 400 , k: 0.005 , h^2: 0.0064 , k/h^2: 0.781 --> mean abs err is: nan
<ipython-input-6-e291cb7cd65f>:68: RuntimeWarning: overflow encountered in double_scalars
x[i] = (d[i] - c[i]*x[i+1])/b[i]
<ipython-input-6-e291cb7cd65f>:68: RuntimeWarning: invalid value encountered in double_scalars
x[i] = (d[i] - c[i]*x[i+1])/b[i]
<ipython-input-6-e291cb7cd65f>:324: RuntimeWarning: invalid value encountered in double_scalars
D[1: N_x] = np.array([(A_prev[i-1] * y_tilde_prev[i - 1] + B_prev[i-1] * y_tilde_prev[i] + C_prev[i-1] * y_tilde_prev[i + 1]) for i in range(1, N_x)],
multiplier: 30 , N_x: 300 , N_t: 3600 , k: 0.001 , h^2: 0.0007 , k/h^2: 0.781 --> mean abs err is: nan
-----------------------------------
```
I cannot understand if I am getting anything wrong, or if these problems are intrinsically numerically unstable...
P.S. The script is not very short, but for reference I report it below.
Thanks a lot for help!
```
from math import exp
import numpy as np
# it returns a tridiagonal matrix of the form:
#|b[0] c[0] |
#|a[0] b[1] c[1] |
#| . . . |
#| . . . |
#| . . . |
#| . . . |
#| a[n-1] b[n] c[n] |
#| a[n] b[n+1]|
def tridiag_matrix(a, b, c, k_1 = -1, k_2 = 0, k_3 = 1):
return np.diag(a, k_1) + np.diag(b, k_2) + np.diag(c, k_3)
# it solves a tridiagonal linear system - i.e. form below:
# b[0]*x[0] + c[0]*x[1] = d[0]
# a[0]*x[0] + b[1]*x[1] + c[1]*x[2] = d[1]
# . . . . .
# . . . . .
# . . . . .
# . . . . .
# . . . . .
# . . . . .
# a[n-1]*x[n-1] + b[n] *x[n] + c[n] *x[n+1] = d[n]
# a[n] *x[n] + b[n+1]*x[n+1] = d[n+1]
def solve_tridiagonal_matrix(a, b, c, d):
b = b.copy()
d = d.copy()
n = len(b)
# we subtract a multiple of the first eq. from the second eq. so that x_0 disappears
# then we subtract a multiple of the second eq. from the first (new) eq. so that x_1 disappears and so on...
# keep track of the coefficients of the new reuslting matrix
# e.g. if we have:
# 0) 2*x0 + 3*x1 = 2
# 1) 4*x0 + 8*x1 + 2*x2 = 5
# 2) 6*x1 + 7*x2 + 1*x3 = 16
# 3) 1*x2 + 2*x3 = 17
# then:
# eq. 0) is unchanged, new eq.1) is obtained subtracting 2*eq.0) from eq.1); new eq.2) is obtained subtraing 3*new_eq.1) from old eq.2) and so on...
# 0-new) 2*x0 + 3*x1 = 2 # unchanged
# 1-new) 2*x1 + 2*x2 = 1 # subtract 2*eq.0 from eq.1 (old)
# 2-new) 1*x2 + 1*x3 = 1 # subtract 3*eq.1-new from eq.2 (old)
# 3-new) 1*x3 = 1 # subtract 1*eq.2-new from eq.3 (old)
# note that at the end of this process:
# i) new array a (below diagonal) becomes 0 (so no calcs needed)
# ii) array c is unchanged (so no calcs needed), while new b, d needs to be calculated iteratively
# iii) the last variable (x3) can be calculated directly, and the previous ones recoursively backward
# (i.e. get x3 from eq.3-new, substitute x3 in eq.2-new and get x_2, and so on...)
# at this point we can solve for x, using iii) above
for i in range(1, n):
# calculate the ratio, w, between the coefficient of x[i-1] in eq. i-1) and the coefficient of x[i-1] in eq. i)
w = a[i-1] / b[i-1]
# new eq.i) is obtained as: eq.i (old) - w * eq.i-1 (new)
# by construction new c[i] = 0, while b[i] and d[i] are calculated:
b[i] = b[i] - w * c[i-1]
d[i] = d[i] - w * d[i-1]
# we obtained the new system in which the lower diagonal is zeroed. We can recoursively calculate the solution x
x = np.zeros(n, dtype = float)
# last (new) equation is b[n-1]*x[n-1] = d[n-1]
# --> we can solve for x[n-1]
x[n-1] = d[n-1] / b[n-1]
for i in range(n-2, -1, -1):
# recoursively we solve each quation, starting from the second-to-last, moving backward to the first
# at step i, the i-th equation reads: b[i]*x[i] + c[i]*x[i+1] = d[i]
# x[i+1] being iteratively calculated at previous step, we can solve for x[i]
x[i] = (d[i] - c[i]*x[i+1])/b[i]
return x
# it solves a linear system of the form:
# (I call it a crank-nickolson form beacause this algorithm involves the solution of a linear system of this form)
# p[0]*x[0] + p[1]*x[1] + p[2]*x[2] + ................ + p[m1-1]*x[m1-1] = d[0]
# a[0]*x[0] + b[0]*x[1] + c[0]*x[2] = d[1]
# a[1]*x[1] + b[1]*x[2] + c[1]*x[3] = d[2]
# . . . . .
# . . . . .
# . . . . .
# . . . . .
# . . . . .
# . . . . .
# . . . . .
# . . . . .
# . . . . .
# . . . . .
# . . . . .
# a[n-1]*x[n-1] + b[n]*x[n] + c[n] *x[n+1] . .
# a[n]*x[n-1] + b[n]*x[n] + c[n]*x[n+1] = d[n]
# q[0]*x[n-m2] + q[1]*x[n-m2+1] + ............. + q[m1-2]*x[n] + q[m2-1]*x[n+1] = d[n+1]
# the below function converts the linear system into a tridiagonal linear system
# it returns the arrays (A, B, C) defining the tridiagonal matrix and the new known vector (D)
def get_trid_arrays(p, q, a, b, c, d):
# e.g. if we have:
# 0) 2*x0 + 3*x1 + 5*x2 + 4*x3 + 5*x4 + 5*x5 = 10
# 1) 1*x0 + 2*x1 + 1*x2 = 0
# 2) 2*x1 + 2*x2 + 3*x3 = -4
# 3) 2*x2 + 1*x3 + 1*x4 = 1
# 4) 1*x3 + 3*x4 + 5*x5 = 2
# 5) 4*x4 + 2*x5 + 3*x6 = 24
# 6) 1*x5 + 4*x6 + 6*x7 = 32
# 7) 1*x2 + 3*x3 + 2*x4 + 6*x5 + 7*x6 + 2*x7 = 13
# STEP 1) subtract eq.4) from eq.0) --> new eq.0) becomes:
# eq.0-new, step 1) 2*x0 + 3*x1 + 5*x2 + 3*x3 + 2*x4 = 8 # x5 disappeared!
# STEP 2) subtract 2*eq.4) from eq. 1-new (step 1) --> new eq.1) becomes:
# eq.0-new, step 2) 2*x0 + 3*x1 + 1*x2 + 1*x3 = 6 # x4 disappeared!
# and so on... until we only x0 and x1 appears in the first equation
# we do the same thing for eq.8, first removing a multiple of eq.2 (x1 disappears), then a multiple of eq.3 (x2 disappears)...
# ... and so on, until only x6 and x7 remains
# the resulting system has a tridiagonal form
# note that intermediate (1 to 6) equations are unaffected
# only 1) and 8) change, so we need only keep track of p, q, d[0], d[last]
n = len(a)
# warning! it is expected that n = len(a) = len(b) = len(c), while len(d) = n+2
m1 = len(p)
p = p.copy()
d = d.copy()
for i in range(m1-1, 1, -1): # i=m1-1,i=m2-2,...,i=2
# in eq.0), we need to make disappear first x[m1-1], then x[m1-2],..., and last x[2]
# (in the example above, p has lenght 6 (from x0 to x5), so we make disappear first x5, then x4, then x3, and last x2
# iteratively, to make disappear x[i], new eq.0 is obtained subtracting a multiple eq. i-1 from previous eq. 0
# (check in the example: e.g. i=5, to make disappear x5, we need to subtract eq.4 (=5-1) from eq. 0)
# the multiple is given by the ratio of the coefficient of x[i] between eq. 0) and eq. i-1)
# in eq. 0), the coefficient of x[i] is p[i] (e.g. in eq.0 coeff of x5 is p[5]
# in eq. i-1), the coefficient of x[i] is c[i-2] (e.g. in eq.4 coeff of x5 is c[3])
# the other two variables in eq. i-1) are x[i-1] and x[i-2], with coefficients b[i-2] and a[i-2], respectively
# (check in the example: e.g. i=5 the other two variables in eq. 4 are x4 and x3, with coefficients b[3] and a[3], respectively
w = p[i] / c[i-2]
p[i] = 0
p[i-1] = p[i-1] - w * b[i-2]
p[i-2] = p[i-2] - w * a[i-2]
d[0] = d[0] - w * d[i-1]
# we need to do something similar for last equation
m2 = len(q)
q = q.copy()
for i in range(m2-2): #i=0,i=1,...,i=m2-3
# the # of variables (= # of equations) is len(d) = n+2 - n inner equations, then the first and last equations
# in last equation, eq. n+1, we need to make disappear first x[n-m2+2], then x[n-m2+2+1],..., x[n-m2+2+i], ..., and last x[n-1] - i ranges from 0 to m2-3
# (check: in the example above, n=6 (eq. 1 to 6), q has lenght m2=6 (from x2 to x7), so we make disappear first x[2] (1=6-6+2), then x[2+1], ..., and last x5=x[2+3],
# 3 = m2-3 = 6+3, so that i ranges from 0 to m2-n-3
# now, the i-th variable of last equation is x[n-m+2+i] and has coefficient of q[i]
# to make disappear this variable we need to subtract a multiple of eq. n-m2+3 from last equation
# (check: in the example above, if i=0, to remove x2 (2=n-m2+2+i=6-6+2+0), we need to subtract eq. 3 (3=n-m2+3=6-3+3) from last equation
# the coefficient of this variable in this eq. is a[n-m2+2+i]
# (check: in the example above, if i=0, the coefficient of x2 in eq.3 is a[2] (2=n-m2+2+i=6-6+2+0))
# (check2: in the example above, if i=1, the coefficient of x3 in eq.4 is a[3] (3=n-m2+2+i=6-6+2+1))
# (check3: if instead last equation had been 2*x4 + 6*x5 + 7*x6 + 2*x7 = 13, then m2 = 4, to remove x4 (4=n-m2+2=6-4+2), we need to subtract eq. 5 from last eq.
# the coeff of x4 in eq. 5 (i=0) is a[4] (4 = n-m2+2+i = 6-4+2). In total, we need to remove x4, x5=x[4+1], so i ranges from 0 to 1 (1=m2-3= 4-3). OK
w = q[i] / a[n-m2+i+2]
q[i] = 0
q[i+1] = q[i+1] - w * b[n-m2+i+2] # n-m2+2+i = n-1 --> i = m2-3
q[i+2] = q[i+2] - w * c[n-m2+i+2]
d[n+1] = d[n+1] - w * d[n-m2+i+3] # d[j] from eq. j, j=n-m2+3 (see comments above)
# now, the resulting tridiagonal matrix is built with a,b,c unchanged and the new p,q:
# |p[0] p[1] |
# |a[0] b[0] c[0] |
# | a[1] b[1] c[1] |
# | . . . |
# | . . . |
# | . . . |
# | a[n-1] b[n-1] c[n-1] |
# | a[n] b[n] c[n] |
# | q[m2-2] q[m2-1] |
# we return below diagonal, diagonal, and upper diagonal arrays
# A: below diagonal
# A = [a[0], a[1], ... , a[n-1], a[n], q[m2-2]]
A = np.zeros(n+1, dtype=float)
A[0:n] = a
A[n] = q[m2-2]
# B: diagonal
# B = [p[0], b[0], b[1], ... , b[n-1], b[n], q[m2-1]]
B = np.zeros(n+2, dtype=float)
B[1:n+1] = b
B[0] = p[0]
B[n+1] = q[m2-1]
# C: upper diagonal
# C = [p[1], c[0], c[1], ... , c[n-1], c[n]]
C = np.zeros(n+1, dtype=float)
C[1:n+1] = c
C[0] = p[1]
D = np.array(d, dtype=float)
return [A, B, C, D]
def solve_crank_nickolson_linear_system(p, q, a, b, c, d):
# convert into tridiagonal form - get trid arrays (A,B,C) and new known vector (D)
[A, B, C, D] = get_trid_arrays(p, q, a, b, c, d)
# solve the new tridiagonal system
x = solve_tridiagonal_matrix(A, B, C, D)
return x
#a = [1, 1, 1, 1, 1, 1]
#b = [1, 2, 3, 4, 5, 6]
#c = [6, 5, 4, 3, 2, 1]
#d = [1, 1, 1, 1, 1, 1, 1, 1]
##dd = np.array([1, 1, 1, 1, 1, 1, 1, 1], dtype = float)
#p = [1, -3, 3, -1]
#q = [1, -3, 3, -1]
#[A, B, C, D] = get_trid_arrays(p, q, a, b, c, d)
#n = len(d)
#m1 = len(p)
#m2 = len(q)
def get_crank_nickolson_matrix(p, q, a, b, c):
n = len(a) +2
m1 = len(p)
m2 = len(q)
mat = np.zeros([n,n], dtype = float)
for i in range(m1):
mat[0, i] = p[i]
for i in range(m2):
mat[n-1, n-m2+i] = q[i]
for i in range(0, n-2):
mat[i+1, i] = a[i]
mat[i+1, i+1] = b[i]
mat[i+1, i+2] = c[i]
return mat
#mat2 = tridiag_matrix(A, B, C)
#np.linalg.det(mat2)
#x = solve_tridiagonal_matrix(A, B, C, D)
#mat = get_crank_nickolson_matrix(p, q, a, b, c)
#check = np.matmul(mat, np.transpose(x))
#################################
def my_fct(N_x, N_t, forw_or_backw = "forward", x_bound_cond = "exact_solution", method_if_backward = "backward_equations"):
# a * y_xx + b * y_x + c * y - y_t = 0
a = 1
b = 1
c = -1
# solution: y = exp(x+t)
x_min = -4
x_max = 4
T_min = 0
T_max = 2
k = (T_max - T_min)/N_t # = d_t
d_t = k
h = (x_max - x_min)/N_x # = d_x
d_x = h
A_next_val = -(2*k*a - k*h*b)
A_prev_val = (2*k*a - k*h*b)
if forw_or_backw == "forward":
h_sq_sign = 1
else:
h_sq_sign = -1
B_next_val = (4*pow(h,2)*h_sq_sign + 4*k*a - 2*pow(h,2)*k*c)
B_prev_val = (4*pow(h,2)*h_sq_sign - 4*k*a + 2*pow(h,2)*k*c)
C_next_val = -(2*k*a + k*h*b)
C_prev_val = (2*k*a + k*h*b)
A_next = np.array([A_next_val for i in range(1, N_x)], dtype=float)
A_prev = np.array([A_prev_val for i in range(1, N_x)], dtype=float)
B_next = np.array([B_next_val for i in range(1, N_x)], dtype=float)
B_prev = np.array([B_prev_val for i in range(1, N_x)], dtype=float)
C_next = np.array([C_next_val for i in range(1, N_x)], dtype=float)
C_prev = np.array([C_prev_val for i in range(1, N_x)], dtype=float)
y = np.zeros([(N_x+1), (N_t+1)])
x = np.array([(x_min + i*d_x) for i in range(N_x+1)], dtype=float)
if forw_or_backw == "forward":
y_T_min = np.array([exp(x_i+T_min) for x_i in x], dtype=float)
y_prev = y_T_min
y[:, 0] = y_T_min
y_tilde = None
else:
y_T_max = np.array([exp(x_i+T_max) for x_i in x], dtype=float)
if method_if_backward == "y_tilde":
y_tilde = np.zeros([(N_x + 1), (N_t + 1)])
y_tilde_prev = y_T_max
y_tilde[:, 0] = y_T_max
elif method_if_backward == "backward_equations":
y_next = y_T_max
y_tilde = None
y[:, N_t] = y_T_max
if x_bound_cond == "exact_solution":
p = [1, 0]
q = [0, 1]
elif x_bound_cond == "third_order_smoothing":
p = [1, -3, 3, -1]
q = [1, -3, 3, -1]
for j in range(1, N_t+1): # j = N_t-1, N_t-2,..., 1, 0
if forw_or_backw == "forward":
D = np.zeros(N_x + 1, dtype=float)
D[1: N_x] = np.array(
[(A_prev[i - 1] * y_prev[i - 1] + B_prev[i - 1] * y_prev[i] + C_prev[i - 1] * y_prev[i + 1]) for i in
range(1, N_x)],
dtype=float) # i-1 = 0, 1,..., N_x-3, N_x-2
if x_bound_cond == "exact_solution":
D[0] = exp(x_min + (T_min + j * d_t))
D[N_x] = exp(x_max + (T_min + j * d_t))
y_next = solve_crank_nickolson_linear_system(p, q, A_next, B_next, C_next, D)
y[:, j] = y_next
y_prev = y_next
else:
D = np.zeros(N_x + 1, dtype=float)
if method_if_backward == "y_tilde":
D[1: N_x] = np.array([(A_prev[i-1] * y_tilde_prev[i - 1] + B_prev[i-1] * y_tilde_prev[i] + C_prev[i-1] * y_tilde_prev[i + 1]) for i in range(1, N_x)],
dtype=float) # i-1 = 0, 1,..., N_x-3, N_x-2
t_j = N_t - j
if x_bound_cond == "exact_solution":
D[0] = exp(x_min + (T_min + t_j * d_t))
D[N_x] = exp(x_max + (T_min + t_j * d_t))
y_tilde_next = solve_crank_nickolson_linear_system(p, q, A_next, B_next, C_next, D)
y_tilde[:, j] = y_tilde_next
y[:, t_j] = y_tilde_next
y_tilde_prev = y_tilde_next
elif method_if_backward == "backward_equations":
D[1: N_x] = np.array([(A_next[i-1] * y_next[i - 1] + B_next[i-1] * y_next[i] + C_prev[i-1] * y_next[i + 1]) for i in range(1, N_x)],
dtype=float) # i-1 = 0, 1,..., N_x-3, N_x-2
t_j = N_t - j
if x_bound_cond == "exact_solution":
D[0] = exp(x_min + (T_min + t_j * d_t))
D[N_x] = exp(x_max + (T_min + t_j * d_t))
y_prev = solve_crank_nickolson_linear_system(p, q, A_prev, B_next, C_next, D)
y[:, t_j] = y_prev
y_next = y_prev
return [y, y_tilde]
def solve_my_crank(multipliers, forw_or_backw = "forward", x_bound_cond = "exact_solution", method_if_backward = "backward_equations"):
print("forward or backward:", forw_or_backw, ", x_bound_cond:", x_bound_cond, ", method_if_backward:", method_if_backward)
# print("N.B. N_x = 10 * multiplier, N_t = 4 * multiplier")
for multiplier in multipliers:
N_x = 10 * multiplier
N_t = 4 * pow(multiplier, 2)
x_min = -4
x_max = 4
T_min = 0
T_max = 2
k = (T_max - T_min) / N_t # = d_t
d_t = k
h = (x_max - x_min) / N_x # = d_x
d_x = h
[y, y_tilde] = my_fct(N_x, N_t, forw_or_backw, x_bound_cond, method_if_backward)
# i = int(N_x/2); j = N_t; x_i = x_min+i*d_x; t_j = T_min+j*d_t; print("i: ", i, ", j: ", j, " --> exp(x_i+t_j): ", exp(x_i+t_j), ", y[i, j]): ", y[i, j])
mean_err = np.mean([np.mean([abs(exp(x_min+i*d_x + T_min+j*d_t) - y[i,j]) for j in range(N_t+1)]) for i in range(N_x+1)])
mean_abs_errs.append(mean_err)
print("multiplier:", multiplier, ", N_x:", N_x, ", N_t:", N_t, ", k:", "{:.3f}".format(k), ", h^2:", "{:.4f}".format(pow(h,2)), ", k/h^2: ", "{:.3f}".format(k/pow(h,2)), " --> mean abs err is:", "{:.4f}".format(mean_err))
print()
print("-----------------------------------")
multipliers = [1, 3, 5, 10, 30]
mean_abs_errs = []
forw_or_backw = "forward"
solve_my_crank(multipliers, forw_or_backw)
forw_or_backw = "forward"
x_bound_cond = "third_order_smoothing"
solve_my_crank(multipliers, forw_or_backw, x_bound_cond)
forw_or_backw = "backward"
x_bound_cond = "exact_solution"
method_if_backward = "backward_equations"
solve_my_crank(multipliers, forw_or_backw, x_bound_cond, method_if_backward)
forw_or_backw = "backward"
x_bound_cond = "exact_solution"
method_if_backward = "y_tilde"
solve_my_crank(multipliers, forw_or_backw, x_bound_cond, method_if_backward)
```Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.