## Reaction-Diffusion Equations and Matlab

## Reaction-Diffusion Equations and Matlab

(OP)

I want to solve the reaction-diffusion problem, in 2D, with Matlab:

Sum_{j=1}^{4} D_{1j} * (drond^{2}z_{j} / drond x^{2} + drond^{2}z_{j} / drond y^{2}) + R_{1} = 0

Sum_{j=1}^{4} D_{2j} * (drond^{2}z_{j} / drond x^{2} + drond^{2}z_{j} / drond y^{2}) + R_{2} = 0

Sum_{j=1}^{4} D_{3j} * (drond^{2}z_{j} / drond x^{2} + drond^{2}z_{j} / drond y^{2}) + R_{3} = 0

Sum_{j=1}^{4} D_{4j} * (drond^{2}z_{j} / drond x^{2} + drond^{2}z_{j} / drond y^{2}) + R_{4} = 0

Boundary conditions: x = 0,1 z_{1} = 1 ; x = 0,1 z_{2,3,4} = 0

So, I have a system with four partial diffrential equations. and I want to plot the solution for every equations.

What it's wrong in my code ? The plots of solutions must be a surface (almost like a bell), but is not.

Sum_{j=1}^{4} D_{1j} * (drond^{2}z_{j} / drond x^{2} + drond^{2}z_{j} / drond y^{2}) + R_{1} = 0

Sum_{j=1}^{4} D_{2j} * (drond^{2}z_{j} / drond x^{2} + drond^{2}z_{j} / drond y^{2}) + R_{2} = 0

Sum_{j=1}^{4} D_{3j} * (drond^{2}z_{j} / drond x^{2} + drond^{2}z_{j} / drond y^{2}) + R_{3} = 0

Sum_{j=1}^{4} D_{4j} * (drond^{2}z_{j} / drond x^{2} + drond^{2}z_{j} / drond y^{2}) + R_{4} = 0

Boundary conditions: x = 0,1 z_{1} = 1 ; x = 0,1 z_{2,3,4} = 0

So, I have a system with four partial diffrential equations. and I want to plot the solution for every equations.

What it's wrong in my code ? The plots of solutions must be a surface (almost like a bell), but is not.

#### CODE --> Matlab

