How to Change Domain

Been working on writing a Finite-Difference scheme for approximating solutions to the Allen-Cahn equation:

u_t = (epsilon)^2*(u_{xx}+u_{yy}) + f(u) where the reaction term f(u) is given as u-u^3.

After lots of patience and searching, I have found the following MATLAB routine for a very similar equation:

    % ------------------------------------------------------------------ %
    % Pattern formation:  a 15 line matlab program                       %
    % PDE form: u_t=D*(u_{xx}+u_{yy})+gamma*q(u) where q(u)='u.*(1-u)';  %
    % ------------------------------------------------------------------ %
    % This KPP equation and its numerical solution are discussed         %
    % in detail in Section 10.4, Chapter 10 of the following book by     %
    %    Xin-She Yang, Introduction to Computational Mathematics,        % 
    %    World Scientific, (2008)                                        %
    %    http://www.worldscibooks.com/mathematics/6867.html              %
    % ------------------------------------------------------------------ %

    % ================================================================== %
    % The solution of this PDE is obtained by the finite difference      %
    % method, assuming dx=dy=dt=1.                                       % 
    % Programmed by X S Yang (Cambridge University) @ 2008               %
    % ================================================================== %
    % Usage: pattern(200)   or  simply >pattern
    % ================================================================== %

    % ----------------------------------------------------
    function rdpde(n)                    % line 1
    % Default the domain size n by n
    if nargin<1, n=200; end                % line 2
    % ----------------------------------------------------
    % Initialize parameters
    % ---- time=100, D=0.2; gamma=0.5; -------------------
    % These values can be changed, but be careful, and 
    % some stabilities analysis is recommended to identify 
    % the right parameter values for pattern formation
    time=100; D=0.001; gamma=1;            % line 3 %D=0.2, gamma=0.5
    % ---- Set initial values of u randomly --------------
    u=rand(n,n);  grad=u*0;                % line 4
    % Vectorization/index for u(i,j) and the loop --------
    I = 2:n-1; J = 2:n-1;                  % line 5
    % ---- Time stepping ---------------------------------
    for step=1:time,                       % line 6
    % Laplace gradient of the equation     % line 7
     grad(I,J)= u(I,J-1)+u(I,J+1)+u(I-1,J)+u(I+1,J);
     u =(1-4*D^2)*u+D^2*grad+gamma*u.*(1-u).*(1+u);   % line 8      %D^2 ??
    % ----- Show results ---------------------------------
     pcolor(u);  shading interp;           % line 9
    % ----- Coloring and showing colorbar ----------------
     colorbar; colormap hsv;               % line 10
     drawnow;                              % line 11
    end                                    % line 12
    % ----- Topology of the final surface ----------------
    surf(u);                               % line 13
    shading interp; colormap jet;          % line 14
    view([-25 70]);                        % line 15
    % ---- end of this program ---------------------------

I have made some modifications to the original script, namely, the 'D' and 'gamma' terms, but I am still unsure if this is the proper way to approach this. Any feedback and/or suggestions on how to modify/generalize this code will be greatly appreciated.

Best Regards,

Jamil Villarreal

Answer 1

[Alan Weiss]

I believe that this equation falls easily into the domain of Partial Differential Equation Toolbox. You seem to have two space dimensions, which the toolbox handles. I did not see any geometry, initial conditions, or boundary conditions, which you would have to specify for PDE Toolbox. But the coefficients are easily put into the PDE Toolbox framework:

d = 1;
c = epsilon^2;
f = @ffun;

Here you need to write a function for the f coefficient whose syntax depends on your toolbox version. For many toolbox versions, you could simply use the string 'u-u.^3' as your f coefficient instead of writing a function.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Jamil Villarreal]

Hello Alan Weiss, and thank you for your feedback.

I have heard about the pdetoolbox before, but have never actually worked with it, so would you mind further explaining how to implement this Allen-Cahn type equation into this framework?

Also, I realized that this program starts from a random state and then shows the changes occurring. I have just written a separate code that uses a Crank-Nicolson approach for the spacial derivatives and a Forward-Difference scheme for time. If possible, I could post it as a separate question since I am having issues with the dimensions of the matrices involved.

Best,

Jamil Villarreal

[Alan Weiss]

The documentation for the toolbox is here. I suggest that you especially look at the types of PDE problems you can solve and how to specify coefficients in string form.

If you then believe that the toolbox might be able to solve your problem, look at the entire set of steps you would need to take and see if you can provide the requisite information. If so, then you would need to purchase a license for the toolbox (unless you already have one; check by issuing the ver command at the MATLAB command line).

FYI, PDE Toolbox uses a finite element method for solving PDEs numerically, not a finite difference scheme.

Alan Weiss

MATLAB mathematical toolbox documentation