Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Numerical solution to two-dimensional reaction diffusion equation

Subject: Numerical solution to two-dimensional reaction diffusion equation

From: Emma Robertson

Date: 17 Jan, 2011 10:25:04

Message: 1 of 10

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

Subject: Numerical solution to two-dimensional reaction diffusion

From: Torsten Hennig

Date: 18 Jan, 2011 07:26:11

Message: 2 of 10

> 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.

Subject: Numerical solution to two-dimensional reaction diffusion

From: Emma Robertson

Date: 18 Jan, 2011 10:29:05

Message: 3 of 10

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.

Subject: Numerical solution to two-dimensional reaction diffusion equation

From: optiplex

Date: 18 Jan, 2011 10:59:54

Message: 4 of 10

ok thanxs for sharing.

Subject: Numerical solution to two-dimensional reaction diffusion

From: Torsten Hennig

Date: 18 Jan, 2011 11:44:20

Message: 5 of 10

> Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote
> in message
> <994710478.5017.1295335738754.JavaMail.root@gallium.ma
> thforum.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.

Say you have stored the initial values for u1,u2 and u3
in the u-vector in the following order:
u1(1,1),u2(1,1),u3(1,1),u1(2,1),u2(2,1),u3(2,1),...,
u1(nx,1),u2(nx,1),u3(nx,1),
u1(1,2),u2(1,2),u3(1,2),u1(2,2),u2(2,2),u3(2,2),...,
u1(nx,2),u2(nx,2),u3(nx,2),
u1(1,3),u2(1,3),u3(1,3),u1(2,3),u2(2,3),u3(2,3),...,
u1(nx,3),u2(nx,3),u3(nx,3),
..
u1(1,ny),u2(1,ny),u3(1,ny),u1(2,ny),u2(2,ny),
u3(2,ny),...,u1(nx,ny),u2(nx,ny),u3(nx,ny)
(so that the u-vector has a total length of
3*nx*ny elements).

Then in myfun proceed as follows:

#Local variables
For j = 1:ny
For i = 1:nx
u1(i,j) = u(3*(j-1)*nx+3*(i-1)+1);
u2(i,j) = u(3*(j-1)*nx+3*(i-1)+2);
u3(i,j) = u(3*(j-1)*nx+3*(i-1)+3);
end
end


#Calculate time derivatives of solution variables
For i = 2:nx-1
For j = 2:ny-1
du1dt(i,j) = (u1(i+1,j)-2*u1(i,j)+u1(i-1,j))/dx^2;
du1dt(i,j) = du1dt(i,j)+(u1(i,j+1)-2*u1(i,j)+u1(i,j-1))/dy^2;
du1dt(i,j) = du1dt(i,j)+u1(i,j)*u2(i,j)-u1(i,j)+u2(i,j)*u3(i,j);
du2dt(i,j) = (u2(i+1,j)-2*u2(i,j)+u2(i-1,j))/dx^2;
du2dt(i,j) = du2dt(i,j)+(u2(i,j+1)-2*u2(i,j)+u2(i,j-1))/dy^2;
du2dt(i,j) = du2dt(i,j)-u1(i,j)*u2(i,j)+u1(i,j)*u3(i,j)-u2(i,j);
end
end

For i = 1:nx
For j = 1:ny
du3dt(i,j) = u1(i,j)*u2(i,j)-u1(i,j)+u3(i,j)*u1(i,j);
end
end


#Write back time derivatives in global variables
dudt=zeros(3*nx*ny);

For j = 2:ny-1
For i = 2:nx-1
dudt(3*(j-1)*nx+3*(i-1)+1)=du1dt(i,j);
dudt(3*(j-1)*nx+3*(i-1)+2)=du2dt(i,j);
end
end

For j = 1:ny
For i = 1:nx
dudt(3*(j-1)*nx+3*(i-1)+3)=du3dt(i,j);
end
end

Boundary values for u1 and u2 are kept constant at
their initial values.

In myode, return dudt to ODE15s.

Best wishes
Torsten.

Subject: Numerical solution to two-dimensional reaction diffusion

From: Emma Robertson

Date: 18 Jan, 2011 15:30:09

Message: 6 of 10


Torsten, thanks for your reply but why is my code wrong? I thought I had put the diffusion
term in correctly using the 'delsq' command which corresponds to the 2-d laplacian:

u1diff=-delsq(numgrid('S',N+2));

Also, I noticed that there are alot of loops in your code which can be computationally costly if you want to run it for a long time or large N...how can the loops be avoided?

Many thanks for your time.

Subject: Numerical solution to two-dimensional reaction diffusion

From: Torsten Hennig

Date: 18 Jan, 2011 15:57:05

Message: 7 of 10

>
> Torsten, thanks for your reply but why is my code
> wrong? I thought I had put the diffusion
> term in correctly using the 'delsq' command which
> corresponds to the 2-d laplacian:
>
> u1diff=-delsq(numgrid('S',N+2));
>

I don't know what delsq does, but it can't
give d^2u1/dx^2 + d^2u1/dy^2 since the argument
of delsq (namely _numgrid('S',N+2)_) does not
contain any information about u1.
 
> Also, I noticed that there are alot of loops in your
> code which can be computationally costly if you want
> to run it for a long time or large N...how can the
> loops be avoided?
>

I can't think of a way how to avoid the loops.
Maybe MATLAB experts know.

> Many thanks for your time.

Best wishes
Torsten.

Subject: Numerical solution to two-dimensional reaction diffusion

From: Emma Robertson

Date: 18 Jan, 2011 16:21:04

Message: 8 of 10

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <2117839184.7572.1295366255486.JavaMail.root@gallium.mathforum.org>...
> >
> > Torsten, thanks for your reply but why is my code
> > wrong? I thought I had put the diffusion
> > term in correctly using the 'delsq' command which
> > corresponds to the 2-d laplacian:
> >
> > u1diff=-delsq(numgrid('S',N+2));
> >
>
> I don't know what delsq does, but it can't
> give d^2u1/dx^2 + d^2u1/dy^2 since the argument
> of delsq (namely _numgrid('S',N+2)_) does not
> contain any information about u1.
>

'delsq' gives the finite difference discretisation of d/dx^2+d/dy^2 and you are right it doesn't have u1 in it but I think you can make it depend on u1 by multiplying 'delsq' by u1(not sure, ?).

> > Also, I noticed that there are alot of loops in your
> > code which can be computationally costly if you want
> > to run it for a long time or large N...how can the
> > loops be avoided?
> >
>

I am trying to run the code you posted at the moment.

> I can't think of a way how to avoid the loops.
> Maybe MATLAB experts know.
>
> > Many thanks for your time.
>
> Best wishes
> Torsten.

Many thanks.

Subject: Numerical solution to two-dimensional reaction diffusion

From: Emma Robertson

Date: 19 Jan, 2011 10:48:05

Message: 9 of 10

Hi All
Can someone please have a look at my original code and see if you can spot any errors in it, I would greatly appreciate it.
Many thanks
Emma

Subject: Numerical solution to two-dimensional reaction diffusion

From: Philipp Steffen

Date: 14 Mar, 2011 11:45:04

Message: 10 of 10

Hi Emma,

this is of interest for me as well. Have you got the code to work in meanwhile? Unfortunately I am also not an expert in Matlab so far. Could you perhaps share your insight into this with me?
Thanks a lot,

Philipp


"Emma Robertson" wrote in message <ih6fh4$loi$1@fred.mathworks.com>...
> Hi All
> Can someone please have a look at my original code and see if you can spot any errors in it, I would greatly appreciate it.
> Many thanks
> Emma

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us