Plotting the heat equation using the explicit method

1 次查看（过去 30 天）
David
David 2015-3-6
评论： David 2015-3-7
Hi, I am supposed to use the explicit method to plot an approximation of the heat equation in Matlab. The heat equation is as follows:
du/dx=d^2u/d^2x (u_t=u_xx).
Initial conditions: u(x,0)=1 if x>0, while if x is equal or greater than zero u(x,0)=0.
The explicit method is the following: u(j,m+1) =r*u(j-1,m) + (1-2*r)*u(j,m)+r*u(j+1,m), where u is the solution, j is which x-value while m is which time value it is.
The matlab-code is the following:
L = 2.; % Length of the bar
T =0.1; % Time
maxm = 2000; % Time steps
dt = T/maxm;
n = 70; % Distance steps
dx = L/n;
r = dt/(dx^2); % Stability parameter, r less or equal to 1/2
for j=1:n+1
x(j)=-10+(j-1)*dx; %Because j actually starts at zero
if x> 0
u(j,1)=1;
else
u(j,1)=0;
end
end
for m=1:maxm % Time Loop
u(1,m)=0;
u(n+1,m)=1;
for j=2:n; % Space Loop
u(j,m+1) =r*u(j-1,m) + (1-2*r)*u(j,m)+r*u(j+1,m);
end
end
figure(2)
plot(x,u(:,1))
I am not sure what is wrong. When I plot this, it equals zero for all x-values. It was supposed to be equal to 1 when x>0, while zero when x was less than zero. Can anyone help me to figure out was is wrong with this Matlab-code?
David
  10 个评论
显示 8更早的评论
Torsten
Torsten 2015-3-6
For each time step, the value of u at the left and at the right boundary has to be fixed - thus u(1,m)=0 and u(n+1,m)=1 is not only possible, it is necessary.
And what do you mean by "(Because it is not supposed to be zero at x=0 when m not equal to zero)." u is not zero at x=0, but at x=-10.
You can easily check whether your implementation is correct by comparing with the analytical solution. It is given by
u(x,t)=0.5*(2-erfc(x/(2*sqrt(t))))
where erfc is the complementary error function (available in MATLAB).
Best wishes
Torsten.
David
David 2015-3-7
The approximation was a little bit steeper than the analytical solution, so that was good. Thank you very much for all your help. I appreciate that very much.
David

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Titus Edelhofer
Titus Edelhofer 2015-3-6
Hi,
if the bar is in the negative x as well, I would rewrite the loop entirely to be more like MATLAB style, something like
L = 2;
% divide into n pieces:
x = linspace(-L/2, L/2, n+1)';
% declare u
u = zeros(n+1, 1);
% set values of u
u(x>=0) = 1;
Titus
  3 个评论
显示 1更早的评论
Titus Edelhofer
Titus Edelhofer 2015-3-6
Yeah, should have been
u = zeros(n+1, maxm);
so that you preallocate u for all time steps. And then it would be
u(x>=0, 1) = 1;
for the first column (initial time step).
Titus
David
David 2015-3-7
Thank you for your help. It is much appreciated.
David

请先登录，再进行评论。

更多回答（1 个）

Titus Edelhofer
Titus Edelhofer 2015-3-6
Hi,
when you set the value for u, you need to replace
if x>0
by
if x(j)>0
Then you will see, that your u looks like you want instead of identical zero.
Titus
  3 个评论
显示 1更早的评论
Titus Edelhofer
Titus Edelhofer 2015-3-6
Maybe you should describe first in more detail what the problem is. Having a single point with a value different from the others doesn't sound good to me. Suppose you would refine the grid, then the proportion of the "bar" having u=0 instead of u=1 would go further down ...?
David
David 2015-3-6
Okay. My problem is that I am supposed use the explicit method to find an approximation for the heat equation with the following initial value: u(x,0)=1 when x>0, and u(x,0)=1 when x is less or equal to zero. That is all the information about the heat equation I am given. I am not given any other boundary conditions. The problem is then to find the correct matlab-code for this task. I do not think the matlab-code I have written above is 100% correct. My question is what I should change in the code for it to be correct according to the heat equation and its initial conditions.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 General Physics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by