# Solving a PDE implicitly by iteration in python

Connected to this question here on Computational Science, I've posted a follow-up question on how to solve a PDE using an implicit scheme like Crank-Nicholson in general in this question on SO.
But I guess it is better to ask here it is correct to solve a PDE with an implicit step in this way in general and ask the more specific coding issues on SO. So I'll try to split up the question I asked on SO.

So I've got a non-linear advection-diffusion(-reaction) PDE which I want to solve for one single timestep using Crank-Nicolson. In this example my 1D-PDE is the following:

$$\frac{T^{n+1} - T^n}{\Delta t} \approx \frac{1}{2} \left(\left(\alpha(T^n) \frac{T^n_{i-1} - 2T^n_{i} + T^n_{i+1}}{\Delta x^2} + v(T^n) \frac{T^n_{i\pm1} - T^n_{i}}{\Delta x}\right) + \left( \alpha(T^{n+1}) \frac{T^{n+1}_{i-1} - 2T^{n+1}_{i} + T^{n+1}_{i+1}}{\Delta x^2} + v(T^{n+1}) \frac{T^{n+1}_{i\pm1} - T^{n+1}_{i}}{\Delta x}\right)\right)$$

Where $\alpha$ is the temperature dependent diffusion coefficient and $v$ is the temperature dependent advection speed. The advection is discretized using an upwind scheme, since when there is advection $|Pe|\geq2$ can be assumed.

My current program code returns the derivatives of the PDE for given input temperatures, summing up the specific advection and diffusion derivatives. Thus my solver can't "see" the diffusion coefficient and the advection speed. Using a matrix notation to solve the PDE is thus not applicable.
Furthermore the derivatives and value (temperature) arrays have the same shape as the calculation grid, meaning that for the simplest case of the 1D PDE with, for example 10 nodes, the derivatives and value arrays are vectors with 10 cells (of course neglecting the additional ghost cells I use for fast stencil solving and boundary conditions).

The derivative terms in my PDE are thus in my program code replaced with the following vector derivative where each element in $derivative$ corresponds to the derivative of its respective temperature element in the value array: $$derivative(T^n) = \alpha(T^n) \frac{T^n_{i-1} - 2T^n_{i} + T^n_{i+1}}{\Delta x^2} + v(T^n) \frac{T^n_{i\pm1} - T^n_{i}}{\Delta x}$$

So in pseudo-code equation my PDE looks like this: $$\frac{T^{n+1} - T^n}{\Delta t} \approx \frac{1}{2} \left(derivative(T^{n}) + derivative(T^{n+1})\right)$$

To solve this, I reorganize this equation so that $$T^{n+1} - T^n - \frac{\Delta t}{2} \left(derivative\left(T^{n}\right) + derivative(T^{n+1})\right) = 0$$

which in python looks like

def crank_nicolson(y, yprev, h):
return (y - yprev - h / 2 * (diff(y) + diff(yprev)))


with the stepsize $h = \Delta t$, $y=T^{n+1}$, $yprev=T^{n}$ and $diff = derivative(T)$.
This can easily be fed to any root finding algorithm like the scipy.optimize.root module or using any other self-written iterative solver for example like:

from scipy import optimize

def diff(T):  # simply example differential function
return np.sqrt(T)

def iterate_cranknicolson(y0, h, rtol, diff):  # simple iterative solver
y = y0
diff_old = diff(y)  # save forward diff
f = np.ones_like(y)
i = 0
while np.any(f > rtol):
y_old = y
y = y0 + h / 2 * (diff_old + diff(y))  # at the first step: euler forward predictor
f = (y / y_old - 1)  # simple stop criterion
i += 1
return y

# some test runs
T0 = np.arange(10, 20).astype(float)  # set starting temperature

cn_hybr = optimize.root(crank_nicolson, T0, args=(T0, 1.), method='hybr', tol=1e-6)
cn_it = iterate_cranknicolson(T0, 1., 1e-6, diff)


## TL,DR. The question(s) summmarized:

• Is this approach to iteratively solve an implicit discretization of a PDE correct? With "this approach" I mean especially reorganizing the PDE to $T^{n+1} - T^n - \frac{\Delta t}{2} \left(derivative\left(T^{n}\right) + derivative(T^{n+1})\right) = 0$ and finding the root of this equation.
• Is there any better way to do it without being able to use matrix notation?
• Is the derivative of $T^n$ the first column of the Jacobian matrix, so that I could use a modified Newtonian iteration, as proposed here and does this work without having a matrix notation at all?
• (not really important, thus secondary: Is it ok to use upwind scheme for advection in Crank-Nicolson or should I stick to central differences?)

Thanks alot in advance! Please tell me if I should refine parts of my question and/or add more detail.

To answer your questions: Is there any better way to do it without being able to use matrix notation? Yes, there are alternative ways, but this is not the best approach. In other words, this method (explained below) work perfectly in simplified homogeneous case using fully implicit scheme than Crank-Nicolson and the way you do it is by taking all the T_i^(n+1) terms from the right-hand side to the left-hand side

(2/dx^2 + 1/dx^2 + 1)T_i^(n+1) = dt[.... rest part of the right-hand sid]

make (2/dx^2 + 1/dx^2 + 1) = K then T_i^(n+1) = dt*[....]/K

Finally solve this system inside loop iteration just as you solving explicit method

*** notice that you need to give a very good initial condition for this method to work very well You can read about some useful ways to give good initial conditions for example you solve the system for one time-step for very small dt ( has some criteria dt < some function) then. then run your iteration and make sure you update T_i^(n) = T_i^(n+1) for the next time step

1. Is the derivative of Tn the first column of the Jacobian matrix, so that I could use a modified Newtonian iteration, as proposed here and does this work without having a matrix notation at all?

The system you are solving is a linear problem; you really do not need Jacobian to use NR to solve the problem. NR only for non-linear systems. However, there are NR solver for linear problem you can use that one, but like I said, you dont need

Finally the matrix approach is the more practical way to solve CN scheme