Crank-Nicolson algorithm for coupled PDEs

Question

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,:).');
end

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,:).');
end

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,:));
end

figure;
imagesc(B_matrix);
figure;
imagesc(B_matrix_old);
figure;
imagesc(B_matrix_euler);

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

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.

Edit:

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.