This line is wrong
T(2:99, 2:99) = y;
because of the error message given.
Method of Lines 2D
29 次查看（过去 30 天）
显示 更早的评论
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.
回答（1 个）
Torsten
2023-2-3
编辑：Torsten
2023-2-3
25 个评论
Tristan
2023-2-3
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
2023-2-3
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
2023-2-3
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
2023-2-3
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
2023-2-3
编辑：Tristan
2023-2-3
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
2023-2-3
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
2023-2-6
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
2023-2-6
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
2023-2-7
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?
Tristan
2023-2-7
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
2023-2-10
编辑：Tristan
2023-2-10
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
2023-2-10
编辑：Torsten
2023-2-10
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
2023-2-13
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?
Tristan
2023-2-13
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
2023-2-13
编辑：Torsten
2023-2-14
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
2023-2-13
编辑：Tristan
2023-2-13
"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
2023-2-13
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
2023-2-13
编辑：Tristan
2023-2-13
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
2023-2-14
编辑：Torsten
2023-2-14
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
2023-2-14
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?
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Partial Differential Equation Toolbox 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改，无法完成操作。请重新加载页面以查看其更新后的状态。
Translated by 
您也可以从以下列表中选择网站：
如何获得最佳网站性能
选择中国网站（中文或英文）以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)