Hi everyone,
I want to numerically solve the following nonlinear two dimensional reaction-diffusion equations by discretising the spatial derivatives using finite difference so that I end up with a set of odes that can be solved using an ode solver -- method of lines.
I have had a go at coding the problem in Matlab but got unrealistic/unexpected results. Can someone please have a look at my code and shed some light? I would appreciate any assistance as I am still learning my way around Matlab. Many thanks.

The three equations read:
du_1/dt=d^2u_1/dx^2 + d^2u_1/dy^2 + u_1u_2-u_1+u_2u_3;
du_2/dt=d^2u_2/dx^2 + d^2u_2/dy^2 - u_1u_2+u_1u_3-1u_2;
du_3/dt=u_1u_2-u_1+u_3u_1;

That's what I have put in the ode file:

function out1 = myode(t,u,flag)

%%parameters
N=50;
L=50;
N2=N*N;
d=2L/N-1;
d2=d^2;
x=-L:d:L;
y=-L:d:L;

%%define the three variables
u1=sparse(N2,N2);
u2=sparse(N2,N2);
u3=sparse(N2,N2);

%%the diffusion term
u1diff=-delsq(numgrid('S',N+2));
u2diff=-delsq(numgrid('S',N+2));

%%set up the boundary terms (fixed at all edges)
u1(1:N2,1)=1;
u1(1:N2,N2)=1;
u1(1,1:N2)=1;
 
u2(1:N2,1)=.98;
u2(1:N2,N2)=.98;
u2(1,1:N2)=.98;
u2(N2,1:N2)=.98;

%% then put the above equations as 
u1eqn= u1diff + the rest of the terms;
u2eqn= u2diff + the rest of the terms;
u3eqn= the RHS;

out11=[u1eqn; u2eqn; u3eqn]';
out1=out11(:);

In the run file I have:

%% set up initial conditions
u1ics=1.*ones(N2,N2);
u2ics=1.*ones(N2,N2);
u3ics=.9.*ones(N2,N2);
u4ics=.9.*ones(N2,N2);

loop over time
    [t,u] = ode15s('spatial_ph_2d_ode',tspan,ics);
    ics = u(end,:)';
    u1= reshape(ics(1:N2),N,N);
    u2= reshape(ics(N2+1:2*N2),N,N);
    u3= reshape(ics(2*N2+1:3*N2),N,N);
pcolor(x,y,u1)
pcolor(x,y,u2)
pcolor(x,y,u3)

end

> Hi everyone,
> I want to numerically solve the following nonlinear
> two dimensional reaction-diffusion equations by
> discretising the spatial derivatives using finite
> difference so that I end up with a set of odes that
> can be solved using an ode solver -- method of lines.
> I have had a go at coding the problem in Matlab but
> got unrealistic/unexpected results. Can someone
> please have a look at my code and shed some light? I
> would appreciate any assistance as I am still
> learning my way around Matlab. Many thanks.
> 
> The three equations read:
> du_1/dt=d^2u_1/dx^2 + d^2u_1/dy^2 +
> u_1u_2-u_1+u_2u_3;
> du_2/dt=d^2u_2/dx^2 + d^2u_2/dy^2 -
> u_1u_2+u_1u_3-1u_2;
> du_3/dt=u_1u_2-u_1+u_3u_1;
> 
> That's what I have put in the ode file:
> 
> function out1 = myode(t,u,flag)
> 
> %%parameters
> N=50;
> L=50;
> N2=N*N;
> d=2L/N-1;
> d2=d^2;
> x=-L:d:L;
> y=-L:d:L;
> 
> %%define the three variables
> u1=sparse(N2,N2);
> u2=sparse(N2,N2);
> u3=sparse(N2,N2);
> 
> %%the diffusion term
> u1diff=-delsq(numgrid('S',N+2));
> u2diff=-delsq(numgrid('S',N+2));
> 
> %%set up the boundary terms (fixed at all edges)
> u1(1:N2,1)=1;
> u1(1:N2,N2)=1;
> u1(1,1:N2)=1;
>  
> u2(1:N2,1)=.98;
> u2(1:N2,N2)=.98;
> u2(1,1:N2)=.98;
> u2(N2,1:N2)=.98;
> 
> %% then put the above equations as 
> u1eqn= u1diff + the rest of the terms;
> u2eqn= u2diff + the rest of the terms;
> u3eqn= the RHS;
> 
> out11=[u1eqn; u2eqn; u3eqn]';
> out1=out11(:);
> 
> In the run file I have:
> 
> %% set up initial conditions
> u1ics=1.*ones(N2,N2);
> u2ics=1.*ones(N2,N2);
> u3ics=.9.*ones(N2,N2);
> u4ics=.9.*ones(N2,N2);
> 
> loop over time
>     [t,u] = ode15s('spatial_ph_2d_ode',tspan,ics);
>     ics = u(end,:)';
>     u1= reshape(ics(1:N2),N,N);
>     u2= reshape(ics(N2+1:2*N2),N,N);
>     u3= reshape(ics(2*N2+1:3*N2),N,N);
> pcolor(x,y,u1)
> pcolor(x,y,u2)
> pcolor(x,y,u3)
> 
> end

