# Deriving coefficients for the implicit scheme
---------------------------------------------------

The ice sheet energy balance model uses an implicit scheme to solve the
heat equation for $N$ layers. It uses the Crank-Nicholson scheme to discretise
the equations. For an equation in one space dimension $x$

$\frac{df}{dt} = F$,

the Crank-Nicholson scheme discretises the equation as

$\frac{f^{i+1}_x - f^{i}_x}{\Delta t} = 0.5\left [ F^{i+1}_x + F^{i}_x \right ]$

where the superscript is time and the subscript is space.


```python
from sympy import *
init_printing()
```


```python
tnew_x = Symbol('T^{i+1}_x')
tnew_xprev = Symbol('T^{i+1}_{x-1}')
tnew_xafter = Symbol('T^{i+1}_{x+1}')

told_x = Symbol('T^{i}_x')
told_xprev = Symbol('T^{i}_{x-1}')
told_xafter = Symbol('T^{i}_{x+1}')

u_x = Symbol('\kappa_x')
u_xprev = Symbol('\kappa_{x-1}')
u_xafter = Symbol('\kappa_{x+1}')

delta_t = Symbol('\Delta t')
delta_x = Symbol('\Delta x')
```


```python
told_x, u_xprev, tnew_xafter, delta_x
```


```python
lhs = (tnew_x - told_x)/delta_t
```


```python
lhs # The time derivative
```


```python
rhs_new = 0.5*(u_x*(tnew_xprev - 2*tnew_x + tnew_xafter)/delta_x**2 + 
               ((tnew_x - tnew_xprev)/(delta_x))*((u_x - u_xprev)/(delta_x)))
```


```python
rhs_old = 0.5*(u_x*(told_xprev - 2*told_x + told_xafter)/delta_x**2 + 
               ((told_x - told_xprev)/(delta_x))*((u_x - u_xprev)/(delta_x)))
```


```python
rhs_new, rhs_old # The two parts of the crank-nicholson RHS.
```


```python
expr = lhs - rhs_new - rhs_old
```


```python
expr
```


```python
poly_form = Poly(expr, tnew_x, tnew_xafter, tnew_xprev, told_x, told_xafter, told_xprev)
```


```python
poly_form
```

The coefficients for the $i+1$ temperature (predicted) are


```python
(poly_form.coeff_monomial(tnew_xprev)*delta_t).simplify(), (poly_form.coeff_monomial(tnew_x)*delta_t).simplify(), (poly_form.coeff_monomial(tnew_xafter)*delta_t).simplify()
```

The coefficients for the $i$ temperature (current) are


```python
-(poly_form.coeff_monomial(told_xprev)*delta_t).simplify(), (poly_form.coeff_monomial(told_x)*-delta_t).simplify(), -(poly_form.coeff_monomial(told_xafter)*delta_t).simplify()
```

These are the coefficients for the diagonals of the sparse matrix. The need to be multiplied by $\frac{1}{\rho C}$ where $\rho, C$ are the vertically varying density and specific heat to represent snow and ice. 

The $\kappa$ values are staggered from the $T$ values, since $T$ represents the interface temperatures and $\kappa$ represents the mid_level values. `K_mid` and `K_bar` in the code therefore interpolate $\kappa$ values to the interface.

bbbf80b08cf1159d74c48907c9f5ddccaec47ab9
ipynb
Jupyter Notebook
docs/Implicit_Scheme_Derivation.ipynb
Ai33L/climt
17fcdc43ec8ed9903dfb566bae99154f282ce85f
BSD-3-Clause
2017-02-12T01:45:50.000Z
2022-03-01T20:19:03.000Z
docs/Implicit_Scheme_Derivation.ipynb
Ai33L/climt
17fcdc43ec8ed9903dfb566bae99154f282ce85f
BSD-3-Clause
2017-07-23T21:19:15.000Z
2022-02-25T16:49:21.000Z
docs/Implicit_Scheme_Derivation.ipynb
Ai33L/climt
17fcdc43ec8ed9903dfb566bae99154f282ce85f
BSD-3-Clause
2017-05-24T11:52:12.000Z
2021-12-09T16:45:19.000Z
Qwen/Qwen-72B
1. YES
2. YES