Revisión del 19:38 18 may 2014

In this work we have studied the modeling of the heat equation, according to Fourier's law discovered in the nineteenth century.

1 Well proposed problem

Thin rod of length L

We will raise the system of equations that satisfies [math]u(x,t)[/math] assuming that the temperature of the rod [math]u(x,t)[/math] satisfies the heat equation [math]u_t-u_{xx}=0[/math]. First, we have a thin, homogeneous and thermally isolated by its lateral surface rod of length [math]L=3[/math]. At its left end the rod is in contact with a material whose temperature is maintained at 0°C, while the right is in contact with the material at 10°C. We also know that at the initial moment, the temperature distribution follows the [math]u(x,0)=u_0(x)[/math] function specified below. Assuming a standard [math]c=\rho=k=1[/math] parameters and there are no heat sources or sinks along the rod, the problem we have to solve is:

[math]
(P)
\begin{cases}
u_t-u_{xx}=0, \qquad x\epsilon(0,3), t\gt0\\
u(0,t)=0, u(3,t)=10, \qquad t\gt0\\
u(x,0)=u_0(x), \qquad x\epsilon[0,3]
\end{cases}
[/math]

[math] u_0(x)= \begin{cases} 10x/3 & \mbox{si} & x \in \mbox{(0,1)}\cup\mbox{(2,3)} \\ 100 & \mbox{si} & x \in \mbox{(1,2)} \end{cases} [/math]

Then we will define what is a well proposed problem is one that meets the following:

• Existence: problem [math](P)[/math] admits a solution [math]u(x,t)[/math].

• Uniqueness: if there is a solution [math](P)[/math] it has to be unique.

• Stability with respect to initial data: We consider the problem: [math] (P_1) \begin{cases} u_t-u_{xx}=0, \qquad x\epsilon(0,3), t\gt0\\ u(0,t)=0, u(3,t)=10, \qquad t\gt0\\ u(x,0)=h(x), \qquad x\epsilon[0,3] \end{cases} [/math]

Be [math]u(x,t)[/math] and [math]u_1(x,t)[/math] [math](P)[/math] and [math](P_1)[/math] solutions respectively. We say that the [math](P)[/math] problem is stable with respect to initial data if we prove the inequality of type [math]sup_{(x,t)\in [0,3]x[0,\infty]}\left|u(x,t)-u_1(x,t)\right|\leq Csup_{(x,t)\in [0,3]x[0,\infty]}\left|u_0(x)-h(x)\right|[/math] with [math]C[/math] absolute, constant independent of the [math](P)[/math] and [math](P_1)[/math] problems.

That [math] (P) [/math] is stable with respect to initial data tells us that if [math] h (x) [/math] is close to [math] u_0 (x) [/math] in the sense that [math] sup_ {(x, t) \ in [0,3] x [0, \infty]} \left | u_0 (x)-h (x) \right | [/math] is small, then the [math] u (x, t) [/math] and [math] u_1 (x, t) [/math] solutions are also nearby.

1.1 Resolution establishing finite difference method

Then the MATLAB code that numerically solves the heat equation posed exposed. It has been solved by the finite difference method with [math] \Delta x = 0.1 [/math] and we have used the method of taking time trapeze [math] \Delta t = \Delta x [/math]. The number of subintervals in which we divide the rod length is [math] Nx = 30 [/math] and time to which we have taken to represent [math] 2 [/math].

