# Adding Non-Linear source term to 2d Implicit MATLAB code

I'm running out of time for this code so any help would be greatly appreciated. I am currently coding the 2D heat/diffusion equation in matlab but i'm having trouble adding in the source term. my equation looks like this:

$$\frac{\partial u}{\partial t} = \alpha_x \frac{\partial^2 u}{\partial x^2} + \alpha_y \frac{\partial^2 u}{\partial y^2} + Au \left( 1 - \frac{u}{K} \right)$$

I'm using Neumann bc's with the flux=0 everywhere. u is the population number, A, the Intrinsic growth rate and K the carrying capacity. The code I'm using is a code I found here that I've edited to use the implicit/BTCS method.

thank you so much.

%%
%Specifying parameters
nx=40;                           %Number of steps in space(x)
ny=50;                           %Number of steps in space(y)
nt=50;                           %Number of time steps
dt=0.01;                         %Width of each time step
dx=2/(nx-1);                     %Width of space step(x)
dy=2/(ny-1);                     %Width of space step(y)
x=0:dx:2;                        %Range of x(0,2) and specifying the grid points
y=0:dy:2;                        %Range of y(0,2) and specifying the grid points
u=zeros(nx,ny);                  %Preallocating u
un=zeros(nx,ny);                 %Preallocating un
visx = 2;                        %Diffusion coefficient/viscocity
visy = 2;                        %Diffusion coefficient/viscocity
% RepRateA = 0.3;                %Intrinsic Growth Rate of Population
% CarCapK = 1.7;                 %Carrying Capacity of Environment

UnW=0;                           %x=0 Neumann B.C (du/dn=UnW)
UnE=0;                           %x=L Neumann B.C (du/dn=UnE)
UnS=0;                           %y=0 Neumann B.C (du/dn=UnS)
UnN=0;                           %y=L Neumann B.C (du/dn=UnN)

%%
%Initial Conditions
for i=1:nx
for j=1:ny
if ((0.5<=y(j))&&(y(j)<=1.5)&&(0.5<=x(i))&&(x(i)<=1.5))
u(i,j)=2;
else
u(i,j)=0;
end
end
end

%%
%B.C vector
bc=zeros(nx-2,ny-2);

bc(1,:)=-UnW/dx;       %Neumann B.Cs
bc(nx-2,:)=UnE/dx;     %Neumann B.Cs
bc(:,1)=-UnS/dy;       %Neumann B.Cs
bc(:,nx-2)=UnN/dy;     %Neumann B.Cs

%B.Cs at the corners:
bc=visx*dt*bc+visy*dt*bc;

%Calculating the coefficient matrix for the implicit scheme
Ex=sparse(2:nx-2,1:nx-3,1,nx-2,nx-2);
% Ax=Ex+Ex'-2*speye(nx-2);
Ax(1,1)=-1;
Ax(nx-2,nx-2)=-1;                 %Neumann B.Cs
Ey=sparse(2:ny-2,1:ny-3,1,ny-2,ny-2);
% Ay=Ey+Ey'-2*speye(ny-2);
Ay(1,1)=-1;
Ay(ny-2,ny-2)=-1;                 %Neumann B.Cs
A=kron(Ay/dy^2,speye(nx-2))+kron(speye(ny-2),Ax/dx^2);
D=speye((nx-2)*(ny-2))-visx*dt*A-visy*dt*A;
%%
%Calculating the field variable for each time step
i=2:nx-1;
j=2:ny-1;

