---
title: 4. Allen-Cahn equation
tags: Finite difference method
---

---

### Problem 1 (閔中)

@mjgu327 再幾個問題: 1. 原題中一切都沒有單位, 所以你另外加上 sec, m 等單位看起來很怪. 2. 原題在 $t=4$, $34$ 以及 $186$ 三個時間點的解究竟長什麼樣子? 在影片中我看到三條不重合的線, 所以解究竟是哪一條? 3. 我還是沒有看到關於時間 $O(\Delta t)$ 的收斂性結果. 關於空間, 需要將使用的 $\Delta x$ 等資訊補上, 而且一般我們會需要將 $\Delta x$ 切細二到三次並計算 R 值, 這樣才知道是否 $R\to 0.25$.

@mjgu327 兩個問題: 1. 看起來有空間方向 dx 的收斂性, 時間上有收斂性嗎?
2. 後半部 R 值隨時間的演進, 算這個有什麼特別的意義嗎? 或是想要表達什麼?

@teshenglin 老師我改好了，再麻煩您批改，謝謝您。

Allen-Cahn equation
$$
u_t = 0.01u_{xx}+u-u^3, \quad x\in[-4, 4], \quad t\in[0, 500],
$$
with initial and boundary conditions
$$
u(x,0) = -1 + 2\left(e^{-35(x+2)^2} + e^{-11 x^2} + e^{-7 (x-2)^2}\right), \quad u(-4, t) = u(4,t)=-1.
$$
Sol:

**(a)**
**(b)**
**(c)**
**(d)**

%Matalb code

clc
clear

format long

L = 8;
T = 500;

Nx1 = 1000;    
dx1 = L/Nx1;
x1 = (-4:dx1:4)';
U1 = -1+2.*(exp(-35*(x1+2).^2)+exp(-11*x1.^2)+exp(-7*(x1-2).^2));
U1(1,1)  = -1;
U1(Nx1,1) = -1;

Nx2 = 2000;    
dx2 = L/Nx2;
x2 = (-4:dx2:4)';
U2 = -1+2.*(exp(-35*(x2+2).^2)+exp(-11*x2.^2)+exp(-7*(x2-2).^2));
U2(1,1)  = -1;
U2(Nx2,1) = -1;

Nx3 = 4000;    
dx3 = L/Nx3;
x3 = (-4:dx3:4)';
U3 = -1+2.*(exp(-35*(x3+2).^2)+exp(-11*x3.^2)+exp(-7*(x3-2).^2));
U3(1,1)  = -1;
U3(Nx3,1) = -1;

dt = 200*dx3^2/(4+100*abs(max((3*U3.^2))-1)*dx3^2);

t = (1:dt:T)';
R_1 = zeros(size(t));
R_2 = zeros(size(t));
R_inf = zeros(size(t));

for  n = 1:length(t);



U1(2:Nx1,1) = U1(2:Nx1,1) + 0.01*dt/dx1^2.*(U1(1:Nx1-1,1)-2*U1(2:Nx1,1)+U1(3:Nx1+1,1)) + ...
                dt*(U1(2:Nx1,1) - U1(2:Nx1,1).^3);
U2(2:Nx2,1) = U2(2:Nx2,1) + 0.01*dt/dx2^2.*(U2(1:Nx2-1,1)-2*U2(2:Nx2,1)+U2(3:Nx2+1,1)) + ...
                dt*(U2(2:Nx2,1) - U2(2:Nx2,1).^3);
U3(2:Nx3,1) = U3(2:Nx3,1) + 0.01*dt/dx3^2.*(U3(1:Nx3-1,1)-2*U3(2:Nx3,1)+U3(3:Nx3+1,1)) + ...
                dt*(U3(2:Nx3,1) - U3(2:Nx3,1).^3);       
            
R_1(n)   = norm(U2(1:2:end)-U3(1:4:end),1)/norm(U1(:)-U2(1:2:end),1);
R_2(n)   = norm(U2(1:2:end)-U3(1:4:end),2)/norm(U1(:)-U2(1:2:end),2);
R_inf(n) = norm(U2(1:2:end)-U3(1:4:end),inf)/norm(U1(:)-U2(1:2:end),inf);   


 if mod(n,1000) == 0;
    figure(1)
    plot(x1,U1,'r'); hold on;
    plot(x2,U2,'b'); hold on;
    plot(x3,U3,'g'); hold off;
    title(['time = ',num2str((n)*dt)])
    axis([-4 4 -2 2])
    xlabel('x')
    ylabel('U')
    legend('U|_N_=_1_0_0_0','U|_N_=_2_0_0_0','U|_N_=_4_0_0_0')

  end
end

figure(2)
plot(t,R_1,'r'); hold on;
plot(t,R_2,'b'); hold on;
plot(t,R_inf,'g'); hold off;
axis([0 T 0 1])
xlabel('t')
ylabel('R')
legend('R_1','R_2','R_\infty')

