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 |

1/17/2011 10:25:04 AM

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

1/18/2011 2:26:11 AM

ok thanxs for sharing.

0 |

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 |

1/18/2011 10:29:05 AM