Adding Non-Linear source term to 2d Implicit MATLAB code - Computational Science Stack Exchange most recent 30 from scicomp.stackexchange.com 2019-11-17T14:54:08Z https://scicomp.stackexchange.com/feeds/question/28281 https://creativecommons.org/licenses/by-sa/4.0/rdf https://scicomp.stackexchange.com/q/28281 -2 Adding Non-Linear source term to 2d Implicit MATLAB code DesperateStudent https://scicomp.stackexchange.com/users/26026 2017-11-20T10:40:08Z 2017-11-30T10:35:13Z <p>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: </p> <p>$$\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)$$</p> <p>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.</p> <p>thank you so much.</p> <pre class="lang-matlab prettyprint-override"><code>%% %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&lt;=y(j))&amp;&amp;(y(j)&lt;=1.5)&amp;&amp;(0.5&lt;=x(i))&amp;&amp;(x(i)&lt;=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 shading interp 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 </code></pre> https://scicomp.stackexchange.com/questions/28281/adding-non-linear-source-term-to-2d-implicit-matlab-code/28304#28304 1 Answer by cbcoutinho for Adding Non-Linear source term to 2d Implicit MATLAB code cbcoutinho https://scicomp.stackexchange.com/users/19791 2017-11-22T16:24:28Z 2017-11-22T18:49:02Z <p>You can solve your problem in two ways: <em>explicitly</em> or <em>implicitly</em>.</p> <p>If you want to solve your system <em>explicitly</em>, 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 <em>small</em> 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. <em>ode45</em>, <em>ode15</em>, etc.). This has the added benefit of using <em>adaptive timesteps</em> because those Matlab routines are able to calculate the size of a timestep based on error measurements.</p> <p>If you want to solve your equation <em>implicitly</em>, then you move into the world of <em>non-linear</em> PDEs, which require a non-linear solver within every time step. This is more expensive than explicit solvers <em>when using the same sized timestep</em>; 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 <em>implicit</em> 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.</p> <hr> <p>A rather simple non-linear implicit solver to implement is called the <em>Newton-Raphson method</em>, and may be useful for your problem.</p> <p>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 <em>linear</em> terms <em>implicitly</em> and <em>non-linear</em> terms <em>explicitly</em>. 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.</p>