Method of Lines 2D

So I am trying to solve the 2D heat equation with no source term using method of lines. So basically I want an ODE at each grid point, and solve those using one of the ODE solvers. My code for trying to run the function is below.

T0(1,:) = 35;
T0(:,1) = 0;
T0(100,:) = 0;
T0(:, 100) = 0;

tspan = [0 20];
y0 = T0(2:99, 2:99);
y0 = reshape(y0, [], 1);
[tsol, Tsol] = ode45( @ (t,y) functionheat(t,y), tspan, y0);

The code for the function is here:

function fval = functionheat(t,y)

T(1,:) = 35;
T(:,1) = 0;
T(100,:) = 0;
T(:, 100) = 0;
T(2:99, 2:99) = y;

dx = 0.1/100;
dy = 0.1/100;
dTdt = zeros(100,100);
for i = 2:(nx-1)
    for j = 2:(ny-1)
        Txx = (T(i+1,j) - 2*T(i,j) + T(i-1,j))/(dx^2);
        Tyy = (T(i,j+1)-2*T(i,j)+T(i,j-1))/(dy^2);
        dTdt(i,j) = Txx + Tyy;
                
    end 
end

fval1 = dTdt(2:99, 2:99);
fval = reshape(fval1, [],1);

Originally I just had fval = dTdt(2:99, 2:99);, but it was giving me this error, that I am still getting after adding a new line. The error is below

"Unable to perform assignment because the size of the left side is 98-by-98 and the size of the right side

is 9604-by-1." What exactly needs to change, I am not sure what is still a 98 by 98 matrix because I have reshaped the derivative array which is ultimately what is being solved.

Answer 1

[Torsten]

This line is wrong

T(2:99, 2:99) = y;

because of the error message given.

[Tristan]

Okay, I just changed that line to:

reshape(T(2:99, 2:99), [],1) = y;

And I am now getting this error, "Index in position 1 is invalid. Array indices must be positive integers or logical values." This was basically something else I tried earlier.

Any suggestions on how to reformulate this?

[Torsten]

y = reshape(y,[98,98]);
T = zeros(100);
T(1,:) = 35;
T(:,1) = 0;
T(100,:) = 0;
T(:, 100) = 0;
T(2:99,2:99) = y;

[Tristan]

Okay great, I tried that and I think I am getting somewhere. I did get this error though:

Error using horzcat

Requested 9604x111808 (8.0GB) array exceeds maximum array size preference (8.0GB). This might cause MATLAB

to become unresponsive.

I am assuming this is because I am solving a very high dimensional problem (9604 ODEs), so it sounds like my computer is just having issues handling all that. So essentially is this telling me that I either need to reduce the dimensionality or work on a better computer? Or is this something else going on?

[Torsten]

You should find the place in your code where you request an array of that size. My guess is that you did something wrong. Especially the second dimension of the array cannot be explained in the context of your problem.

If you further want to use ODE15S to solve the 2d heat equation, you should use sparse matrices. Especially supplying the constant Jacobian for the solver in sparse matrix format is mandatory. It will speed up your calculations enormously and save a great amount of RAM.

[Tristan]

Agreed. The 9604 dimension makes sense as that is obviously the grid locations 98x98 since boundaries are fixed. Could the second dimension be the number of time steps or something? Not sure why it would be so high, I am just thinking.

So you are recommending switching from ode45 to ode15s? Because as of now I am using ode45.

[Torsten]

Could the second dimension be the number of time steps or something? Not sure why it would be so high, I am just thinking.

I don't know how many solutions you save within your tspan. If you chose 111808 output times, you are correct.

So you are recommending switching from ode45 to ode15s? Because as of now I am using ode45.

Systems of ODEs resulting from the discretization of second-order derivatives (like those of the heat equation) are stiff. If you follow the advices concerning sparse computations and Jacobian supply, I'd recommend switching to ode15s because it will speed up your computations.

[Tristan]

I did not specify an amount of solution within tspan, so I guess it was left as default. Perhaps I can try to reduce it and see.

Could you give any guidance on how exactly I determine the Jacobian that I need to be supplying? I didn't realize this set of ODEs are stiff.

[Torsten]

Most probably, you only define start and end time for the integrator. This means that the integrator will save every basic integration step and return it to your calling program. Since the system is stiff and you use ode45, 111808 small time steps are needed. This is the reason for your large output matrix.

So change the solver and supply the Jacobian.

The equations you solve read

dT(i,j)/dt = 
(T(i-1,j)-2*T(i,j)+T(i+1,j))/dx^2 + (T(i,j-1)-2*T(i,j)+T(i,j+1))/dy^2 =: 
f(T(i-1,j),T(i+1,j),T(i,j-1),T(i,j+1),T(i,j))

So you see that the derivatives of the (i,j)-th equation are all zero except for the derivatives with respect to

T(i-1,j),T(i+1,j),T(i,j-1),T(i,j+1),T(i,j)

which are

1/dx^2, 1/dx^2, 1/dy^2, 1/dy^2, -2*(1/dx^2+1/dy^2)

respectively.

You already chose an order of the variables T(i,j) by building the column vector

reshape(T, [], 1);

So I think you should be able to build the Jacobian now.

[Tristan]

I see so in this case, the Jacobian would be a 5x5 matrix correct?

[Torsten]

The Jacobian matrix is a (98x98) x (98x98) = 9604 x 9604 matrix for your example.

But each row of this matrix only has maximum 5 elements different from 0.

[Tristan]

Okay, I think I somewhat understand. I don't have an immediate idea, but I can try something. I guess I was thinking it would be 98x98 rather than (98x98)x(98x98), so I am a bit confused why there are that many elements?

[Torsten]

The Jacobian of a system of n equations in n unknowns has size nxn.

You have a system of 98^2 equations in 98^2 unknowns, thus ...

[Tristan]

Okay right, I am following now.

So how can we say where the position of the 5 nonzero elements are in each row?

Would the first row of the Jacobian basically be the derivatives of dTdt(1,1), then the next row derivatives of dTdt(1,2), etc?

So the nonzero elements would always be at the location on the row that corresponds to the i,j location, then the i-1,j i+1,j i,j+1 and i,j-1 locations on the grid?

I am just trying to brainstorm how I can formulate this in a loop.

[Tristan]

One other question, if I don't necessarily care about computational time, I can ignore specifying the Jacobian correct?

In other words, using the stiff solver has still yielded solution of the set of ODEs, but I am not sure if it looks correct. The Jacobian being specified would only have impact on time and not accuracy, correct?

[Torsten]

Specifying the Jacobian would have impact on speed and - more important - RAM required for larger problems, but not on accuracy.

You can imagine that it makes a difference if - without the Jacobian option - MATLAB assumes a full (98*98)x(98*98) matrix or - with the Jacobian option- works with a sparse matrix of size (98*98)x(98*98) with only 98*98*5 elements different from 0.

[Tristan]

Yes, thank you that makes sense. For now, however, I am not too concerned with RAM and speed. I have seen an implementation in Python of what I am trying to do, but done a bit differently. I want to compare their results to mine. They use some kind of animation plotting.

Currently my code is:

Tuse = zeros(100, 100);
Tuse(1,:) = 10;
T0 = zeros(100,100);
dx = 0.1/100;
dy = 0.1/100;


tspan = [0 5];
y0 = T0(2:99,2:99);
y0 = reshape(y0, [], 1);
[tsol, Tsol] = ode15s( @ (t,y) functionheat(t,y), tspan, y0);
[m,n] = size(tsol);
for i = 1:m
    Tgrid = reshape(Tsol, 98,98,[]);
    Tuse(2:99, 2:99) = Tgrid(:,:,m);
    surf(Tuse(:,:))
end 

Right now, it seems that the surface plot returns just the first time instance. I have seen MATLAB has a heat map function but it requires input as a table. Do you have any reccommendation for visualizing the solution as time goes on, potentially as an animation?

[Torsten]

Please supply the complete code you are using at the moment (including "functionheat").

This might be what you want:

https://de.mathworks.com/help/matlab/creating_plots/animating-a-surface.html

[Tristan]

function fval = functionheat(t,y)

y = reshape(y,[98,98]);
T = zeros(98);
T(1,:) = 10;
T(:,1) = 0;
T(100,:) = 0;
T(:, 100) = 0;
T(2:99,2:99) = y;
nx = 100;
ny = 100;
dx = 0.1/100;
dy = 0.1/100;
dTdt = zeros(100,100);
alpha = 2;
for i = 2:(nx-1)
    for j = 2:(ny-1)
        Txx = (T(i+1,j) - 2*T(i,j) + T(i-1,j))/(dx^2);
        Tyy = (T(i,j+1)-2*T(i,j)+T(i,j-1))/(dy^2);
        dTdt(i,j) = alpha*(Txx + Tyy);
                
    end 
end
dTdt = dTdt(2:99,2:99);
dTdtuse = reshape(dTdt, [], 1);
fval = dTdtuse;
end

The above is the function, and the code from my previous reply is the rest. I think that could work, however I would prefer if there was a way to plot more of a heat map where the values are shown in a 2D plane and vary in color across that 2D space without rising out of the plane in 3D space.

[Torsten]

I didn't check if your reshaping from matrix to vector and vice versa is correct. You should test it independent of this code with small matrices and column vectors.

The result looks strange.

Tuse = zeros(11);
Tuse(1,:) = 10;
T0 = zeros(11);
nx = 11;
ny = 11;
dx = 0.1/(nx-1);
dy = 0.1/(ny-1);
tspan = [0,0.1,100];
y0 = T0(2:10,2:10);
y0 = reshape(y0, [], 1);
[tsol, Tsol] = ode15s( @ (t,y) functionheat(t,y), tspan, y0);
x = 0:dx:0.1;
y = 0:dy:0.1;
[X,Y] = meshgrid(x,y);
[m,n] = size(tsol);
for i=1:m
  Tgrid = reshape(Tsol,9,9,[]);
  Tuse(2:10,2:10) = Tgrid(:,:,m);
  contourf(X,Y,Tuse);
  pause(2)
end
colorbar

function fval = functionheat(t,y)

y = reshape(y,[9,9]);
T = zeros(11);
T(1,:) = 10;
T(:,1) = 0;
T(11,:) = 0;
T(:, 11) = 0;
T(2:10,2:10) = y;
nx = 11;
ny = 11;
dx = 0.1/(nx-1);
dy = 0.1/(ny-1);
dTdt = zeros(11);
alpha = 2;
for i = 2:nx-1
    for j = 2:ny-1
        Txx = (T(i+1,j) - 2*T(i,j) + T(i-1,j))/(dx^2);
        Tyy = (T(i,j+1)-2*T(i,j)+T(i,j-1))/(dy^2);
        dTdt(i,j) = alpha*(Txx + Tyy);
                
    end 
end
dTdt = dTdt(2:10,2:10);
dTdtuse = reshape(dTdt, [], 1);
fval = dTdtuse;
end

[Tristan]

"I didn't check if your reshaping from matrix to vector and vice versa is correct. You should test it independent of this code with small matrices and column vectors."

When you say it are you referring to the testing you showed talking about the heat map or are you saying I need to test the reshaping? I am pretty sure you were saying the former, just want to make sure.

That result is strange. Any idea why that would be the case? I think the bottom part looks okay and makes physical sense, but I am not sure what those random diamond things are and why there are appearing. I don't think this would have cause that but shouldn't dx and dy be 0.1/11 since it seems you are using an 11x11 grid? Perhaps the number is too low and somehow that causes an issue?

[Torsten]

When you say it are you referring to the testing you showed talking about the heat map or are you saying I need to test the reshaping? I am pretty sure you were saying the former, just want to make sure.

I mean all parts of your code where you reshape from matrix to column vector or from column vector to matrix. This must be done consistently.

I don't think this would have cause that but shouldn't dx and dy be 0.1/11 since it seems you are using an 11x11 grid? Perhaps the number is too low and somehow that causes an issue?

No. dx = dy = 0.1/10 for 11 points. The number of points is always by 1 greater than the number of subintervals. Count it.

And no. You won't see such effects for heat conduction problems although the resolution is coarse. But you can easily test it - your former case had (should have) a grid of 101 points.

[Tristan]

Ah, I see what you meant. Yes, I think the reshaping is fine, but I can definitely go back and verify again.

I see what you mean about the number of points being greater than subintervals. I meant to have a 100x100 grid and all my code was written with that in mind, so I made changes there.

I will try to plot doing something similar with my setup and get back, MATLAB is acting weird for me right now and isn't running things right. Do you have any hypothesis for the result you got with a more coarse resolution? I'm not sure how those spots can even become nonzero if the inital condition at those points were 0.

[Torsten]

If you like, you can remove the loops for reshaping. But be careful !

nx = 11;
ny = 11;
T0 = zeros(nx,ny);
T0(1,:) = 10;
dx = 0.1/(nx-1);
dy = 0.1/(ny-1);
tspan = [0,500,1000];
for iy = 1:ny
  for ix = 1:nx
    y0((iy-1)*nx + ix) = T0(ix,iy);
  end
end
[tsol, Tsol] = ode15s( @ (t,y) functionheat(t,y,nx,ny,dx,dy), tspan, y0);

x = 0:dx:0.1;
y = 0:dy:0.1;
[X,Y] = ndgrid(x,y);

Tgrid = zeros(nx,ny,numel(tsol));
for it = 1:numel(tsol)
  for iy = 1:ny
    for ix = 1:nx
      Tgrid(ix,iy,it) = Tsol(it,(iy-1)*nx + ix);
    end
  end
end

for it = 1:numel(tsol)
  for ix = 1:nx
    for iy = 1:ny
      Tuse(ix,iy) = Tgrid(ix,iy,it);    
    end
  end
  contourf(X,Y,Tuse)
  drawnow
  %pause(10)
end
colorbar

function fval = functionheat(t,y,nx,ny,dx,dy)

  for iy = 1:ny
    for ix = 1:nx
      T(ix,iy) = y((iy-1)*nx + ix);
    end
  end

  dTdt = zeros(nx,ny);
  alpha = 2;
  for ix = 2:nx-1
    for iy = 2:ny-1
        Txx = (T(ix+1,iy) - 2*T(ix,iy) + T(ix-1,iy))/(dx^2);
        Tyy = (T(ix,iy+1) - 2*T(ix,iy) + T(ix,iy-1))/(dy^2);
        dTdt(ix,iy) = alpha*(Txx + Tyy);
    end
  end

  fval = zeros(nx*ny,1);
  for iy = 1:ny
    for ix = 1:nx
      fval((iy-1)*nx + ix) = dTdt(ix,iy);
    end
  end

end

[Tristan]

Oh okay, wow! This looks great.

So my previous reshaping was wrong? Or what change do you think was able to fix those pods away from the bounday?

[Torsten]

I like programming more than debugging ... So I didn't check. But as you can see, same coarse resolution, no spikes.