f



Numerical solution to two-dimensional reaction diffusion equation

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
0
Emma
1/17/2011 10:25:04 AM
comp.soft-sys.matlab 211264 articles. 26 followers. lunamoonmoon (257) is leader. Post Follow

3 Replies
2380 Views

Similar Articles

[PageSpeed] 22

> 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.
0
Torsten
1/18/2011 2:26:11 AM
ok thanxs for sharing.
0
optiplex
1/18/2011 5:59:54 AM
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.
0
Emma
1/18/2011 10:29:05 AM
Reply: