How Can I Vectorize Heat Equation Diffusion Stepping in Matlab?

[robby991]

Hi, I had developed some Matlab code to solve the heat equation using central difference. I was having issues with my code because I was running out of memory trying to keep the algorithm stable. I need to reconfigure my code so to use single vectors in the time/space stepping, rather than initializing a whole matrix for plotting, so the last step in time and space is a single vector at t=tmax which is used for plotting. A simplified code of how I am currently solving the heat equation is shown below. I am confused how to solve this by vectorizing. Any suggestions would be appreciated.

nx = 30;                      
nt = 100000;             

M = zeros(nt,nx);          

M(:,nx) = M0;                 

for j=1:nt-1  
         
    for i=2:nx-1
      M(j+1,i) = M(j,i) + (dts/dxs^2)*Ds*(M(j,i+1) - 2*M(j,i) + M(j,i-1))-((Vmax*dts*M(j,i))/(Km+M(j,i))); 
     end

       M(j+1,1)=M(j,1)+dts*Ds*2*(M(j,2)-M(j,1))./dxs.^2-((Vmax*dts*M(j,i))/(Km+M(j,i))); %neumann BC
    
end

 

[jvicens]

Thank you for sharing your code and the issues you are facing. It seems like the main issue you are encountering is running out of memory while trying to keep the algorithm stable. I would suggest trying to vectorize your code as much as possible to reduce the amount of memory needed. Here are some suggestions on how to do so:

1. Use single vectors for time and space stepping instead of initializing a whole matrix for plotting. This will significantly reduce the memory needed.

2. Instead of using a for loop for the time stepping, you can use a vectorized approach by using the "dot" operator. For example, you can replace your for loop with the following code:

M(2:end,2:end) = M(1:end-1,2:end) + (dts/dxs^2)*Ds*(M(1:end-1,3:end) - 2*M(1:end-1,2:end) + M(1:end-1,1:end-1))-((Vmax*dts*M(1:end-1,2:end))./(Km+M(1:end-1,2:end)));

This will perform the same calculations as your for loop, but in a vectorized manner which will reduce the memory needed and also improve the performance of your code.

3. You can also try to vectorize the boundary conditions to further reduce the memory needed. For example, instead of using a for loop for the neumann boundary condition, you can use a vectorized approach similar to the one shown above.

I hope these suggestions help you in solving the heat equation using central difference in a more memory-efficient manner. If you have any further questions or need more assistance, please don't hesitate to reach out to me. Best of luck with your project!