for it=0:nt
un=u;
h=surf(x,y,u','EdgeColor','none');       %plotting the field variable
axis ([0 2 0 2 0 2])
title({['2-D Diffusion with alphax = ',num2str(visx)];['2-D Diffusion with alphay = ',num2str(visy)];['time (\itt) = ',num2str(it*dt)]})
xlabel('Spatial co-ordinate (x) \rightarrow')
ylabel('{\leftarrow} Spatial co-ordinate (y)')
zlabel('Transport property profile (u) \rightarrow')
drawnow;
refreshdata(h)

%Implicit method:
U=un;
U(1,:)=[];
U(end,:)=[];
U(:,1)=[];
U(:,end)=[];
U=reshape(U+bc,[],1);
U=D\U;
U=reshape(U,nx-2,ny-2);
u(2:nx-1,2:ny-1)=U;
%Neumann:
u(1,:)=RepRateA*u(2,:)*(1-u(2,:).'/CarCapK)-UnW*dx;
u(nx,:)=RepRateA*u(nx-1,:)*(1-u(nx-1,:).'/CarCapK)+UnE*dx;
u(:,1)=RepRateA*u(:,2).'*(1-u(:,2)/CarCapK)-UnS*dy;
u(:,ny)=RepRateA*u(:,ny-1).'*(1-u(:,ny-1)/CarCapK)+UnN*dy;
end

• Not sure why it looks like this, I apologise. first time poster. hope the layout doesn't deter anyone – DesperateStudent Nov 20 '17 at 10:41
• What's the issue you're having? I don't understand what you mean by 'having trouble adding in the source code' – cbcoutinho Nov 21 '17 at 20:32
• @cbcoutinho I'm not sure where to code the term Au(1-u/K) nor how exactly i should go about doing it. also, if i were to change the diffusion coefficients (visx and visy) to variables or functions, how would i go about doing that? – DesperateStudent Nov 22 '17 at 8:15

You can solve your problem in two ways: explicitly or implicitly.

If you want to solve your system explicitly, then the solution is quite simple: the only term you would be solving for is the $u_{i+1}$ term in the time derivative. All the other instances of $u$ are known quantities ($u_{i}$), which means they can be placed on the RHS - because they are knowns! To ensure numerical stability, the size of the time step should be small enough that the difference between $u_i$ and $u_{i+1}$ is small. The definition of small is not trivial in some cases. Most textbooks on PDEs will give you more information on what stable means w.r.t explicit solvers of PDEs. If you really want to get fancy you could set up your PDE into a system of ODEs in time and use one of the Matlab supplied ODE solvers (e.g. ode45, ode15, etc.). This has the added benefit of using adaptive timesteps because those Matlab routines are able to calculate the size of a timestep based on error measurements.

If you want to solve your equation implicitly, then you move into the world of non-linear PDEs, which require a non-linear solver within every time step. This is more expensive than explicit solvers when using the same sized timestep; however, implicit solvers are commonly more stable than explicit ones and can sometimes warrant the use of larger timesteps to speed up the solver. I noticed you wrote implicit in your code above: if that's what you decide to do then you would need to write a non-linear solver for your problem.

A rather simple non-linear implicit solver to implement is called the Newton-Raphson method, and may be useful for your problem.

If you don't want to implement a non-linear implicit solver, you may instead decide to implement a hybrid approach where you treat all linear terms implicitly and non-linear terms explicitly. What is best for solving your problem can only be determined by looking more closely at literature on the subject or experimentation. Since you're already working with the BTCS method, I would recommend extending it by including your source term explicitly into the RHS. This would result in a 'lag' of your source term with respect to time, but may be reasonable for your problem.

• Thank you so much. i would prefer to code it directly into the RHS but i've fiddled around with that on both U sections and the u section of the RHS, and both times i've received errors. so i'm not certain where exactly this should go. When added to line 72 as un=RepRateA*u*(1-u.'/CarCapK) the error comes at line 89 {U=reshape(U+bc,[],1)} and says that "Matrix dimensions must agree". – DesperateStudent Nov 23 '17 at 8:09
• When i code for U separately, i get the same error in line 89 as a mentioned above. Error in ==> U(:,1)=[RepRateAU(:,2)*(1-U(:,2).'/CarCapK)]; U(:,end)=[RepRateAU(:,end-1)*(1-U(:,end-1).'/CarCapK)] etc but when I add the code in at u, (The last bit of code in the program) u(:,1) = RepRateAu(:,2)*(1-u(:,2).'/CarCapK)-UnSdy u(:,ny)=RepRateAu(:,ny-1)*(1-u(:,ny-1).'/CarCapK)+UnNdy etc instead of an error i just have the results go haywire, especially at the boundaries – DesperateStudent Nov 23 '17 at 8:50
• I don't the code you copied above is complete - there's no end statement associated with the for loop. I also don't see the line that you say is giving you problems un=RepRateA*u.... Make the code example above match the one you are using and it will make it easier to help you – cbcoutinho Nov 23 '17 at 12:14
• As i previously stated, i wasn't sure where to include the term Au(1-u/K) so ln the example code i gave, i left it out. – DesperateStudent Nov 26 '17 at 14:30
• after the for it=0:nt the next lione of code says un=u. when i change that to un=RepRateA*u*(1-u.'/CarCapK) then the error "Matrix dimensions must agree" pops up. – DesperateStudent Nov 26 '17 at 14:37