**(a)**
**(b)**
**(c)**
**(d)**

%Matalb code
clc
clear

format long

L  = 8;
Nx = 2000;
dx = L/Nx;
x  = (-4:dx:4)';
T = 500;

U1 = -1+2.*(exp(-35*(x+2).^2)+exp(-11*x.^2)+exp(-7*(x-2).^2));
U1(1,1)  = -1;
U1(Nx,1) = -1;
U2 = U1;
U3 = U1;

dt = 200*dx^2/(4+100*abs(max((3*U1.^2))-1)*dx^2);

R_1 = zeros(fix(T/dt),1);
R_2 = zeros(fix(T/dt),1);
R_inf = zeros(fix(T/dt),1);

for  n = 1:fix(T/dt);


U1(2:Nx,1) = U1(2:Nx,1) + 0.01*dt/dx^2.*(U1(1:Nx-1,1)-2*U1(2:Nx,1)+U1(3:Nx+1,1)) + ...
                dt*(U1(2:Nx,1) - U1(2:Nx,1).^3);

for i = 1:2
U2(2:Nx,1) = U2(2:Nx,1) + 0.01*(dt/2)/dx^2.*(U2(1:Nx-1,1)-2*U2(2:Nx,1)+U2(3:Nx+1,1)) + ...
                (dt/2)*(U2(2:Nx,1) - U2(2:Nx,1).^3);
end

for j = 1:4
U3(2:Nx,1) = U3(2:Nx,1) + 0.01*(dt/4)/dx^2.*(U3(1:Nx-1,1)-2*U3(2:Nx,1)+U3(3:Nx+1,1)) + ...
                (dt/4)*(U3(2:Nx,1) - U3(2:Nx,1).^3);
end            

R_1(n)   = norm(U2-U3,1)/norm(U1-U2,1);
R_2(n)   = norm(U2-U3,2)/norm(U1-U2,2);
R_inf(n) = norm(U2-U3,inf)/norm(U1-U2,inf);   

 if mod(n,200) == 0;
    figure(1)
    plot(x,U1,'r'); hold on;
    plot(x,U2,'b'); hold on;
    plot(x,U3,'g'); hold off;
    title(['time = ',num2str((n+0.4)*dt)])
    axis([-4 4 -2 2])
    xlabel('x')
    ylabel('U')
    legend('U|dt_1','U|dt_2','U|dt_3')
 end
 
end

---

### Problem 2 (睿珩)

@NKFUST 1. 證明你做的是對的, 2. 所以 t=10000 解長怎樣?

Allen-Cahn equation
$$
u_t = 0.01u_{xx}+u-u^3, \quad x\in[-4, 4], \quad t\in[0, 10000],
$$
with initial and boundary conditions
$$
u(x,0) = -1 + 10\left(e^{-35(x+2)^2} + e^{-11 x^2} + e^{-7 (x-2)^2}\right), \quad u(-4, t) = u(4,t)=-1.
$$

Solution:

Divide $x$ into $N$ equal parts.$U(-4,t)=U_{0}^{m}=-1$, $U(4,t)=U_{N}^{m}=-1$.

forward difference:
$$
\frac{\partial^2{U_{n}^{m}}}{\partial{x^2}}=\frac{U_{n-1}^{m}-2U_{n}^{m}+U_{n+1}^{m}}{\Delta{x^2}}+O(\Delta{x^2})
$$
$$
\frac{\partial{U_{n}^{m}}}{\partial{t}}=\frac{U_{n}^{m+1}-U_{n}^{m}}{\Delta{t}}+O(\Delta{t})
$$
Algorithm formula:
$$
U_{n}^{m+1}=U_{n}^{m}+\frac{0.01}{\Delta{x^2}}(U_{n-1}^{m}-2U_{n}^{m}+U_{n+1}^{m})+\Delta{t}[U_{n}^{m}-(U_{n}^{m})^3]+O(\Delta{t},\Delta{x^2})
$$

Matlab Code:

clear;

L = 8;
Nx = 1000;
T = 10000;
Nt = 10000000;

dx = L/Nx;
dt = T/Nt;
x = [-4:dx:4]';

U = -1+10*(exp(-35*(x+2).^2)+exp(-11*x.^2)+exp(-7*(x-2).^2));
U(1,1)  = -1;
U(Nx,1) = -1;

for  n = 1:Nt
    U(2:Nx,1) = U(2:Nx,1) + 0.01*dt/dx^2.*(U(1:Nx-1,1)-2*U(2:Nx,1)+U(3:Nx+1,1)) + ...
                dt*(U(2:Nx,1) - U(2:Nx,1).^3);
    figure(1);
    plot(x,U,'r');
    title([num2str(n*dt*1000),' milliseconds']);
    axis([-4 4 -10 10]);
    xlabel('x');
    ylabel('U');
end

Number of x meshes is 1000, number of t meshes is 10000000, t = 10000 :

