I appreciate this might be an easy question, but I've managed to get myself quite thoroughly confused
I'm trying to solve a system of physics equations that look as follows
$$ \frac{\partial \mathbf{E}}{\partial t}=-\alpha \mathbf{j} \\ \frac{\partial \mathbf{j}}{\partial t}=\beta\mathbf{E}+\gamma\mathbf{j}\times\mathbf{B} $$ $\alpha, \beta, \gamma$ are scalar constants and $\mathbf{B}$ is a constant vector with 3 components. $\mathbf{E}, \mathbf{j}$ also both have 3 components.
I would like to solve these equations implicitly. I can easily work out the 1st order implicit Euler method, it ends up with something like $$ \begin{pmatrix} \mathbf{E} \\ \mathbf{j} \end{pmatrix}^{n}=\underline{\mathbf{M}}\cdot\begin{pmatrix} \mathbf{E} \\ \mathbf{j} \end{pmatrix}^{n+1} $$ where $\mathbf{M}$ is just some 6x6 matrix I can invert through some linear algebra routine. However, I am finding that this implicit Euler method for solving this problem isn't sufficiently accurate and I would like a second order method. I just can't work out how to apply a second order method to this problem.
I've been thinking about using Crank-Nicholson, which for $\frac{\partial y}{\partial t} = f(t, y)$ can be written like
$$ y_{n+1}=y_n + \frac{h}{2}(k_1 + k_2) \\ k_1 = f(t_n, y_n) \\ k_2 = f(t_n + h, y_n + \frac{h}{2}(k_1 + k_2)) $$ But I can't get my head around how you solve the last step for $k_2$ for my problem when it isn't so much a function of $f(t, y)$ as much as a sort of matrix product like
$$ \frac{\partial}{\partial t}\begin{pmatrix} \mathbf{E} \\ \mathbf{j} \\ \end{pmatrix} = \underline{\mathbf{M}}\cdot\begin{pmatrix} \mathbf{E} \\ \mathbf{j} \end{pmatrix} $$ I don't have a good conception of what $k_1$ or $k_2$ actually look like for my problem.

I think they should look something like $$ \mathbf{k_1} = \underline{\mathbf{M}}\cdot\mathbf{y}^n \\ \mathbf{k_2} = \underline{\mathbf{M}}\cdot(\mathbf{y}^n + \frac{h}{2}\underline{\mathbf{M}}\cdot\mathbf{y}^n + \frac{h}{2}\mathbf{k_2}) $$ where $\mathbf{y}$ is the vector of $\mathbf{E}, \mathbf{j}$. But is this remotely correct?


1 Answer 1


In the linear case you get $$k_1+k_2=2My+\frac h2 M(k_1+k_2),$$ so that isolating the sum you get $$k_1+k_2=2(I-hM/2)^{-1}My.$$ For the step this gives $$ y^{n+1}=y^n+\frac h2(k_1+k_2)=y^n+(I-hM/2)^{-1}hMy^n \\ =(I-hM/2)^{-1}[(I-hM/2)+hM]y^n=(I-hM/2)^{-1}(I+hM/2)y^n $$ as one also gets directly from the trapezoidal formula $$ \frac{y^{n+1}-y^n}h=\frac{My^n+My^{n+1}}2 $$

Practically one could most easily solve for the mean slope $k_m=\frac{k_1+k_2}2$, as then $$ k_m=\frac{f(t^n,y^n)+f(t^{n+1},y^n+hk_m)}2. $$ If $J$ is a sufficiently accurate approximation of the Jacobian around $(t^n,y^n)$, then $f(t^{n+1},y^n+hk_m)-hJk_m$ will be nearly flat as function of $k_m$ and small $h$, so that the fixed-point iteration $$ k_m^{new}=(I-hJ/2)^{-1}\frac{f(t^n,y^n)+f(t^{n+1},y^n+hk_m^{old})-hJk_m^{old}}2 $$ Here "flat" implies a small Lipschitz constant of $O(h^2)$ for this iteration. This is usually sufficient to get good enough approximation in 1-3 iterations, as the error in the fixed-point equation decreases by a factor $O(h^2)$ in each step. So starting with $k_m=k_1$, after one step the error is $O(h^4)$ and thus smaller than the $O(h^3)$ error for the method step.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.