clear all
%solve the heat equation ut-uxx = 0 with
% u(0,t)=0 u(l,t)=10 u(x,0)=the piecewise title
L=3;T=2;
Nx=30;hx=L/Nx;
ht=hx;
K=(2*diag(ones(1,Nx-1))-diag(ones(1,Nx-2),1)-diag(ones(1,Nx-2),-1))/(hx^2);
x=hx:hx:L-hx;
F=zeros(Nx-1,1);F(Nx-1)=10/(hx^2);
t=0:ht:T;
u0=10*x/3;u0(10:20)=100;
uu=u0';
U(1,:)=[0 u0 10];
%method of trapezoids
for n=1:(length(t)-1)
 uu=(eye(Nx-1)+(ht/2)*K)\(uu+ht*(-K*uu+F+F)/2);
 U(n+1,:)=[0 uu' 10];
end
X=0:hx:L;
[xx,tt]=meshgrid(X,t);
m=U(:,16);
figure(1)
surf(xx,tt,U)
xlabel('space')
ylabel('time')
zlabel('temperature')
figure(2)
plot(t,m)
xlabel('time')
ylabel('temperature')

In this graph we have shown the 3D surface of the solution of the heat equation posed. As shown, although the trapezoidal method is an implicit method, not well approximated by the points of discontinuity of the initial condition and require less in the discretization step to remove these "peaks" that appear on the surface.

Solution of the heat equation

In the graph shown below the temperature behavior is shown in the middle of the rod with time. This is also obtained from the upper MATLAB code. Comparing the two graphs shows that the latter is a cut [math] x = 1.5 [/math] of the above.

u(t)

2 Resolution with different methods

We will solve the problem initially posed by the implicit and explicit methods by Euler and Runge-Kutta of order 4, following the same steps as with the method of the Trapezium.

2.1 Implicit Euler

clear all
%solve the heat equation ut-uxx = 0 with
% u(0,t)=0 u(l,t)=10 u(x,0)=the piecewise title
L=3;T=2;
Nx=30;hx=L/Nx;
ht=hx;
K=(2*diag(ones(1,Nx-1))-diag(ones(1,Nx-2),1)-diag(ones(1,Nx-2),-1))/(hx^2);
x=hx:hx:L-hx;
F=zeros(Nx-1,1);F(Nx-1)=10/(hx^2);
t=0:ht:T;
u0=10*x/3;u0(10:20)=100;
uu=u0';
U(1,:)=[0 u0 10];
%implicit Euler method
for n=1:(length(t)-1)
 uu=(eye(Nx-1)+(ht)*K)\(uu+ht*F);
 U(n+1,:)=[0 uu' 10];
end
X=0:hx:L;
[xx,tt]=meshgrid(X,t);
m=U(:,16);
figure(1)
surf(xx,tt,U)
xlabel('space')
ylabel('time')
zlabel('temperature')
figure(2)
plot(t,m)

Surface with the implicit Euler method

[math] x = 1.5 [/math] Graph with the implicit Euler method

2.2 Explicit Euler

clear all
% Explicit Euler method
L=3;T=2;
Nx=30;hx=L/Nx;
ht=(hx^2)/2;% must do so, if not it does not work
K=(2*diag(ones(1,Nx-1))-diag(ones(1,Nx-2),1)-diag(ones(1,Nx-2),-1))/(hx^2);
x=hx:hx:L-hx;
F=zeros(Nx-1,1);F(Nx-1)=10/(hx^2);
t=0:ht:T;
u0=10*x/3;u0(10:20)=100;
uu=u0';
U(1,:)=[0 u0 10];
for n=1:(length(t)-1)
 uu=uu+ht*(-K*uu+F);
 U(n+1,:)=[0 uu' 10];
end
X=0:hx:L;
[xx,tt]=meshgrid(X,t);
m=U(:,16);
figure(3)
surf(xx,tt,U)
xlabel('space')
ylabel('time')
zlabel('Temperature')
figure(4)
plot(t,m)

Surface with explicit Euler method

Graph in [math]x=1.5[/math] with explicit Euler method

2.3 Runge-Kutta

clear all
% Runge Kutta method
L=3;T=2;
Nx=30;hx=L/Nx;
ht=(hx^2)/2;% if we do not add ht it does not work
K=(2*diag(ones(1,Nx-1))-diag(ones(1,Nx-2),1)-diag(ones(1,Nx-2),-1))/(hx^2);
x=hx:hx:L-hx;
F=zeros(Nx-1,1);F(Nx-1)=10/(hx^2);
t=0:ht:T;
u0=10*x/3;u0(10:20)=100;
uu=u0';
U(1,:)=[0 u0 10];
for n=1:(length(t)-1)
 k1=-K*uu+F;
 k2=-K*(uu+k1*ht/2)+F;
 k3=-K*(uu+k2*ht/2)+F;
 k4=-K*(uu+k3*ht)+F;
 uu=uu+(ht/6)*(k1+2*k2+2*k3+k4);
 U(n+1,:)=[0 uu' 10];
end
X=0:hx:L;
[xx,tt]=meshgrid(X,t);
m=U(:,16);
figure(5)
surf(xx,tt,U)
xlabel('space')
ylabel('time')
zlabel('Temperature')
figure(6)
plot(t,m)

Surface with Runge Kutta method

Graph in [math]x=1.5[/math] Runge Kutta method

We can see that the method that works best is the implicit Euler, whereas explicit Euler and Runge-Kutta, being explicit methods require a rodent control into smaller intervals, and still not a good approximation is achieved as can be seen in the graph of the explicit Euler method. Therefore, for what follows, we use implicit Euler method.

3 Stationary state

It is said that a physical system is in stationary state when its characteristics do not vary with time. In this section we address this stationary state, which consists of neglecting the time and see what happens to our problem without taking into account the time variable.

[math]u_t(x,t)\approx 0; \ u_{xx}=0; \ u_x=c_1(t) \ u=c_1(t)x+c_2(t)[/math]

Substituting the boundary conditions:

[math]u(0)=0; \ c_2(t)=0; \ u(3)=10; \ c_1(t)=\frac{10}{3}[/math]

Thus the stationary solution is [math]u(x,t)=\frac{10x}{3}[/math], which is related to the initial condition. It seems logical that once the temperature in the center of the rod has dissipated, the ends having this constant temperature varies linearly between the two.

clear all
%solve the heat equation ut-uxx = 0 with
% u(0,t)=0 u(l,t)=10 u(x,0)=the piecewise title
L=3;T=10;
Nx=30;hx=L/Nx;
ht=hx;
K=(2*diag(ones(1,Nx-1))-diag(ones(1,Nx-2),1)-diag(ones(1,Nx-2),-1))/(hx^2);
x=hx:hx:L-hx;
F=zeros(Nx-1,1);F(Nx-1)=10/(hx^2);
t=0:ht:T;
u0=10*x/3;u0(10:20)=100;
uu=u0';
U(1,:)=[0 u0 10];
%método de Euler implícito
for n=1:(length(t)-1)
 uu=(eye(Nx-1)+(ht)*K)\(uu+ht*F);
 U(n+1,:)=[0 uu' 10];
end
X=0:hx:L;
[xx,tt]=meshgrid(X,t);
%solución estacionaria u(x)=10/3*x
V=10*xx/3;
figure(1)
hold on
surf(xx,tt,U)
xlabel('space')
ylabel('time')
zlabel('Temperature')
h1=surf(xx,tt,V),set(h1,'FaceColor','magenta','FaceAlpha',0.5,'EdgeColor','w')
hold off
t0=U(1,:);t1=U(1/ht+1,:);t2=U(2/ht+1,:);t10=U(10/ht+1,:);
figure(2)
subplot(2,2,1)
hold on
plot(X,t0)
plot(X,V(1,:),'r')
xlabel('x')
ylabel('temperature with t=0')
hold off
subplot(2,2,2)
hold on
plot(X,t1)
plot(X,V(1+1/ht,:),'r')
xlabel('x')
ylabel('temperature with t=1')
hold off
subplot(2,2,3)
hold on
plot(X,t2)
plot(X,V(1+2/ht,:),'r')
xlabel('x')
ylabel('temperature with t=2')
hold off
subplot(2,2,4)
hold on
plot(X,t10)
plot(X,V(1+10/ht,:),'r')
xlabel('x')
ylabel('temperature with t=10')
hold off
figure(3)
subplot(2,2,1)
e1=abs(t0-V(1,:));me1=max(e1);...
sprintf('The maximum difference with the stationary solution in t=0 is %d .',me1)
plot(X,e1,'g'),xlabel('x'),...
ylabel('Difference with the stationary solution in t=0') 
subplot(2,2,2)
e2=abs(t1-V(1+1/ht,:));me2=max(e2);...
sprintf('The maximum difference with the stationary solution in t=1 is %d .',me2)
plot(X,e2,'g'),xlabel('x'),...
ylabel('Difference with the stationary solution t=1') 
subplot(2,2,3)
e3=abs(t2-V(1+2/ht,:));me3=max(e3);...
sprintf('The maximum difference with the stationary solution in t=2 is %d .',me3)
plot(X,e3,'g'),xlabel('x'),...
ylabel('Difference with the stationary solution in t=2') 
subplot(2,2,4)
e4=abs(t10-V(1+10/ht,:));me4=max(e4);...
sprintf('The maximum difference with the stationary solution in t=10 is %d .',me4)
plot(X,e4,'g'),xlabel('x'),ylabel('Difference with the stationary solution in t=10')

Real surfaces and the stationary solution (in pink)

This second graph shows as as we move in time (t older) solving our heat equation (blue) is more assimilated to the (red) stationary solution.

Comparing solutions in [math]t=0,1,2,10[/math]

We now show the difference between the previous two solutions stationary real and represented throughout the rod for different values of time. We see how to increasingly large time difference between the two is narrowing, observing the order of magnitude in the ordinate.

Differences with the stationary solution [math]t=0,1,2,10[/math] order

4 Neumann type boundary condition

Now let's consider a different boundary condition at the right end. Instead of assuming a constant temperature at that end as above, we will place on it an insulating piece. This isolate causes no loss of heat at the right end, that is, the flow temperature is null. This condition is of Neumann type, unlike the previous ones were Dirichlet. So, we keep the condition at the left end and apply the new on the far right, which is

In this situation, the temperature of the rod is given by the following problem

So far in the stationary state for large times the [math] function u (x, t) [/math] that models the temperature of the rod is solution of the following boundary value problem (we call it that because the differential equation depends only [math] x [/math] in the stationary state)

We solve the differential equation to obtain the same result as above

Applying the boundary conditions

This result shows that after a large enough time to consider a steady state in the rod, it acquires a uniform zero temperature. The behavior of the rod is consistent with the boundary conditions, as its final temperature matches that remains constant in the far left.

4.1 Finite difference method

The following image shows an approximation of the problem is shown by the method of finite differences. It can be seen as a high value of the temperature in the rod can be considered constant and uniform throughout, reaching the stationary value [math] u (x, t) = 0 [/math].

Solution of the problem with Neumann type boundary condition

Specifically, from a time [math] t = 26.4[/math] we can consider that the temperature reaches stationary value with an error of 0.05, that is, at that moment the difference between the calculated and the thermal distribution stationary takes that value.

Below is reflected Matlab code which approximates the temperature of the rod by using the finite difference method as the implicit Euler that provides a better approximation with a step size [math] h = 0.1 [/math] in time and space, and [math]t \in \mbox{[0,30]}[/math]. Furthermore, in the instant code approximation value differs 0.05 steady in all parts of the rod is calculated.

clear all
%Sixth paragraph of labor
%solve the heat equation ut-uxx = 0 with
% U (0, t) = 0 ux (L, t) = 0 u (x, 0) = the function piecewise title
L=3;T=30;
Nx=30;hx=L/Nx;
ht=2*hx;
K=(2*diag(ones(1,Nx))-diag(ones(1,Nx-1),1)-diag(ones(1,Nx-1),-1));
K(Nx,Nx-1)=-2;K=K/(hx^2);
x=hx:hx:L;
F=zeros(Nx,1);
t=0:ht:T;
u0=10*x/3;u0(10:20)=100;
uu=u0';
U(1,:)=[0 u0];
%metoodo de Euler implicito
for n=1:(length(t)-1)
 uu=(eye(Nx)+(ht)*K)\(uu+ht*F);
 U(n+1,:)=[0 uu'];
end
p=0.05*ones(1,length(x)+1);
for i=1:length(t)
if min(U(i,:)<=p)==1 
    break 
end
end
SOL=t(i)
X=0:hx:L;
[xx,tt]=meshgrid(X,t);
surf(xx,tt,U)
title('Solution of the problem with Neumann type boundary condition')
xlabel('Space')
ylabel('Time')
zlabel('Temperature')
p=U(:,Nx+1);

4.2 Fourier Method

We propose now the same problem using the Fourier method. Thus, we seek solutions [math] u (x, t) = \ varphi (x) T (t) [/math] form, where [math] \ varphi (x) [/math] must satisfy the following problem eigenvalue

The solution of the differential equation is [math]\varphi(x)=a\cos(\sqrt{\lambda}x)+b\sin(\sqrt{\lambda}x)[/math]

Applying the boundary conditions we obtain the eigenvalues [math]\mu_k[/math] and eigenfunctions [math]\varphi_k(x)[/math] of the problem

The solution of this differential equation $\boxed{T_k(t)=C_ke^{-\mu_kt}= C_ke^{-(k-{1\over2})^2{\pi^2\over9}t}}$ and therefore [math]u_k(x,t)= \varphi_k(x)T_k(t)=C_ke^{-(k-{1\over2})^2{\pi^2\over9}t}\sin(k-{1\over2}){\pi\over3}x[/math]

If we express the solution of the problem as $ \ boxed {u (x, t) = \sum_ {k = 1} ^ Nu_k (x, t) = \sum_ {k = 1} ^ N ^ C_ke {- (k-{1 \over2}) ^ 2 {\pi ^ 2 \over9} t} \sin (k-{1 \over2}) {\pi \over3} x} $ and make satisfying the initial condition we obtain [math]u(x,0)=\sum_{k=1}^NC_k\sin(k-{1\over2}){\pi\over3}x[/math]. Thus, by uniqueness of the Fourier coefficients, the coefficients [math] C_K [/math] match those of the Fourier series with respect to the eigenfunctions [math] \varphi_k (x) [/math] Function apart that determines the initial condition (expressed at the beginning of the article).

The problem is thus limited to the calculation of these coefficients according to the expression

The degree of accuracy of the approximation with this method depends on the number of elements in the Fourier series, ie, the value of [math] N [/math]. We study the temperature of the rod taking [math] N = 1,3,5,10,20 [/math] values, as you can see in the picture below. It can be seen, especially in the initial condition, that as the value of [math] N [/math] is the approximate increase function is closer to the real.

Solutions with different number of terms of the Fourier series

These results can be better taking a compare [math] t = 0.5 [/math] fixed instant, and representing each function it in the same graph, as we see below. For this case we also added the approximation with 2 terms of the Fourier series to reflect that up to 3 terms approaches can be distinguished, but once this number of elements in the series approximations are virtually indistinguishable when compared on the same graph.

Graph of temperature in t=0.5 with N terms of the Fourier series

The following Matlab code which approximates the temperature of the rod by the Fourier method, with N = 1,3,5,10,20 terms of the series, and step size is reflected [math]h=0,1[/math] in time and space, and [math]t \in \mbox{[0,10]}[/math]

clear all
%Seventh section of the paper
% solve the heat equation ut-uxx = 0 with
% U (0, t) = 0 ux (L, t) = 0 u (x, 0) = the function piecewise title
% Resolution with fourier
L=3;T=10;
Nx=30;hx=L/Nx;
ht=hx;
x=0:hx:L;t=0:ht:T;
[xx,tt]=meshgrid(x,t);
q=[1,2,3,5,10,20];
a=[1,3,5,10,20];
f=10*x/3;f(11:21)=100;
i=0;
for Q=q
u=0;
for k=1:Q
    p=sin((k-1/2)*(pi/3)*x);
    c=trapz(x,f.*p)/trapz(x,p.*p);
    u=u+c*exp(-((k-1/2)^2)*((pi/3)^2)*tt).*sin((k-1/2)*(pi/3)*xx);
end
b=find(Q==a);
if b<6
i=i+1;
figure(1)
subplot(2,3,i)
ca=num2str(Q);r=strcat(['Solution with ',ca,' terms of the Fourier series']);
if Q==1
r=strcat(['Solution with ',ca,' term of the Fourier series']);
end
surf(xx,tt,u),xlabel('Space'),ylabel('Time'),zlabel('Temperature'),title(r)
end
d(find(q==Q),:)=u(0.5/ht+1,:);
end
figure(2)
hold on
title('Temperature in t=0.5 with N terms of the Fourier series')
plot(x,d(1:6,:))
xlabel('Space'),ylabel('Time')
legend('N=1','N=2','N=3','N=5','N=10','N=20','Location','best')
hold off

5 Losses along the rod

Now we will study the case where there are heat sources or sinks along the rod. Specifically, if there is heat loss through the air having a constant temperature of 16 degrees. With the boundary conditions that we had initially keeping the left and right at 0 and 10 degrees respectively ends, the problem would be:

[math] \begin{cases} u_t-u_{xx}+u-16=0,\ x \in \mbox{[0,3]} & t\gt0 \\ u(0,t)=0; u(3,t)=10; & u(x,0)=g(x) ; \end{cases} [/math]

Where [math]g(x)[/math] is defined as in the first problem:

[math] u_0= \begin{cases} 10x/3 & \mbox{si} & x \in \mbox{(0,1)} \cup \mbox{(2,3)} \\ 100 & \mbox{si} & x \in \mbox{(1,2)} \end{cases} [/math]

To solve the problem by finite differences with a term in [math] u (x, t) [/math] need to rethink the discretization in space:

[math] \begin{cases} u’_n(t)+\frac{-u_{n-1}(t)+2u_n(t)-u_{n+1}(t)}{h^2}+u_n(t)=16&n=1,2,…,N&t\gt0 \\ u_0(t)=0 \\ u_N(t)=10 \\ u_n(0)=g(x) \end{cases} [/math]

Thus, the resulting matrices are as follows

[math] K= \begin{pmatrix} 2 & -1 & 0 & 0 & … & 0 & 0 \\ -1 & 2 & 0 & 0 & … &0 & 0 \\ … & … & … & … & … & … & … \\ 0 & 0 & 0 & 0 & … & -1 & 2 \end{pmatrix} \frac{1}{h^2}+ \begin{pmatrix} 1 & 0 & 0 & 0 & … & 0 & 0 \\ 0 & 1 & 0 & 0 & … &0 & 0 \\ … & … & … & … & … & … & … \\ 0 & 0 & 0 & 0 & … & 0 & 1 \end{pmatrix} [/math]; [math] F= \begin{pmatrix} 16+\frac{0}{h^2} \\ 16 \\ … \\ 16+\frac{10}{h^2} \\ \end{pmatrix} [/math]; [math] U= \begin{pmatrix} u_1 \\ u_2 \\ … \\ u_{N-1} \\ \end{pmatrix} [/math]; [math] U^0= \begin{pmatrix} g(x_1) \\ g(x_2) \\ … \\ g(x_{N-1}) \end{pmatrix} [/math]

With this we can move to numerically solve the problem, but we must recalculate the steady state of the despising rod [math]u_t(x,t)[/math].

[math]u_t(x,t)\approx 0; \ u_{xx}-u+16=0; \ u(0)=0 \ u(3)=10[/math]

This is nothing more than an ordinary differential equation of 2nd order with constant coefficients, whose solution is:

[math]a(t)e^x+b(t)e^{-x}+16[/math]

Where [math] a (t) [/math] and [math] b (t) [/math] are constants to be obtained to replace and solve the system with the boundary conditions, which in our case have let him Matlab resolved (see lines 22 and 23 of the code). Now, we turn to numerically solve the problem, with [math]h=0,1[/math] and [math]t \in \mbox{[0,10]}[/math]:

clear all
%solve the heat equation ut-uxx + u-16 = 0 with
% U (0, t) = 0 u (l, t) = 10 u (x, 0) = the function to set pieces
L=3;T=10;
Nx=30;hx=L/Nx;
ht=hx;
K=(2*diag(ones(1,Nx-1))-diag(ones(1,Nx-2),1)-diag(ones(1,Nx-2),-1))/(hx^2);
K=K+eye(Nx-1);
x=hx:hx:L-hx;
F=16*ones(Nx-1,1);F(Nx-1)=16+10/(hx^2);
t=0:ht:T;
u0=10*x/3;u0(10:20)=100;
uu=u0';
U(1,:)=[0 u0 10];
%metoodo de Euler implicito
for n=1:(length(t)-1)
 uu=(eye(Nx-1)+(ht)*K)\(uu+ht*F);
 U(n+1,:)=[0 uu' 10];
end
X=0:hx:L;
[xx,tt]=meshgrid(X,t);
a=[1 1;exp(3) exp(-3)];b=[-16;-6];d=a\b;
V=d(1)*exp(xx)+d(2)*exp(-xx)+16*ones(T/ht+1,Nx+1);
figure(1)
hold on
surf(xx,tt,U)
xlabel('space')
ylabel('time')
zlabel('temperature')
h1=surf(xx,tt,V),set(h1,'FaceColor','red','FaceAlpha',0.5,'EdgeColor','w')
xlabel('space')
ylabel('time')
zlabel('temperature')
E=abs(U-V);
e1=E(46:101,Nx/6+1);
e2=E(46:101,Nx/3+1);
e3=E(46:101,Nx/2+1);
e4=E(46:101,2*Nx/3+1);
e5=E(46:101,5*Nx/6+1);
figure(3)
hold on
plot(t(46:101),e1)
plot(t(46:101),e2,'r')
plot(t(46:101),e3,'k')
plot(t(46:101),e4,'m')
plot(t(46:101),e5,'c'),legend('Difference in x=0.5','Difference in x=1',...
'Difference in x=1.5','Difference in x=2','Difference in x=2.5')
plot(t(46:101),0.001*ones(1,length(t)-45),'g'),xlabel('time'),...
ylabel('Difference with the stationary solution')
hold off

Graph of [math]u(x,t)[/math] in [math]x \in \mbox{[0,3]};t \in \mbox{[0,10]}[/math]. In pink is the stationary solution.

As shown in the graph, the actual solution and stationary are virtually identical for [math] t\gt 2 [/math]. Furthermore we see that near [math] t = 5 [/math] the stationary solution is above the real solution. This is best seen in the following graph in which are represented the difference between the actual solution and the stationary for different values of [math] x [/math], namely [math]x = 0.5,1,1.5,2,2.5 [/math].

Graph with the difference between the real and the stationary solution for[math]t \in \mbox{[4.5,10]}[/math]

At that point cut the difference between the two is minimal, below [math] 10 ^ {-3 } [/math], while from there this difference increases slightly and remained constant when time is increasing. This was expected, since for large times the solution of the equation is stationary, and therefore effects of time we have is the difference between two constants. The existing small difference between the two may be because logically, does not disclose exact thing but we are solving the equation numerically.

5.1 Changing the boundary conditions

Suppose now that we change the boundary conditions, so that now the left end of the rod to be in contact with a material whose temperature varies according to the function [math] 10sen (t) [/math], and by the end right there is a flow of constant heat input [math] 1 [/math]. These conditions translate as [math] u (0, t) = 10sen (t) [/math] and [math] u_x (3, t) = 1 [/math].

As the second condition of Neumann type, size of the matrices increases one unit, to be known [math] u_n [/math] term. Also [math] K [/math] and [math] F [/math] change their terms, becoming

[math] K= \begin{pmatrix} 2 & -1 & 0 & 0 & … & 0 & 0 \\ -1 & 2 & 0 & 0 & … &0 & 0 \\ … & … & … & … & … & … & … \\ 0 & 0 & 0 & 0 & … & -2 & 2 \end{pmatrix} \frac{1}{h^2}+ \begin{pmatrix} 1 & 0 & 0 & 0 & … & 0 & 0 \\ 0 & 1 & 0 & 0 & … &0 & 0 \\ … & … & … & … & … & … & … \\ 0 & 0 & 0 & 0 & … & 0 & 1 \end{pmatrix} [/math]; [math] F= \begin{pmatrix} 16+\frac{10sen(t_n)}{h^2} \\ 16 \\ … \\ 16+\frac{2}{h} \\ \end{pmatrix} [/math]

We solve the implicit Euler method with [math]h=0,1[/math] and [math]t \in \mbox{[0,10]}[/math]:

clear all
%solve the heat equation ut-uxx + u-16 = 0 with
% U (0, t) = 10 * sin (t) ux (L, t) = 1 u (x, 0) = the function piecewise title
L=3;T=10;
Nx=30;hx=L/Nx;
ht=hx;
K=(2*diag(ones(1,Nx))-diag(ones(1,Nx-1),1)-diag(ones(1,Nx-1),-1));K(Nx,Nx-1)=-2;K=K/(hx^2);
K=K+eye(Nx);
x=hx:hx:L;
t=0:ht:T;
u0=10*x/3;u0(10:20)=100;
uu=u0';
g=10*sin(t);
U(1,:)=[g(1) u0];
%implicit Euler method
for n=1:(length(t)-1)
   F=16*ones(Nx,1);F(Nx)=16+2/hx;F(1)=16+g(n)/(hx^2);
   uu=(eye(Nx)+(ht)*K)\(uu+ht*F);
   U(n+1,:)=[g(n+1) uu'];
end
X=0:hx:L;
[xx,tt]=meshgrid(X,t);
figure(1)
surf(xx,tt,U)
xlabel('space')
ylabel('time')
zlabel('Temperature')

Graph of [math]u(x,t)[/math] in [math]x \in \mbox{[0,3]};t \in \mbox{[0,10]}[/math]

In view of the graph we can conclude that as previously when we had a Neumann type condition, the heat "escapes" from the left end at a variable temperature with time while entering from the right. This can be seen if the program is run on the graph by observing the slight slope that the surface of temperatures has at the left side.