Assumed I have the following two coupled equations $$\begin{split} \partial_tA&=a_0AB\\ \partial_tB&=b_0AB \end{split} $$ but I am not sure how to calculate them. One approach is a crank-nicolson approach: $$\begin{split} \left(1-\frac{dt\cdot a_0B_0}{2}\right)A_1&=\left(1+\frac{dt\cdot a_0B_0}{2}\right)A_0\\ \left(1-\frac{dt\cdot b_0A_0}{2}\right)B_1&=\left(1+\frac{dt\cdot b_0A_0}{2}\right)B_0 \end{split} $$ or a modified crank-nicolson approach, using the result already obtained: $$\begin{split} \left(1-\frac{dt\cdot a_0B_0}{2}\right)A_1&=\left(1+\frac{dt\cdot a_0B_0}{2}\right)A_0\\ \left(1-\frac{dt\cdot b_0A_1}{2}\right)B_1&=\left(1+\frac{dt\cdot b_0A_1}{2}\right)B_0 \end{split} $$ Finally, I tested the usual forward euler method: $$\begin{split} A_1&=A_0+dt\cdot a_0AB\\ B_1&=B_0+dt\cdot b_0AB \end{split} $$ When choosing $a_0=1$, $b_0=2$ and $A_0=B_0=1$, all three methods result in different results, according to my matlab script:

clear all;
x_max = 1;
x_num = 1000;
a0_val = 1;
b0_val = 2;
dx = x_max/20000;
steps = 1/dx;

%euler equations
A_equation = @(A_vec, B_vec) a0_val*A_vec.*B_vec;
B_equation = @(A_vec, B_vec) b0_val*A_vec.*B_vec;

A_matrix_A = @(B_vec) speye(x_num)-speye(x_num).*a0_val*dx/2.*B_vec;
B_matrix_A = @(B_vec) speye(x_num)+speye(x_num).*a0_val*dx/2.*B_vec;

A_matrix_B = @(A_vec) speye(x_num)-speye(x_num).*b0_val*dx/2.*A_vec;
B_matrix_B = @(A_vec) speye(x_num)+speye(x_num).*b0_val*dx/2.*A_vec;

A_matrix = zeros(steps, x_num);
B_matrix = zeros(steps, x_num);
A_matrix_old = zeros(steps, x_num);
B_matrix_old = zeros(steps, x_num);
A_matrix_euler = zeros(steps, x_num);
B_matrix_euler = zeros(steps, x_num);

B0_vec = ones(x_num, 1);
A0_vec = ones(x_num, 1);
A_matrix(1,:) = A0_vec;
B_matrix(1,:) = B0_vec;
A_matrix_old(1,:) = A0_vec;
B_matrix_old(1,:) = B0_vec;
A_matrix_euler(1,:) = A0_vec;
B_matrix_euler(1,:) = B0_vec;

%Do steps
for i = 2:steps
    B_matrix(i,:) = A_matrix_B(A_matrix(i-1,:))\(B_matrix_A(A_matrix(i-1,:))*B_matrix(i-1,:).'+A_matrix(i-1,:).');
    A_matrix(i,:) = A_matrix_A(B_matrix(i,:))\(B_matrix_A(B_matrix(i,:))*A_matrix(i-1,:).'+B_matrix(i,:).');

for i = 2:steps
    B_matrix_old(i,:) = A_matrix_B(A_matrix_old(i-1,:))\(B_matrix_A(A_matrix_old(i-1,:))*B_matrix_old(i-1,:).'+A_matrix_old(i-1,:).');
    A_matrix_old(i,:) = A_matrix_A(B_matrix_old(i-1,:))\(B_matrix_A(B_matrix_old(i-1,:))*A_matrix_old(i-1,:).'+B_matrix_old(i-1,:).');

for i = 2:steps
    A_matrix_euler(i,:) = A_matrix_euler(i-1,:)+dx*A_equation(A_matrix_euler(i-1,:), B_matrix_euler(i-1,:));
    B_matrix_euler(i,:) = B_matrix_euler(i-1,:)+dx*B_equation(A_matrix_euler(i-1,:), B_matrix_euler(i-1,:));


Thus I was wondering which method is correct in this situation?


1 Answer 1


You should reformulate your problem. Let's define the vector $u$ as $u=\left(\begin{array}{c}A\\ B\end{array}\right)$ Then you can write your coupled system as $$\frac{\partial u}{\partial t}=\left(\begin{array}{c}a_0\\ b_0\end{array}\right)\left(u^T \left(\begin{array}{cc}0&0.5\\0.5&0\end{array}\right) u\right)$$

Now we can apply Crank Nicolson as defined in the wikipedia article:

$$\frac{u^{n+1}-u^n}{\Delta t}=\frac{1}{2}\left(\left(\begin{array}{c}a_0\\ b_0\end{array}\right)\left({u^{n+1}}^T \left(\begin{array}{cc}0&0.5\\0.5&0\end{array}\right) u^{n+1}\right)+\left(\begin{array}{c}a_0\\ b_0\end{array}\right)\left({u^n}^T \left(\begin{array}{cc}0&0.5\\0.5&0\end{array}\right) u^n\right)\right)$$

As one can see, this is a nonlinear function in $u$. So for solving this step of the Crank Nicolson scheme one should use an iterative scheme like Newton-Kantorovich with the Frechet derivative of the expression.


Be $M=\left(\begin{array}{cc}0&0.5\\0.5&0\end{array}\right)$ and $c=\left(\begin{array}{c}a_0\\b_0\end{array}\right)$. To solve the nonlinear equation for each time step, we apply the Newton-Kantorovich iteration: $$N_u\Delta=-N(u)$$, here $N_u\Delta$ is the Frechet derivative of $$N(u)=\frac{u-u^n}{\Delta t}-\frac{1}{2}c\left(\left({u}^T M u\right)+\left({u^n}^T M u^n\right)\right)$$. I will use the greek index $\alpha$ to illustrate the fix point iteration. The Frechet derivative is defined as $$\lim_{\epsilon\rightarrow 0} \frac{N(u+\epsilon\Delta) -N(u)}{\epsilon}$$, hence $$N_u(u)\Delta=\left(\left(\begin{array}{cc}\frac{1}{\Delta t}&0\\0&\frac{1}{\Delta t}\end{array}\right)-c\left(u^T M\right)\right) \Delta$$.

Now for each time step, you do the fix point iteration starting with $u^{0,n+1}=u^n$ by solving $$N_u(u^{\alpha,n+1})\Delta=-N(u^{\alpha,n+1})$$ for $\Delta$, till $\Delta$ is small and updating $$u^{\alpha+1,n+1}=u^{\alpha,n+1}+\Delta$$ for each iteration step.

When the fix point iteration is converged (e.g.$|\Delta|<\epsilon$), you have found the next time step $u^{n+1}$ you can use in your Crank Nicolson scheme.


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.