clc; clear; close all; a=0;b=1;n=5; h = (b-a)/(n+1); % step I=eye((n-1)*(n-1)); A11=zeros((n-1)*(n-1));A12=zeros((n-1)*(n-1));A13=zeros((n-1)*(n-1)); A14=zeros((n-1)*(n-1)); A21=zeros((n-1)*(n-1));A22=zeros((n-1)*(n-1));A23=zeros((n-1)*(n-1));A24=zeros((n-1)*(n-1)); A31=zeros((n-1)*(n-1));A32=zeros((n-1)*(n-1));A33=zeros((n-1)*(n-1));A34=zeros((n-1)*(n-1)); A41=zeros((n-1)*(n-1));A42=zeros((n-1)*(n-1));A43=zeros((n-1)*(n-1));A44=zeros((n-1)*(n-1)); b1=zeros((n-1)*(n-1),1);b2=zeros((n-1)*(n-1),1);b3=zeros((n-1)*(n-1),1);b4=zeros((n-1)*(n-1),1); k1=1;k2=0.4;k3=0.5;k4=0.2;k5=0.1; % diffusion matrix D11=1; D12=-0.43; D13=-0.19; D14=-0.265; D21=-0.14; D22=1.475; D23=-0.025; D24=0.235; D31=0.105; D32=0.07; D33=1.475; D34=-0.19; D41=-0.085; D42=-0.245; D43=0.27; D44=1.17; [x,y] = meshgrid(a:h:b); % uniform mesh % Matrix A K1D = spdiags(ones(n,1)*[-1 2 -1],-1:1,n-1,n-1); I1D = speye(size(K1D)); Delta = kron(K1D,I1D)+kron(I1D,K1D); % A1i A11=D11*Delta-(k1+k3)*I*h^2; A12=D12*Delta; A13=D13*Delta; A14=D14*Delta-k5*I*h^2; %A2i A21=D21*Delta-k1*I*h^2; A22=D22*Delta+k2*I*h^2; A23=D23*Delta-k4*I*h^2; A24=D24*Delta; %A3i A31=D31*Delta; A32=D32*Delta-k2*I*h^2; A33=D33*Delta+k4*I*h^2; A34=D34*Delta; %A4i A41=D41*Delta-k3*I*h^2; A42=D42*Delta; A43=D43*Delta; A44=D44*Delta+k5*I*h^2; % b b1(1)=(-1/h^2)*D11; b2(1)=(-1/h^2)*D21; b3(1)=(-1/h^2)*D31; b4(1)=(-1/h^2)*D41; b1(n-1)=(-1/h^2)*D11; b2(n-1)=(-1/h^2)*D21; b3(n-1)=(-1/h^2)*D31; b4(n-1)=(-1/h^2)*D41; % aprox solution A=[A11 A12 A13 A14;A21 A22 A23 A24;A31 A32 A33 A34;A41 A42 A43 A44]; bn=[b1;b2;b3;b4]; ub=reshape(bn,(n-1)*(n-1)*4,1); v=inv(A)*ub; % split the solution vector pas=size(A11); v1=v(1:n+2); v2=v(pas+1:2*pas); v3=v(2*pas+1:3*pas); v4=v(3*pas+1:4*pas); vv1=zeros(n+2,n+2); vv2=zeros(n+2,n+2); vv3=zeros(n+2,n+2); vv4=zeros(n+2,n+2); % boundary conditions vv1=[ones(n+2,1) v1 ones(n+2,1) ones(n+2,1) ones(n+2,1) ones(n+2,1) ones(n+2,1)]; vv1=[zeros(n+2,1) v1 zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1)]; vv1=[zeros(n+2,1) v1 zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1)]; vv1=[zeros(n+2,1) v1 zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1)]; % plot figure set(gcf,'DefaultAxesFontSize',8,'PaperPosition', [0 0 3.5 3.5]) subplot(2,2,1); surf(x,y,vv1); title('Numerical solution for Ecuation1'); subplot(2,2,2); surf(x,y,vv2); title('Numerical solution for Ecuation2'); subplot(2,2,3); surf(x,y,vv3) title('Numerical solution for Ecuation3'); subplot(2,2,4); surf(x,y,vv4); title('Numerical solution for Ecuation4');

## RE: Reaction-Diffusion Equations and Matlab

vv1=[ones(n+2,1) v1 ones(n+2,1) ones(n+2,1) ones(n+2,1) ones(n+2,1) ones(n+2,1)];

vv1=[zeros(n+2,1) v1 zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1)];

vv1=[zeros(n+2,1) v1 zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1)];

vv1=[zeros(n+2,1) v1 zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1)];

makes no sense to me, and you don't seem to do anything much with vv2 etc

Cheers

Greg Locock

New here? Try reading these, they might help FAQ731-376: Eng-Tips.com Forum Policies http://eng-tips.com/market.cfm?

## RE: Reaction-Diffusion Equations and Matlab

Thanks for your reply.

I think my issues start when I try to split the big vector v, into 4 small vectors (one for each equation:v1,v2,v3,v4).

After that, the boundary conditions(vv1,vv2,vv3,vv4) is not correctly implemented (it's need to be created dynamically).

In conclusion, there is a size issue between:

- [x,y] = meshgrid(a:h:b);

- split the vector v

- surf(x,y,vv1);

But, I don't have any idea to resolv the issue.

Cheers,

George Sarbu

## RE: Reaction-Diffusion Equations and Matlab

v2=v(pas+1:2*pas);

pas is a 1x2 matrix, what do you think that line is supposed to accomplish? I suspect that you really want

v2=v(pas(1)+1:2*pas(1));

doubtless there are similar problems with array sizes later on

Cheers

Greg Locock

New here? Try reading these, they might help FAQ731-376: Eng-Tips.com Forum Policies http://eng-tips.com/market.cfm?

## RE: Reaction-Diffusion Equations and Matlab

## CODE --> Matlab

Now, if I take h very small: h = (b-a)/(n+10), the surf function run well. But, of course the plot is not like a bell :)

## RE: Reaction-Diffusion Equations and Matlab

Cheers

Greg Locock

New here? Try reading these, they might help FAQ731-376: Eng-Tips.com Forum Policies http://eng-tips.com/market.cfm?