Convergence of t. Number of t meshes in U1 is 10000000, U2 is 20000000 and U4 is 40000000.
Algorithm validation in convergence of t.

clear;

L = 8;
T = 10000;

Nx = 1000;
Nt = 10000000;

dx = L/Nx;
dt = T/Nt;
x = [-4:dx:4]';

U = -1+10*(exp(-35*(x+2).^2)+exp(-11*x.^2)+exp(-7*(x-2).^2));
U(1,1)  = -1;
U(Nx,1) = -1;
%
Nx2 = 1000;
Nt2 = 20000000;

dx2 = L/Nx2;
dt2 = T/Nt2;
x2 = [-4:dx2:4]';

U2 = -1+10*(exp(-35*(x2+2).^2)+exp(-11*x2.^2)+exp(-7*(x2-2).^2));
U2(1,1)  = -1;
U2(Nx2,1) = -1;
%
Nx4 = 1000;
Nt4 = 40000000;

dx4 = L/Nx4;
dt4 = T/Nt4;
x4 = [-4:dx4:4]';

U4 = -1+10*(exp(-35*(x4+2).^2)+exp(-11*x4.^2)+exp(-7*(x4-2).^2));
U4(1,1)  = -1;
U4(Nx4,1) = -1;

for  n = 1:Nt
    U(2:Nx,1) = U(2:Nx,1) + 0.01*dt/dx^2.*(U(1:Nx-1,1)-2*U(2:Nx,1)+U(3:Nx+1,1)) + ...
                dt*(U(2:Nx,1) - U(2:Nx,1).^3);
    
    U2(2:Nx2,1) = U2(2:Nx2,1) + 0.01*dt2/dx2^2.*(U2(1:Nx2-1,1)-2*U2(2:Nx2,1)+U2(3:Nx2+1,1)) + ...
                dt2*(U2(2:Nx2,1) - U2(2:Nx2,1).^3);
    
    U4(2:Nx4,1) = U4(2:Nx4,1) + 0.01*dt4/dx4^2.*(U4(1:Nx4-1,1)-2*U4(2:Nx4,1)+U4(3:Nx4+1,1)) + ...
                dt4*(U4(2:Nx4,1) - U4(2:Nx4,1).^3);
    
    figure(1);
    plot(x,U,'r',x2,U2,'g',x4,U4,'b');
    title([num2str(n*dt*1000),' milliseconds']);
    axis([-4 4 -10 10]);
    xlabel('x');
    ylabel('U');
    legend({'U1','U2','U4'})
end

Convergence of x. Number of x meshes in U1 is 100, U2 is 200 and U4 is 400.
Algorithm validation in convergence of x.

clear;

L = 8;
T = 10000;

Nx = 100;
Nt = 10000000;

dx = L/Nx;
dt = T/Nt;
x = [-4:dx:4]';

U = -1+10*(exp(-35*(x+2).^2)+exp(-11*x.^2)+exp(-7*(x-2).^2));
U(1,1)  = -1;
U(Nx,1) = -1;
%
Nx2 = 200;
Nt2 = 10000000;

dx2 = L/Nx2;
dt2 = T/Nt2;
x2 = [-4:dx2:4]';

U2 = -1+10*(exp(-35*(x2+2).^2)+exp(-11*x2.^2)+exp(-7*(x2-2).^2));
U2(1,1)  = -1;
U2(Nx2,1) = -1;
%
Nx4 = 400;
Nt4 = 10000000;

dx4 = L/Nx4;
dt4 = T/Nt4;
x4 = [-4:dx4:4]';

U4 = -1+10*(exp(-35*(x4+2).^2)+exp(-11*x4.^2)+exp(-7*(x4-2).^2));
U4(1,1)  = -1;
U4(Nx4,1) = -1;

for  n = 1:Nt
    U(2:Nx,1) = U(2:Nx,1) + 0.01*dt/dx^2.*(U(1:Nx-1,1)-2*U(2:Nx,1)+U(3:Nx+1,1)) + ...
                dt*(U(2:Nx,1) - U(2:Nx,1).^3);
    
    U2(2:Nx2,1) = U2(2:Nx2,1) + 0.01*dt2/dx2^2.*(U2(1:Nx2-1,1)-2*U2(2:Nx2,1)+U2(3:Nx2+1,1)) + ...
                dt2*(U2(2:Nx2,1) - U2(2:Nx2,1).^3);
    
    U4(2:Nx4,1) = U4(2:Nx4,1) + 0.01*dt4/dx4^2.*(U4(1:Nx4-1,1)-2*U4(2:Nx4,1)+U4(3:Nx4+1,1)) + ...
                dt4*(U4(2:Nx4,1) - U4(2:Nx4,1).^3);
    
    figure(1);
    plot(x,U,'r',x2,U2,'g',x4,U4,'b');
    title([num2str(n*dt*1000),' milliseconds']);
    axis([-4 4 -10 10]);
    xlabel('x');
    ylabel('U');
    legend({'U1','U2','U4'});
end

---