Was reading up on solving non-linear PDE's, and came across Newton's method as an approach.

Newton's method for solving non-linear equations,

f(x) = 0

involves first starting with a guess x_0 and then improving that guess using

x_1=x_0-\frac{f(x_0)}{f'(x_0)}

This comes from Taylor expanding the function f about x_0

f(x_1) = f(x_0) + f'(x_0) (x_1-x_0) = 0

and then solving for x_1

For a system of coupled equations, we have

\vec{f}(\mathbf{x}_1) = \vec{f}(\mathbf{x}_0) + D\vec{f}(\mathbf{x}_0) (\mathbf{x}_1-\mathbf{x}_0) = 0

where D\vec{f} is a matrix of partial derivatives

To get an better guess, we need to first solve

D\vec{f}(\mathbf{x}_0) \Delta \mathbf{x} = -\vec{f}(\mathbf{x}_0)

and then we can update

\mathbf{x}_1=\mathbf{x}_0+\Delta \mathbf{x}

In terms of the diffusion equation, we have a non-linear system of equations for the new temperature T'

\Xi_{\vec{k}}(\vec{T'}) = 0

where \vec{k} is a vector index that maps to a scalar index in the matrix equation.

Newton's method would say we update our new value of T' using

\vec{T''}=\vec{T'}+\vec{\Delta T'}

where

D_{\vec{k}}\Xi_{\vec{k'}}(\vec{T'}) \Delta T'_{\vec{k}} = -\Xi_{\vec{k'}}(\vec{T'})

What would \Xi(T') and \Xi(T') look like for anisotropic diffusion?

Our starting point for solving the non-linear equation implicitly is the following

\partial_t T = \displaystyle \sum_{\parallel, \perp} A \partial_i B_{ij} \partial_j T_*^{\lambda+1}

where T_* = T + \phi ( T' - T ) = (1-\phi) T + \phi T' = \psi T + \phi T'

we can identify then that

\Xi(T') = \partial_t T - \displaystyle \sum_{\parallel, \perp} A \partial_i B_{ij} \partial_j \left ( \psi T + \phi T' \right ) ^{\lambda+1}

or in discretized form

\Xi_{\vec{k}}(\vec{T'}) = \frac{T_{\vec{k}}' - T_{\vec{k}}}{\Delta t} - \displaystyle \sum_{\parallel, \perp} \left [ \alpha_0 \left ( \psi T_{\vec{k}} + \phi T'_{\vec{k}} \right ) ^{\lambda+1} + \displaystyle \sum_{\pm j} \alpha_{\pm j} \left (\psi T_{\vec{k}\pm \hat{j}} + \phi T'_{\vec{k}\pm \hat{j}} \right )^{\lambda+1} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} \left ( \psi T_{\vec{k}\pm \hat{i} \pm \hat{j}} + \phi T'_{\vec{k}\pm \hat{i} \pm \hat{j}} \right )^{\lambda+1} \right ]

and

D_{\vec{k}}\Xi_{\vec{k'}}(T') = \frac{\delta_{\vec{k}}^{\vec{k'}}}{\Delta t} - \phi \left ( \lambda + 1 \right ) \displaystyle \sum_{\parallel, \perp} \left [ \alpha_0 \delta^{\vec{k'}}_{\vec{k}}\left ( \psi T_{\vec{k'}} + \phi T'_{\vec{k'}} \right ) ^{\lambda} + \displaystyle \sum_{\pm j} \alpha_{\pm j} \delta^{\vec{k'}\pm \hat{j}}_{ \vec{k}} \left (\psi T_{\vec{k'}\pm \hat{j}} + \phi T'_{\vec{k'}\pm \hat{j}} \right )^{\lambda} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} \delta^{\vec{k'}\pm \hat{i} \pm \hat{j}}_{\vec{k}} \left ( \psi T_{\vec{k'}\pm \hat{i} \pm \hat{j}} + \phi T'_{\vec{k'}\pm \hat{i} \pm \hat{j}} \right )^{\lambda} \right ]

This makes the matrix equation for \Delta T'_{\vec{k}}

\left [ \frac{\delta_{\vec{k}}^{\vec{k'}}}{\Delta t} - \displaystyle \sum_{\parallel, \perp} \phi \left ( \lambda + 1 \right ) \left [ \alpha_0 \delta^{\vec{k'}}_{\vec{k}}\left ( \psi T_{\vec{k'}} + \phi T'_{\vec{k'}} \right ) ^{\lambda} + \displaystyle \sum_{\pm j} \alpha_{\pm j} \delta^{\vec{k'}\pm \hat{j}}_{ \vec{k}} \left (\psi T_{\vec{k'}\pm \hat{j}} + \phi T'_{\vec{k'}\pm \hat{j}} \right )^{\lambda} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} \delta^{\vec{k'}\pm \hat{i} \pm \hat{j}}_{\vec{k}} \left ( \psi T_{\vec{k'}\pm \hat{i} \pm \hat{j}} + \phi T'_{\vec{k'}\pm \hat{i} \pm \hat{j}} \right )^{\lambda} \right ] \right ] \Delta T'_{\vec{k}} = -\Xi_{\vec{k'}}(\vec{T'})

and summing over \vec{k}

\frac{\Delta T'_{\vec{k'}}}{\Delta t} - \displaystyle \sum_{\parallel, \perp} \phi \left ( \lambda + 1 \right ) \left [ \alpha_0 \left ( \psi T_{\vec{k'}} + \phi T'_{\vec{k'}} \right ) ^{\lambda} \Delta T'_{\vec{k'}} + \displaystyle \sum_{\pm j} \alpha_{\pm j} \left (\psi T_{\vec{k'}\pm \hat{j}} + \phi T'_{\vec{k'}\pm \hat{j}} \right )^{\lambda} \Delta T'_{\vec{k'} \pm \hat{j}} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} \left ( \psi T_{\vec{k'}\pm \hat{i} \pm \hat{j}} + \phi T'_{\vec{k'}\pm \hat{i} \pm \hat{j}} \right )^{\lambda} \Delta T'_{\vec{k'}\pm \hat{i} \pm \hat{j}} \right ] =-\Xi_{\vec{k'}}(\vec{T'})

or

\Delta T'_{\vec{k}}- \Delta t \phi \left ( \lambda + 1 \right ) \displaystyle \sum_{\parallel, \perp} \left [ \alpha_0 \left ( \psi T_{\vec{k}} + \phi T'_{\vec{k}} \right ) ^{\lambda} \Delta T'_{\vec{k}} + \displaystyle \sum_{\pm j} \alpha_{\pm j} \left (\psi T_{\vec{k}\pm \hat{j}} + \phi T'_{\vec{k}\pm \hat{j}} \right )^{\lambda} \Delta T'_{\vec{k} \pm \hat{j}} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} \left ( \psi T_{\vec{k}\pm \hat{i} \pm \hat{j}} + \phi T'_{\vec{k}\pm \hat{i} \pm \hat{j}} \right )^{\lambda} \Delta T'_{\vec{k}\pm \hat{i} \pm \hat{j}} \right ]

= -T'_{\vec{k}} + T_{\vec{k}} + \Delta t \displaystyle \sum_{\parallel, \perp} \left [ \alpha_0 \left ( \psi T_{\vec{k}} + \phi T'_{\vec{k}} \right ) ^{\lambda+1} + \displaystyle \sum_{\pm j} \alpha_{\pm j} \left (\psi T_{\vec{k}\pm \hat{j}} + \phi T'_{\vec{k}\pm \hat{j}} \right )^{\lambda+1} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} \left ( \psi T_{\vec{k}\pm \hat{i} \pm \hat{j}} + \phi T'_{\vec{k}\pm \hat{i} \pm \hat{j}} \right )^{\lambda+1} \right ]

or we can write it as

T''_{\vec{k}} = T'_{\vec{k}} + \Delta T'_{\vec{k}} = T_{\vec{k}} + \Delta t \displaystyle \sum_{\parallel, \perp} \left [ \alpha_0 \left ( \psi T_{\vec{k}} + \phi T'_{\vec{k}} \right ) ^{\lambda+1} + \displaystyle \sum_{\pm j} \alpha_{\pm j} \left (\psi T_{\vec{k}\pm \hat{j}} + \phi T'_{\vec{k}\pm \hat{j}} \right )^{\lambda+1} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} \left ( \psi T_{\vec{k}\pm \hat{i} \pm \hat{j}} + \phi T'_{\vec{k}\pm \hat{i} \pm \hat{j}} \right )^{\lambda+1} \right ] +\Delta t \displaystyle \sum_{\parallel, \perp} \phi \left ( \lambda + 1 \right ) \left [ \alpha_0 \left ( \psi T_{\vec{k}} + \phi T'_{\vec{k}} \right ) ^{\lambda} \Delta T'_{\vec{k}} + \displaystyle \sum_{\pm j} \alpha_{\pm j} \left (\psi T_{\vec{k}\pm \hat{j}} + \phi T'_{\vec{k}\pm \hat{j}} \right )^{\lambda} \Delta T'_{\vec{k} \pm \hat{j}} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} \left ( \psi T_{\vec{k}\pm \hat{i} \pm \hat{j}} + \phi T'_{\vec{k}\pm \hat{i} \pm \hat{j}} \right )^{\lambda} \Delta T'_{\vec{k}\pm \hat{i} \pm \hat{j}} \right ]

Now if we start with T'=T and set \Delta T' = T''-T we end up with

T''_{\vec{k}} = T_{\vec{k}} + \Delta t \displaystyle \sum_{\parallel, \perp} \left [ \alpha_0 T_{\vec{k}}^{\lambda+1} + \displaystyle \sum_{\pm j} \alpha_{\pm j}T_{\vec{k}\pm \hat{j}}^{\lambda+1} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T_{\vec{k}\pm \hat{i} \pm \hat{j}}^{\lambda+1} \right ] +\Delta t \displaystyle \sum_{\parallel, \perp} \phi \left ( \lambda + 1 \right ) \left [ \alpha_0 T_{\vec{k}}^{\lambda} \left(T''_{\vec{k}} - T_{\vec{k}} \right ) + \displaystyle \sum_{\pm j} \alpha_{\pm j} T_{\vec{k}\pm \hat{j}}^{\lambda} \left(T''_{\vec{k}\pm \hat{j}} - T_{\vec{k}\pm \hat{j}} \right ) + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T_{\vec{k}\pm \hat{i} \pm \hat{j}}^{\lambda} \left ( T''_{\vec{k}\pm \hat{i} \pm \hat{j}} - T_{\vec{k}\pm \hat{i} \pm \hat{j}} \right ) \right ]

which can be rearranged as

T''_{\vec{k}} - \Delta t \displaystyle \sum_{\parallel, \perp} \phi \left ( \lambda + 1 \right ) \left [ \alpha_0 T_{\vec{k}}^{\lambda} T''_{\vec{k}} + \displaystyle \sum_{\pm j} \alpha_{\pm j} T_{\vec{k}\pm \hat{j}}^{\lambda} T''_{\vec{k}\pm \hat{j}} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T_{\vec{k}\pm \hat{i} \pm \hat{j}}^{\lambda}T''_{\vec{k}\pm \hat{i} \pm \hat{j}} \right ] = T_{\vec{k}} + \Delta t \displaystyle \sum_{\parallel, \perp} \left ( 1 - \phi \left ( \lambda + 1 \right ) \right ) \left [ \alpha_0 T_{\vec{k}}^{\lambda+1} + \displaystyle \sum_{\pm j} \alpha_{\pm j}T_{\vec{k}\pm \hat{j}}^{\lambda+1} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T_{\vec{k}\pm \hat{i} \pm \hat{j}}^{\lambda+1} \right ]

which is the linearized implicit equation we are solving. Not surprisingly, a single-step newton method is what we are doing by Taylor expanding T*^{\lambda+1}.