Plotting the heat equation using the explicit method

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

[Torsten]

What do you expect if you set u=0 at both boundaries ?

Best wishes

Torsten.

[David]

That is a good point. Since I am not given any specific boundaries, I have now deleted them in the matlab-code. Then I got the following error-message:Undefined function or variable 'm'.

Error in heat_equation_explicit (line 30) u(j,m+1) =r*u(j-1,m) + (1-2*r)*u(j,m)+r*u(j+1,m);

Do you know how I can "fix" this error in the code?

[Torsten]

You still use u=0 at both boundaries. You will see this if you look at u(:,1).

For the error message, I don't have an explanation. Maybe the wrong semicolon behind "for j=2:n".

Best wishes

Torsten.

[David]

Okay, thank you for your reply. There is now no error messages. When I write u(:,1) I get 0 as the first value, then 1 at all the others. Do you see if there anything I should write differently in order for the code to represent the heat equation with initial condition:

 u(x,0)=1 when x>0, and

u(x,0)=1 when x is less or equal to zero? It has to be written by using the explicit method.

David

[Torsten]

Why do you use a certain length for a bar ?

According to the problem formulation, you should solve the heat equation over (-oo;oo). Because the solution will be antisymmetric around x=0, you should choose the solution interval to be [-10:10], e.g. with boundary value 0 at x=-10 and 1 at x=10.

Best wishes

Torsten.

[David]

Yes, that is absolutely correct, but how do you choose negative x-values as well? When I try to do that by setting j to go from negative n to n, I get an error message which says: "index must be a positive integer or logical."

About the boundary values. Are you sure I am "allowed" to choose boundary values 0 at x=-10 and 1 at x=10? Since I am only given two initial conditions from the text, but no boundary conditions when t not equal to zero.

David

[Torsten]

>Yes, that is absolutely correct, but how do you choose negative x-values as >well? When I try to do that by setting j to go from negative n to n, I get >an error message which says: "index must be a positive integer or logical."

Choose L=20 and the loop as

for j=1:n+1
    x(j)=-10+(j-1)*dx; %Because j actually starts at zero
    if x(j)> 0
        u(j,1)=1;
    else
        u(j,1)=0;
    end

end

>About the boundary values. Are you sure I am "allowed" to choose boundary values 0 at x=-10 and 1 at x=10? Since I am only given two initial conditions from the text, but no boundary conditions when t not equal to zero.

You will have to choose the x-interval long enough in positive and negative direction such that during the time of integration, the solution at the boundary remains constant. And don't forget to set the boundary values for u in time step m:

for m=1:maxm % Time Loop 
u(1,m)=0;
u(n+1,m)=1;
for j=2:n
...

Best wishes

Torsten.

[David]

Okay, thank you again. This helped a lot. Is it possible to choose u(1,m)=0 while using an x-value which is for example -10? (Because it is not supposed to be zero at x=0 when m not equal to zero). I have now changed the matlab-code with the changes in your post, do you believe it now is correct with respect to the heat equation with the initial conditions? (Is it for example okay to use the "if-method"? )

David

[Torsten]

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]

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

Answer 1

[Titus Edelhofer]

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

[David]

Thank you again for your reply. This looked very nice. When I write this in Matlab, how can I get this in without being replaced by the other u-values when using the for loop as written in the matlab-code on the top?

Does u=zeros(n+1,1) mean that it is only for the initial values?

[Titus Edelhofer]

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]

Thank you for your help. It is much appreciated.

David

Answer 2

[Titus Edelhofer]

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

[David]

Okay, thank you. Is the code correct according to the fact that the initial condition was that u=0 when t=0 and x>0, and u=1 when t=0 and x<0? The code does not feel correct to me. Should I replace the if-conditions, and use something else. It is only an initial condition.

[Titus Edelhofer]

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]

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.