Hello Emma,

I miss the essential routine 'spatial_ph_2d_ode'.
If you think the above routine 'myode' returns
the time derivatives of u in the grid points, it can't
be true. In the calculation of the array 'out1', you
do not use 'u' to get the diffusion terms u1diff and 
u2diff.

Best wishes
Torsten.

ok thanxs for sharing.

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <994710478.5017.1295335738754.JavaMail.root@gallium.mathforum.org>...
> > Hi everyone,
> > I want to numerically solve the following nonlinear
> > two dimensional reaction-diffusion equations by
> > discretising the spatial derivatives using finite
> > difference so that I end up with a set of odes that
> > can be solved using an ode solver -- method of lines.
> > I have had a go at coding the problem in Matlab but
> > got unrealistic/unexpected results. Can someone
> > please have a look at my code and shed some light? I
> > would appreciate any assistance as I am still
> > learning my way around Matlab. Many thanks.
> > 
> > The three equations read:
> > du_1/dt=d^2u_1/dx^2 + d^2u_1/dy^2 +
> > u_1u_2-u_1+u_2u_3;
> > du_2/dt=d^2u_2/dx^2 + d^2u_2/dy^2 -
> > u_1u_2+u_1u_3-1u_2;
> > du_3/dt=u_1u_2-u_1+u_3u_1;
> > 
> > That's what I have put in the ode file:
> > 
> > function out1 = myode(t,u,flag)
> > 
> > %%parameters
> > N=50;
> > L=50;
> > N2=N*N;
> > d=2L/N-1;
> > d2=d^2;
> > x=-L:d:L;
> > y=-L:d:L;
> > 
> > %%define the three variables
> > u1=sparse(N2,N2);
> > u2=sparse(N2,N2);
> > u3=sparse(N2,N2);
> > 
> > %%the diffusion term
> > u1diff=-delsq(numgrid('S',N+2));
> > u2diff=-delsq(numgrid('S',N+2));
> > 
> > %%set up the boundary terms (fixed at all edges)
> > u1(1:N2,1)=1;
> > u1(1:N2,N2)=1;
> > u1(1,1:N2)=1;
> >  
> > u2(1:N2,1)=.98;
> > u2(1:N2,N2)=.98;
> > u2(1,1:N2)=.98;
> > u2(N2,1:N2)=.98;
> > 
> > %% then put the above equations as 
> > u1eqn= u1diff + the rest of the terms;
> > u2eqn= u2diff + the rest of the terms;
> > u3eqn= the RHS;
> > 
> > out11=[u1eqn; u2eqn; u3eqn]';
> > out1=out11(:);
> > 
> > In the run file I have:
> > 
> > %% set up initial conditions
> > u1ics=1.*ones(N2,N2);
> > u2ics=1.*ones(N2,N2);
> > u3ics=.9.*ones(N2,N2);
> > u4ics=.9.*ones(N2,N2);
> > 
> > loop over time
> >     [t,u] = ode15s('myode',tspan,ics);
> >     ics = u(end,:)';
> >     u1= reshape(ics(1:N2),N,N);
> >     u2= reshape(ics(N2+1:2*N2),N,N);
> >     u3= reshape(ics(2*N2+1:3*N2),N,N);
> > pcolor(x,y,u1)
> > pcolor(x,y,u2)
> > pcolor(x,y,u3)
> > 
> > end
> 
> Hello Emma,
>
> I miss the essential routine 'spatial_ph_2d_ode'.

Thanks Torsten for replying to my message.
sorry, the run file should read 
 [t,u] = ode15s('myode',tspan,ics) instead.

> If you think the above routine 'myode' returns
> the time derivatives of u in the grid points, it can't
> be true. 

Can you please elaborate? I know there is something wrong with my code and I have spent alot of time trying to find it but no success.
> In the calculation of the array 'out1', you do not use 'u' to get the diffusion terms u1diff and 
> u2diff.
So, I should have something like this instead
u1diff=-delsq(numgrid('S',N+2)).*u1;
u2diff=-delsq(numgrid('S',N+2)).*u2 ?

Many thanks for your assistance.

PS: the code runs out of memory for any values of N above 5.

> Best wishes
> Torsten.