MateWiki - Contribuciones del usuario [es]

Heat equation (Grupo 1B)

2014-05-20T15:42:10Z

Marino Rivera:

{{ TrabajoED |Heat equation. Grupo 1-B | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandro Andrés Martínez

David Ayala Díez

Claudia Cózar Coarasa

Lorena de la Fuente Sanz

Marino Rivera Muñoz

José Manuel Torres Serrano }}

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

= Well proposed problem =

[[Archivo:Nueva imagen.png|thumb|300px|left| 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>0\\
u(0,t)=0, u(3,t)=10, \qquad t>0\\
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>0\\
u(0,t)=0, u(3,t)=10, \qquad t>0\\
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.

==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>.

{{matlab|codigo=
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.

{|
|-
| [[Archivo: untitled11.jpg|thumb|750px|left|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.

{|
|-
| [[Archivo: gráfica3a.jpg|thumb|750px|left| u(t) ]]
|}

= 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.

==Implicit Euler==

{{matlab|codigo=

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)

}}

{|
|-
|[[Archivo: untitled00010.jpg|thumb|750px|left|Surface with the implicit Euler method]]

|}

{|
|-
|[[Archivo: Grafica_Euler_impl%C3%ADcito.jpg |thumb|750px|left|<math> x = 1.5 </math> Graph with the implicit Euler method]]

|}

==Explicit Euler==

{{matlab|codigo=

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)

}}

{|
|-
|[[Archivo: untitled9.jpg|thumb|750px|left|Surface with explicit Euler method ]]

|}

{|
|-
|[[Archivo: Euler1.3.jpg |thumb|750px|left|Graph in <math>x=1.5</math> with explicit Euler method]]

|}

==Runge-Kutta==

{{matlab|codigo=

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)

}}

{|
|-
|[[Archivo: untitled8.jpg|thumb|750px|left|Surface with Runge Kutta method]]

|}

{|
|-
|[[Archivo: Grafica_Runge_Kutta.jpeg|thumb|750px|left|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.

=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.

{{matlab|codigo=

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')

}}

{|
|-
|[[Archivo: untitled7.jpg|thumb|750px|left|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.

{|
|-
|[[Archivo: untitled5.jpg|thumb|750px|left|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.

{|
|-
|[[Archivo: untitled6.jpg|thumb|750px|left|Differences with the stationary solution <math>t=0,1,2,10</math> order]]

|}

= 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

<center><math>-ku_x(3,t)=0 \rightarrow </math> $\boxed{u_x(3,t)=0}$</center>

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

<center><math>
\begin{cases}
u_t-u_{xx}=0, \qquad x\epsilon(0,3), t>0\\
u(0,t)=0, u_x(3,t)=0, \qquad t>0\\
u(x,0)=u_0(x), \qquad x\epsilon[0,3]
\end{cases}
</math></center>

So far in the stationary state for large times the function <math>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)

<center><math>
\begin{cases}
u_{xx}=0, \qquad x\epsilon(0,3), t>0\\
u(0,t)=0, u_x(3,t)=0, \qquad t>0
\end{cases}
</math></center>

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

<center><math>u_{xx}=0 \rightarrow u(x,t)=C_1(t)x+C_2(t)</math></center>

Applying the boundary conditions

<center><math>C_1(t)=C_2(t)=0 \rightarrow </math> $\boxed{u(x,t)=0}$</center>

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.

==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>.

{|
|-
|[[Archivo:AP6a.png|thumb|750px|left|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.

{{matlab|codigo=
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);
}}

==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

<center><math>
\begin{cases}
\varphi’’(x)+\lambda\varphi(x)=0\\
\varphi(0)=0, \varphi’(3)=0
\end{cases}
</math></center>

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
<center><math>\varphi(0)=0 \rightarrow a=0</math>;</center>
<center><math>\varphi’(3)=0 \rightarrow \sqrt{\lambda}\,b\cos(\sqrt{\lambda}\,3)=0 \rightarrow \sqrt{\lambda}\,3=(k-{1\over2})\,\pi \rightarrow </math> $\boxed{\lambda=\mu_k=(k-{1\over2})^2{\pi^2\over9}}$;</center>

<center>So $\boxed{\varphi_k(x)=\sin(k-{1\over2}){\pi\over3}x}$, <math> \qquad k=1,2,3…N</math></center>

If we apply that <math>u_k(x,t)=\varphi_k(x)T_k(t)</math> satisfies the differential equation <math> u_t-u_{xx}=0</math>, we obtain the differential equation determines <math>T_k(t)</math> <center> <math>u_t-u_{xx}=0 \rightarrow \varphi_k(x)T_k’(t)-\varphi_k’’(x)T_k(t)=\varphi_k(x)T_k’(t)-(-\mu_k)\varphi_k(x)T_k(t)=0 \rightarrow</math> $\boxed{T_k’(t)+\mu_kT_k(t)=0}$</center>

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

<center>$\boxed{C_k={\int_{0}^{3}u(x,0)\varphi_k(x)dx\over\int_{0}^{3}\varphi_k^2(x)dx}}$</center>

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.
{|
|-
|[[Archivo:Untitled4.jpg|thumb|750px|left|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.

{|
|-
| [[Archivo:Untitled3.jpg|thumb|750px|left|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>

{{matlab|codigo=
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
}}

=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>0 \\
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>0 \\
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>:

{{matlab|codigo=
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

}}

{|
|-
|[[Archivo:Untitled002.jpg|thumb|750px|left|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> 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>.

{|
|-
|[[Archivo:Dif7.jpg|thumb|750px|left|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.

== 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>:

{{matlab|codigo=
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')

}}

{|
|-
|[[Archivo:untitled.jpg|thumb|750px|left|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.

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Heat equation (Grupo 1B)

2014-05-20T15:39:12Z

Marino Rivera:

{{ TrabajoED |Heat equation. Grupo 1-B | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandro Andrés Martínez

David Ayala Díez

Claudia Cózar Coarasa

Lorena de la Fuente Sanz

Marino Rivera Muñoz

José Manuel Torres Serrano }}

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

= Well proposed problem =

[[Archivo:Nueva imagen.png|thumb|300px|left| 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>0\\
u(0,t)=0, u(3,t)=10, \qquad t>0\\
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>0\\
u(0,t)=0, u(3,t)=10, \qquad t>0\\
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.

==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>.

{{matlab|codigo=
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.

{|
|-
| [[Archivo: untitled11.jpg|thumb|750px|left|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.

{|
|-
| [[Archivo: gráfica3a.jpg|thumb|750px|left| u(t) ]]
|}

= 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.

==Implicit Euler==

{{matlab|codigo=

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)

}}

{|
|-
|[[Archivo: untitled00010.jpg|thumb|750px|left|Surface with the implicit Euler method]]

|}

{|
|-
|[[Archivo: Grafica_Euler_impl%C3%ADcito.jpg |thumb|750px|left|<math> x = 1.5 </math> Graph with the implicit Euler method]]

|}

==Explicit Euler==

{{matlab|codigo=

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)

}}

{|
|-
|[[Archivo: untitled9.jpg|thumb|750px|left|Surface with explicit Euler method ]]

|}

{|
|-
|[[Archivo: Euler1.3.jpg |thumb|750px|left|Graph in <math>x=1.5</math> with explicit Euler method]]

|}

==Runge-Kutta==

{{matlab|codigo=

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)

}}

{|
|-
|[[Archivo: untitled8.jpg|thumb|750px|left|Surface with Runge Kutta method]]

|}

{|
|-
|[[Archivo: Grafica_Runge_Kutta.jpeg|thumb|750px|left|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.

=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.

{{matlab|codigo=

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')

}}

{|
|-
|[[Archivo: untitled7.jpg|thumb|750px|left|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.

{|
|-
|[[Archivo: untitled5.jpg|thumb|750px|left|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.

{|
|-
|[[Archivo: untitled6.jpg|thumb|750px|left|Differences with the stationary solution <math>t=0,1,2,10</math> order]]

|}

= 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

<center><math>-ku_x(3,t)=0 \rightarrow </math> $\boxed{u_x(3,t)=0}$</center>

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

<center><math>
\begin{cases}
u_t-u_{xx}=0, \qquad x\epsilon(0,3), t>0\\
u(0,t)=0, u_x(3,t)=0, \qquad t>0\\
u(x,0)=u_0(x), \qquad x\epsilon[0,3]
\end{cases}
</math></center>

So far in the stationary state for large times the function <math>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)

<center><math>
\begin{cases}
u_{xx}=0, \qquad x\epsilon(0,3), t>0\\
u(0,t)=0, u_x(3,t)=0, \qquad t>0
\end{cases}
</math></center>

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

<center><math>u_{xx}=0 \rightarrow u(x,t)=C_1(t)x+C_2(t)</math></center>

Applying the boundary conditions

<center><math>C_1(t)=C_2(t)=0 \rightarrow </math> $\boxed{u(x,t)=0}$</center>

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.

==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>.

{|
|-
|[[Archivo:AP6a.png|thumb|750px|left|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.

{{matlab|codigo=
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);
}}

==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

<center><math>
\begin{cases}
\varphi’’(x)+\lambda\varphi(x)=0\\
\varphi(0)=0, \varphi’(3)=0
\end{cases}
</math></center>

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
<center><math>\varphi(0)=0 \rightarrow a=0</math>;</center>
<center><math>\varphi’(3)=0 \rightarrow \sqrt{\lambda}\,b\cos(\sqrt{\lambda}\,3)=0 \rightarrow \sqrt{\lambda}\,3=(k-{1\over2})\,\pi \rightarrow </math> $\boxed{\lambda=\mu_k=(k-{1\over2})^2{\pi^2\over9}}$;</center>

<center>So $\boxed{\varphi_k(x)=\sin(k-{1\over2}){\pi\over3}x}$, <math> \qquad k=1,2,3…N</math></center>

If we apply that <math>u_k(x,t)=\varphi_k(x)T_k(t)</math> satisfies the differential equation <math> u_t-u_{xx}=0</math>, we obtain the differential equation determines <math>T_k(t)</math> <center> <math>u_t-u_{xx}=0 \rightarrow \varphi_k(x)T_k’(t)-\varphi_k’’(x)T_k(t)=\varphi_k(x)T_k’(t)-(-\mu_k)\varphi_k(x)T_k(t)=0 \rightarrow</math> $\boxed{T_k’(t)+\mu_kT_k(t)=0}$</center>

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

<center>$\boxed{C_k={\int_{0}^{3}u(x,0)\varphi_k(x)dx\over\int_{0}^{3}\varphi_k^2(x)dx}}$</center>

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.
{|
|-
|[[Archivo:Untitled4.jpg|thumb|750px|left|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.

{|
|-
| [[Archivo:Untitled3.jpg|thumb|750px|left|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>

{{matlab|codigo=
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
}}

=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>0 \\
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>0 \\
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>:

{{matlab|codigo=
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

}}

{|
|-
|[[Archivo:Untitled002.jpg|thumb|750px|left|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> 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>.

{|
|-
|[[Archivo:Untitled01.jpg|thumb|750px|left|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.

== 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>:

{{matlab|codigo=
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')

}}

{|
|-
|[[Archivo:untitled.jpg|thumb|750px|left|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.

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Heat equation (Grupo 1B)

2014-05-18T17:17:08Z

Marino Rivera:

{{ TrabajoED |Heat equation. Grupo 1-B | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandro Andrés Martínez

David Ayala Díez

Claudia Cózar Coarasa

Lorena de la Fuente Sanz

Marino Rivera Muñoz

José Manuel Torres Serrano }}

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

= Well proposed problem =

[[Archivo:Nueva imagen.png|thumb|300px|left| 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>0\\
u(0,t)=0, u(3,t)=10, \qquad t>0\\
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>0\\
u(0,t)=0, u(3,t)=10, \qquad t>0\\
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.

==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>.

{{matlab|codigo=
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.

{|
|-
| [[Archivo: untitled11.jpg|thumb|750px|left|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.

{|
|-
| [[Archivo: gráfica3a.jpg|thumb|750px|left| u(t) ]]
|}

= 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.

==Implicit Euler==

{{matlab|codigo=

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)

}}

{|
|-
|[[Archivo: untitled00010.jpg|thumb|750px|left|Surface with the implicit Euler method]]

|}

{|
|-
|[[Archivo: Grafica_Euler_impl%C3%ADcito.jpg |thumb|750px|left|<math> x = 1.5 </math> Graph with the implicit Euler method]]

|}

==Explicit Euler==

{{matlab|codigo=

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)

}}

{|
|-
|[[Archivo: untitled9.jpg|thumb|750px|left|Surface with explicit Euler method ]]

|}

{|
|-
|[[Archivo: Euler1.3.jpg |thumb|750px|left|Graph in <math>x=1.5</math> with explicit Euler method]]

|}

==Runge-Kutta==

{{matlab|codigo=

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)

}}

{|
|-
|[[Archivo: untitled8.jpg|thumb|750px|left|Surface with Runge Kutta method]]

|}

{|
|-
|[[Archivo: Grafica_Runge_Kutta.jpeg|thumb|750px|left|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.

=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.

{{matlab|codigo=

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')

}}

{|
|-
|[[Archivo: untitled7.jpg|thumb|750px|left|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.

{|
|-
|[[Archivo: untitled5.jpg|thumb|750px|left|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.

{|
|-
|[[Archivo: untitled6.jpg|thumb|750px|left|Differences with the stationary solution <math>t=0,1,2,10</math> order]]

|}

= 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

<center><math>-ku_x(3,t)=0 \rightarrow </math> $\boxed{u_x(3,t)=0}$</center>

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

<center><math>
\begin{cases}
u_t-u_{xx}=0, \qquad x\epsilon(0,3), t>0\\
u(0,t)=0, u_x(3,t)=0, \qquad t>0\\
u(x,0)=u_0(x), \qquad x\epsilon[0,3]
\end{cases}
</math></center>

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)

<center><math>
\begin{cases}
u_{xx}=0, \qquad x\epsilon(0,3), t>0\\
u(0,t)=0, u_x(3,t)=0, \qquad t>0
\end{cases}
</math></center>

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

<center><math>u_{xx}=0 \rightarrow u(x,t)=C_1(t)x+C_2(t)</math></center>

Applying the boundary conditions

<center><math>C_1(t)=C_2(t)=0 \rightarrow </math> $\boxed{u(x,t)=0}$</center>

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.

==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>.

{|
|-
|[[Archivo:AP6a.png|thumb|750px|left|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.

{{matlab|codigo=
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);
}}

==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

<center><math>
\begin{cases}
\varphi’’(x)+\lambda\varphi(x)=0\\
\varphi(0)=0, \varphi’(3)=0
\end{cases}
</math></center>

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
<center><math>\varphi(0)=0 \rightarrow a=0</math>;</center>
<center><math>\varphi’(3)=0 \rightarrow \sqrt{\lambda}\,b\cos(\sqrt{\lambda}\,3)=0 \rightarrow \sqrt{\lambda}\,3=(k-{1\over2})\,\pi \rightarrow </math> $\boxed{\lambda=\mu_k=(k-{1\over2})^2{\pi^2\over9}}$;</center>

<center>Así $\boxed{\varphi_k(x)=\sin(k-{1\over2}){\pi\over3}x}$, <math> \qquad k=1,2,3…N</math></center>

If we apply that <math> u_k (x, t) = \ varphi_k (x) T_k (t) </ math> satisfies the differential equation<math> u_t-u {xx} = 0 </ math>, we obtain the differential equation determines <math>T_k(t)</math><center><math>u_t-u{xx}=0 \rightarrow \varphi_k(x)T_k’(t)-\varphi_k’’(x)T_k(t)=\varphi_k(x)T_k’(t)-(-\mu_k)\varphi_k(x)T_k(t)=0 \rightarrow </math> $\boxed{T_k’(t)+\mu_kT_k(t)=0}$</center>

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

<center>$\boxed{C_k={\int_{0}^{3}u(x,0)\varphi_k(x)dx\over\int_{0}^{3}\varphi_k^2(x)dx}}$</center>

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.
{|
|-
|[[Archivo:Untitled4.jpg|thumb|750px|left|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.

{|
|-
| [[Archivo:Untitled3.jpg|thumb|750px|left|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>

{{matlab|codigo=
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
}}

=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>0 \\
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>0 \\
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>:

{{matlab|codigo=
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

}}

{|
|-
|[[Archivo:Untitled002.jpg|thumb|750px|left|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> 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>.

{|
|-
|[[Archivo:Untitled01.jpg|thumb|750px|left|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.

== 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>:

{{matlab|codigo=
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')

}}

{|
|-
|[[Archivo:untitled.jpg|thumb|750px|left|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.

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Heat equation (Grupo 1B)

2014-05-18T17:11:07Z

Marino Rivera: Marino Rivera movió la página Heat equation a Heat equation (Grupo 1B)

{{ TrabajoED |Heat equation. Grupo 1-B | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandro Andrés Martínez

David Ayala Díez

Claudia Cózar Coarasa

Lorena de la Fuente Sanz

Marino Rivera Muñoz

José Manuel Torres Serrano }}

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

= Well proposed problem =

[[Archivo:Nueva imagen.png|thumb|300px|left| 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>0\\
u(0,t)=0, u(3,t)=10, \qquad t>0\\
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>0\\
u(0,t)=0, u(3,t)=10, \qquad t>0\\
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.

==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>.

{{matlab|codigo=
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.

{|
|-
| [[Archivo: untitled11.jpg|thumb|750px|left|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.

{|
|-
| [[Archivo: gráfica3a.jpg|thumb|750px|left| u(t) ]]
|}

= 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.

==Implicit Euler==

{{matlab|codigo=

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)

}}

{|
|-
|[[Archivo: untitled00010.jpg|thumb|750px|left|Surface with the implicit Euler method]]

|}

{|
|-
|[[Archivo: Grafica_Euler_impl%C3%ADcito.jpg |thumb|750px|left|<math> x = 1.5 </math> Graph with the implicit Euler method]]

|}

==Explicit Euler==

{{matlab|codigo=

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)

}}

{|
|-
|[[Archivo: untitled9.jpg|thumb|750px|left|Surface with explicit Euler method ]]

|}

{|
|-
|[[Archivo: Euler1.3.jpg |thumb|750px|left|Graph in <math>x=1.5</math> with explicit Euler method]]

|}

==Runge-Kutta==

{{matlab|codigo=

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)

}}

{|
|-
|[[Archivo: untitled8.jpg|thumb|750px|left|Surface with Runge Kutta method]]

|}

{|
|-
|[[Archivo: Grafica_Runge_Kutta.jpeg|thumb|750px|left|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.

=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.

{{matlab|codigo=

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')

}}

{|
|-
|[[Archivo: untitled7.jpg|thumb|750px|left|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.

{|
|-
|[[Archivo: untitled5.jpg|thumb|750px|left|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.

{|
|-
|[[Archivo: untitled6.jpg|thumb|750px|left|Differences with the stationary solution <math>t=0,1,2,10</math> order]]

|}

= 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

<center><math>-ku_x(3,t)=0 \rightarrow </math> $\boxed{u_x(3,t)=0}$</center>

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

<center><math>
\begin{cases}
u_t-u_{xx}=0, \qquad x\epsilon(0,3), t>0\\
u(0,t)=0, u_x(3,t)=0, \qquad t>0\\
u(x,0)=u_0(x), \qquad x\epsilon[0,3]
\end{cases}
</math></center>

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)

<center><math>
\begin{cases}
u_{xx}=0, \qquad x\epsilon(0,3), t>0\\
u(0,t)=0, u_x(3,t)=0, \qquad t>0
\end{cases}
</math></center>

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

<center><math>u_{xx}=0 \rightarrow u(x,t)=C_1(t)x+C_2(t)</math></center>

Applying the boundary conditions

<center><math>C_1(t)=C_2(t)=0 \rightarrow </math> $\boxed{u(x,t)=0}$</center>

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.

==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>.

{|
|-
|[[Archivo:AP6a.png|thumb|750px|left|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.

{{matlab|codigo=
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);
}}

==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

<center><math>
\begin{cases}
\varphi’’(x)+\lambda\varphi(x)=0\\
\varphi(0)=0, \varphi’(3)=0
\end{cases}
</math></center>

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
<center><math>\varphi(0)=0 \rightarrow a=0</math>;</center>
<center><math>\varphi’(3)=0 \rightarrow \sqrt{\lambda}\,b\cos(\sqrt{\lambda}\,3)=0 \rightarrow \sqrt{\lambda}\,3=(k-{1\over2})\,\pi \rightarrow </math> $\boxed{\lambda=\mu_k=(k-{1\over2})^2{\pi^2\over9}}$;</center>

<center>Así $\boxed{\varphi_k(x)=\sin(k-{1\over2}){\pi\over3}x}$, <math> \qquad k=1,2,3…N</math></center>

If we apply that <math> u_k (x, t) = \ varphi_k (x) T_k (t) </ math> satisfies the differential equation<math> u_t-u {xx} = 0 </ math>, we obtain the differential equation determines <math>T_k(t)</math><center><math>u_t-u{xx}=0 \rightarrow \varphi_k(x)T_k’(t)-\varphi_k’’(x)T_k(t)=\varphi_k(x)T_k’(t)-(-\mu_k)\varphi_k(x)T_k(t)=0 \rightarrow </math> $\boxed{T_k’(t)+\mu_kT_k(t)=0}$</center>

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

<center>$\boxed{C_k={\int_{0}^{3}u(x,0)\varphi_k(x)dx\over\int_{0}^{3}\varphi_k^2(x)dx}}$</center>

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.
{|
|-
|[[Archivo:Untitled4.jpg|thumb|750px|left|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.

{|
|-
| [[Archivo:Untitled3.jpg|thumb|750px|left|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>

{{matlab|codigo=
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
}}

=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>0 \\
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>0 \\
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>:

{{matlab|codigo=
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

}}

{|
|-
|[[Archivo:Untitled002.jpg|thumb|750px|left|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> 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>.

{|
|-
|[[Archivo:Untitled01.jpg|thumb|750px|left|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.

== 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>:

{{matlab|codigo=
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')

}}

{|
|-
|[[Archivo:untitled.jpg|thumb|750px|left|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.

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Heat equation

2014-05-18T17:11:07Z

Marino Rivera: Marino Rivera movió la página Heat equation a Heat equation (Grupo 1B)

#REDIRECCIÓN [[Heat equation (Grupo 1B)]]

Heat equation (Grupo 1B)

2014-05-18T17:02:06Z

Marino Rivera: Página creada con «{{ TrabajoED |Heat equation. Grupo 1-B | Ecuaciones Diferenciales|Curso 2013-14 | Sandro Andrés Martínez ...»

{{ TrabajoED |Heat equation. Grupo 1-B | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandro Andrés Martínez

David Ayala Díez

Claudia Cózar Coarasa

Lorena de la Fuente Sanz

Marino Rivera Muñoz

José Manuel Torres Serrano }}

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

= Well proposed problem =

[[Archivo:Nueva imagen.png|thumb|300px|left| 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>0\\
u(0,t)=0, u(3,t)=10, \qquad t>0\\
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>0\\
u(0,t)=0, u(3,t)=10, \qquad t>0\\
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.

==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>.

{{matlab|codigo=
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.

{|
|-
| [[Archivo: untitled11.jpg|thumb|750px|left|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.

{|
|-
| [[Archivo: gráfica3a.jpg|thumb|750px|left| u(t) ]]
|}

= 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.

==Implicit Euler==

{{matlab|codigo=

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)

}}

{|
|-
|[[Archivo: untitled00010.jpg|thumb|750px|left|Surface with the implicit Euler method]]

|}

{|
|-
|[[Archivo: Grafica_Euler_impl%C3%ADcito.jpg |thumb|750px|left|<math> x = 1.5 </math> Graph with the implicit Euler method]]

|}

==Explicit Euler==

{{matlab|codigo=

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)

}}

{|
|-
|[[Archivo: untitled9.jpg|thumb|750px|left|Surface with explicit Euler method ]]

|}

{|
|-
|[[Archivo: Euler1.3.jpg |thumb|750px|left|Graph in <math>x=1.5</math> with explicit Euler method]]

|}

==Runge-Kutta==

{{matlab|codigo=

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)

}}

{|
|-
|[[Archivo: untitled8.jpg|thumb|750px|left|Surface with Runge Kutta method]]

|}

{|
|-
|[[Archivo: Grafica_Runge_Kutta.jpeg|thumb|750px|left|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.

=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.

{{matlab|codigo=

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')

}}

{|
|-
|[[Archivo: untitled7.jpg|thumb|750px|left|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.

{|
|-
|[[Archivo: untitled5.jpg|thumb|750px|left|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.

{|
|-
|[[Archivo: untitled6.jpg|thumb|750px|left|Differences with the stationary solution <math>t=0,1,2,10</math> order]]

|}

= 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

<center><math>-ku_x(3,t)=0 \rightarrow </math> $\boxed{u_x(3,t)=0}$</center>

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

<center><math>
\begin{cases}
u_t-u_{xx}=0, \qquad x\epsilon(0,3), t>0\\
u(0,t)=0, u_x(3,t)=0, \qquad t>0\\
u(x,0)=u_0(x), \qquad x\epsilon[0,3]
\end{cases}
</math></center>

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)

<center><math>
\begin{cases}
u_{xx}=0, \qquad x\epsilon(0,3), t>0\\
u(0,t)=0, u_x(3,t)=0, \qquad t>0
\end{cases}
</math></center>

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

<center><math>u_{xx}=0 \rightarrow u(x,t)=C_1(t)x+C_2(t)</math></center>

Applying the boundary conditions

<center><math>C_1(t)=C_2(t)=0 \rightarrow </math> $\boxed{u(x,t)=0}$</center>

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.

==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>.

{|
|-
|[[Archivo:AP6a.png|thumb|750px|left|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.

{{matlab|codigo=
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);
}}

==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

<center><math>
\begin{cases}
\varphi’’(x)+\lambda\varphi(x)=0\\
\varphi(0)=0, \varphi’(3)=0
\end{cases}
</math></center>

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
<center><math>\varphi(0)=0 \rightarrow a=0</math>;</center>
<center><math>\varphi’(3)=0 \rightarrow \sqrt{\lambda}\,b\cos(\sqrt{\lambda}\,3)=0 \rightarrow \sqrt{\lambda}\,3=(k-{1\over2})\,\pi \rightarrow </math> $\boxed{\lambda=\mu_k=(k-{1\over2})^2{\pi^2\over9}}$;</center>

<center>Así $\boxed{\varphi_k(x)=\sin(k-{1\over2}){\pi\over3}x}$, <math> \qquad k=1,2,3…N</math></center>

If we apply that <math> u_k (x, t) = \ varphi_k (x) T_k (t) </ math> satisfies the differential equation<math> u_t-u {xx} = 0 </ math>, we obtain the differential equation determines <math>T_k(t)</math><center><math>u_t-u{xx}=0 \rightarrow \varphi_k(x)T_k’(t)-\varphi_k’’(x)T_k(t)=\varphi_k(x)T_k’(t)-(-\mu_k)\varphi_k(x)T_k(t)=0 \rightarrow </math> $\boxed{T_k’(t)+\mu_kT_k(t)=0}$</center>

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

<center>$\boxed{C_k={\int_{0}^{3}u(x,0)\varphi_k(x)dx\over\int_{0}^{3}\varphi_k^2(x)dx}}$</center>

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.
{|
|-
|[[Archivo:Untitled4.jpg|thumb|750px|left|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.

{|
|-
| [[Archivo:Untitled3.jpg|thumb|750px|left|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>

{{matlab|codigo=
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
}}

=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>0 \\
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>0 \\
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>:

{{matlab|codigo=
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

}}

{|
|-
|[[Archivo:Untitled002.jpg|thumb|750px|left|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> 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>.

{|
|-
|[[Archivo:Untitled01.jpg|thumb|750px|left|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.

== 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>:

{{matlab|codigo=
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')

}}

{|
|-
|[[Archivo:untitled.jpg|thumb|750px|left|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.

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Archivo:AP7a.png

2014-05-15T10:11:43Z

Marino Rivera:

Archivo:AP7b.png

2014-05-15T10:02:11Z

Marino Rivera:

Archivo:AP6a.png

2014-05-15T10:01:39Z

Marino Rivera:

Archivo:AP7A.png

2014-05-15T09:53:04Z

Marino Rivera:

Archivo:AP6A.png

2014-05-15T09:52:49Z

Marino Rivera:

Archivo:AP7B.png

2014-05-15T09:51:05Z

Marino Rivera:

Archivo:AP6A.jpg

2014-05-15T09:48:52Z

Marino Rivera:

Boundary layer in laminar fluids (Grupo 1B)

2014-05-01T18:06:07Z

Marino Rivera: Deshecha la revisión 10751 de Marino Rivera (disc.)

Boundary layer in laminar fluids (Grupo 1B)

2014-05-01T18:03:31Z

Marino Rivera:

{{ TrabajoED |Boundary layer in laminar fluids. Grupo 1-B | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandro Andrés Martínez

David Ayala Díez

Claudia Cózar Coarasa

Lorena de la Fuente Sanz

Marino Rivera Muñoz

José Manuel Torres Serrano }}

In this numerical project we have studied what happens when we introduce a flat plate in a laminar fluid whose speed and viscosity are constant.

= Condición de frontera tipo Neumann =

Ahora vamos a considerar una condición de frontera diferente en el extremo derecho. En vez de suponer una temperatura constante en ese extremo como anteriormente, vamos a colocar en él una pieza aislante. Este aislante provoca que no haya pérdida de calor por el extremo derecho, es decir, que el flujo de temperatura sea nulo. Esta condición es de tipo Neumann, a diferencia de las anteriores que eran de tipo Dirichlet. Así, mantenemos la condición en el extremo izquierdo y aplicamos la nueva en el extremo derecho, que resulta

<center><math>-ku_x(3,t)=0 \rightarrow u_x(3,t)=0</math></center>

En esta situación, la temperatura de la varilla viene dada por el siguiente problema

<center><math>
\begin{cases}
u_t-u_{xx}=0, \qquad x\epsilon(0,3), t>0\\
u(0,t)=0; u_x(3,t)=0, \qquad t>0\\
u(x,0)=u_0(x), \qquad x\epsilon[0,3]
\end{cases}
</math></center>

Por lo tanto, ahora en el estado estacionario para tiempos grandes la función <math>u(x,t)</math> que modeliza la temperatura de la varilla es solución del siguiente problema de contorno (lo llamamos así pues la ecuación diferencial depende únicamente de <math>x</math> en el estado estacionario)

<center><math>
\begin{cases}
u_{xx}=0, \qquad x\epsilon(0,3), t>0\\
u(0,t)=0, u_x(3,t)=0, \qquad t>0
\end{cases}
</math></center>

Resolvemos la ecuación diferencial obteniendo el mismo resultado que anteriormente

<center><math>u_{xx}=0 \rightarrow u(x,t)=C_1(t)x+C_2(t)</math></center>

Aplicando las condiciones de frontera: <math>C_1(t)=C_2(t)=0 \rightarrow u(x,t)=0</math>

Este resultado refleja que al cabo de un tiempo lo suficientemente grande como para considerar un estado estacionario en la varilla, esta adquiere una temperatura uniforme nula. El comportamiento de la varilla es coherente con las condiciones de frontera, pues la temperatura final de la varilla es la que se mantiene constante en el extremo izquierdo, que permite flujo de temperatura.

==Método de diferencias finitas==

En la siguiente imagen se muestra una aproximación del problema mediante el método de diferencias finitas. En ella se puede apreciar como para un valor elevado del tiempo la temperatura en la varilla se puede considerar constante y uniforme en toda ella, alcanzando el valor estacionario <math>u(x,t)=0</math>.

(IMAGEN)

En concreto, a partir de un tiempo <math>t=…</math> podemos considerar que la temperatura alcanza dicho valor estacionario con un error del 5%, es decir, en ese instante la diferencia entre la distribución térmica calculada y la estacionaria es… (ESTA PARTE HAY QUE COMPLETARLA CUANDO PREGUNTEMOS A CARLOS A QUÉ SE REFIERE EL 5%).

A continuación se refleja el código de Matlab que aproxima la temperatura de la varilla por diferencias finitas usando el método de Euler implícito por ser el que aporta una mejor aproximación.

==Método de Fourier==

Planteamos ahora el mismo problema usando el método de Fourier. Por tanto, buscamos soluciones de la forma <math>u(x,t)=\varphi(x)T(t)</math>, donde <math>\varphi(x)</math> satisface el siguiente problema de autovalores:

<math>
\begin{cases}
\varphi’’(x)+\lambda\varphi(x)=0 ;\\
\varphi(0)=0; \varphi’(3)=0
\end{cases}
</math>

La solución de la ecuación diferencial es <math>\varphi(x)=a\cos(\sqrt{\lambda}x)+b\sin(\sqrt{\lambda}x)</math>

Aplicando las condiciones de frontera obtenemos los autovalores <math>\mu_k</math> y autofunciones <math>\varphi_k(x)</math> del problema
<center><math>\varphi(0)=0 \rightarrow a=0</math>;</center>
<center><math>\varphi’(3)=0 \rightarrow \sqrt{\lambda}b\cos(\sqrt{\lambda}\,3)=0 \rightarrow \sqrt{\lambda}\,3=(k-{1\over2})\pi \rightarrow \lambda=\mu_k=(k-{1\over2})^2{\pi^2\over9}</math>;</center>
<center>

Así <math>\boldsymbol{\varphi_k(x)=\sin(k-{1\over2}){\pi\over3}x} \qquad k=1,2,3…N</math></center>

Si aplicamos que <math>u_k(x,t)=\varphi_k(x)T_k(t)</math> satisface la ecuación diferencial <math>u_t-u{xx}=0</math>, obtenemos la ecuación diferencial que determina <math>T_k(t)</math>
<center><math>u_t-u{xx}=0 \rightarrow \varphi_k(x)T_k’(t)-\varphi_k’’(x)T_k(t)=\varphi_k(x)T_k’(t)-(-\mu_k)\varphi_k(x)T_k(t)=0 \rightarrow T_k’(t)+\mu_kT_k(t)=0</math></center>

La solución de esta ecuación diferencial es <math>T_k(t)=C_ke^{-\mu_kt}= C_ke^{-(k-{1\over2})^2{\pi^2\over9}t}</math> y por tanto <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>

Si expresamos la solución del problema como <math>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</math> y hacemos que satisfaga la condición inicial obtenemos <math>u(x,0)=\sum_{k=1}^NC_k\sin(k-{1\over2}){\pi\over3}x</math>. De esta forma, por unicidad de los coeficientes de Fourier, los coeficientes <math>C_k</math> coinciden con los de la serie de Fourier con respecto a las autofunciones <math>\varphi_k(x)</math> de la función a trozos que determina la condición inicial (expresada al principio del artículo).

El problema se limita así al cálculo de dichos coeficientes, de acuerdo con la expresión

<center><math>C_k={\int_{0}^{3}u(x,0)\varphi_k(x)dx\over\int_{0}^{3}\varphi_k^2(x)dx}</math></center>

El grado de exactitud de la aproximación con este método depende del número de elementos de la serie de Fourier, es decir, del valor de <math>k</math>. Estudiamos la temperatura de la varilla tomando los valores <math>k=1,3,5,10,20</math>, como se puede ver en la siguiente imagen. En ella se aprecia, sobre todo en la condición inicial, que a medida que se aumente el valor de <math>k</math> la función aproximada se acerca más a la real.

Estos resultados se pueden comparar mejor tomando un instante fijo <math>t=0.5</math> y representando en él cada función en una misma gráfica, como podemos ver a continuación.

First, we must assume that the fluid velocity before reaching the plate is constant, as in remote areas after passing the plate.

In our case we take this constant as <math>2</math>, such that

<math>\overrightarrow{u} = u_0\cdot\overrightarrow{i}, u_0 = 2 </math>
[[Archivo:Flujo.png|thumb|200px|left|Laminar fluid with velocity <math> u_0 </math>]]

Then, we must define the fluid stream function

<math>\psi(x,y) = \sqrt[]{ \nu \cdot\ u_0 \cdot x} f(\eta)</math>

Where we take the viscosity <math> \nu </math> as a unit value and

<math>\eta = y \sqrt[]{ \frac{u_0}{\nu x}}</math>

<math> f(\eta) </math> Satisfies the Blasius equation, and therefore we will raise the initial value problem associated with this equation with the following initial conditions

<math>
\begin{cases}
f’’’(\eta)+\frac{1}{2}f(\eta)f’’(\eta)=0 ; \\
f(0)=f’(0)=0, \lim_{\eta \to \infty}f’(\eta)= 1 ;
\end{cases}
</math>

However, on time of programming we cannot introduce a conditional limit, so we have to replace them by the condition <math> f’’(0)=k </math>, and vary the values of k to find the one that satisfies the limit.

We cannot solve the differential equation like such, to apply the numerical methods, we need to pass it to a system of equations:

<math>
f(\eta)=y_1,f’(\eta)=y_2,f’’(\eta)=y_3\\
\begin{cases}
y_1’=y_2;\\
y_2’=y_3;\\
y_3’=-\frac{1}{2}y_1y_3;\\
y_1(0)=y_2(0)=0; y_3(0)=k
\end{cases}
</math>

Once the system has been formulated , we start to solve it.

== Resolution with the modified Euler method ==
Then is exposed the Matlab code that numerically that solves the Blasius equation for different values of <math>k</math> , with <math> k \in \mbox{(0,1;1)}</math> and <math>dk=0,01</math> ,with <math> \eta \in \mbox{(0,20)}</math> and with <math>h=0,05</math> .

{{matlab|codigo=
%Resolution of Blasius equation(with modified Euler)
clear all
%Initial conditions
t0=0;
tN=20;
h=0.05;
N=(tN-t0)/h;
F2=zeros(91,401);%We create the matrix F2 where we will store the different
%solutions of f2 for each value of k
for k=0.1:0.01:1
y=[0;0;k];
y1=y(1);
y2=y(2);
y3=y(3);
for n=1:N
A=[0 1 0;0 0 1;(-y(3)/2) 0 0]; %To simplify and solve using matrices, we create
%the matrix A in the loop with different values of f3
z=y+h*A*y;
y=y+(h/2)*(A*y+A*z);
y1(n+1)=y(1);
y2(n+1)=y(2);
y3(n+1)=y(3);
end
t=[t0:h:tN];
num=int8(100*(k-0.1+0.01));
%F2 has as rows approximations of y2 for the different values of k
F2(num,:)=y2;
end
k1=[0.1:0.01:1]; %Vector to represent the values of f2 in 20
f20=F2(:,401);
f20=f20';
o=ones(1,91);%We represented f = 1 for better viewing
figure(1)
hold on
plot(k1,f20,'+')
plot(k1,o,'r')
xlabel('k')
ylabel('f´(20)')
legend('f´(20)','y=1')
hold off

}}

{|
|-
|[[Archivo:Graficaf'(20)M.jpg|thumb|500px|left|Graph of <math>f’(20)</math> for each <math>k</math> ]]

|}
As noted in the graph the value for which the function is closer to <math>1</math> is <math>k=0,33</math>.

==Resolution with 4th order Runge-Kutta method==
Then is exposed the Matlab code that numerically that solves the Blasius equation for different values of <math>k</math> , with <math> k \in \mbox{(0,1;1)}</math> and <math>dk=0,01</math> ,with <math> \eta \in \mbox{(0,20)}</math> and with <math>h=0,05</math> .

{{matlab|codigo=
% Resolution of Blasius equation(with Runge-Kutta)
clear all
t0=0;
tN=20;
h=0.05;
N=(tN-t0)/h;
F2=zeros(91,401);
for k=0.1:0.01:1
y=[0;0;k];
y1=y(1);
y2=y(2);
y3=y(3);
for n=1:N
A=[0 1 0;0 0 1;(-y(3)/2) 0 0];
k1=A*y;
k2=A*(y+(h/2)*k1);
k3=A*(y+(h/2)*k2);
k4=A*(y+h*k3);
y=y+(h/6)*(k1+2*k2+2*k3+k4);
y1(n+1)=y(1);
y2(n+1)=y(2);
y3(n+1)=y(3);
end
t=[t0:h:tN];
num=int8(100*(k-0.1+0.01));
%F2 has as rows approximations of y2 for the different values of k
F2(num,:)=y2;
end
t=[t0:h:tN];
k1=[0.1:0.01:1];
f20=F2(:,401);
f20=f20';
o=ones(1,91);
figure(1)
hold on
plot(k1,f20,'+')
plot(k1,o,'r')
xlabel('k')
ylabel('f´(20)')
legend('f´(20)','y=1')
hold off

}}
{|
|-
| [[Archivo: Graficaf'(20)RK.jpg|thumb|500px|left|Graph of <math>f’(20)</math> for each <math>k</math>]]
|}

==Resolution with Euler method==
Then is exposed the Matlab code that numerically that solves the Blasius equation for different values of <math>k</math> , with <math> k \in \mbox{(0,1;1)}</math> and <math>dk=0,01</math> ,with <math> \eta \in \mbox{(0,20)}</math> and with <math>h=0,05</math> .

{{matlab|codigo=
%Resolution of Blasius equation(with Euler)

clear all
t0=0;
tN=20;
h=0.05;
N=(tN-t0)/h;
F2=zeros(91,401);
for k=0.1:0.01:1
y=[0;0;k];
y1=y(1);
y2=y(2);
y3=y(3);
for n=1:N
A=[0 1 0;0 0 1;(-y(3)/2) 0 0];
y=y+h*A*y;
y1(n+1)=y(1);
y2(n+1)=y(2);
y3(n+1)=y(3);
end
t=[t0:h:tN];
num=int8(100*(k-0.1+0.01));
%F2 has as rows approximations of y2 for the different values of k
F2(num,:)=y2;
end
t=[t0:h:tN];
k1=[0.1:0.01:1];
f20=F2(:,401);
f20=f20';
o=ones(1,91);
figure(1)
hold on
plot(k1,f20,'+')
plot(k1,o,'r')
xlabel('k')
ylabel('f´(20)')
legend('f´(20)','y=1')
hold off

}}
{|
|-
| [[Archivo: Graficaf'(20).jpg|thumb|500px|left|Graph of <math>f’(20)</math> for each <math>k</math>]]
|}
== Conclusion==
Comparing the graphs, we see that the difference between modified Euler and 4th order Runge Kutta methods is minimal and the value for each parameter is <math>k=0,33</math> , on the other hand, by using the Euler method (less accurate than the above) the value of the parameter is <math>k=0,32</math>, although graphically this difference is hardly seen.

==Graph of <math>f´(\eta)</math>==

The graph of <math>f'(\eta)</math>has been realized in <math>\eta \in \mbox {(0,20)}</math> for the value of parameter <math>k</math> obtained in the modified Euler method, ie <math>k=0.33</math>.

To get this graph we add into the modified Euler method, the following MATLAB code:

{{matlab|codigo=

%We plot f2 for k = 0.33 which is in row 24 of the matrix
f2=F2(24,:);
figure(2)
plot(t,f2,'*')
xlabel('\eta')
ylabel('f´(\eta)')
}}

{|

|-

| [[Archivo: Graficaf'eta.jpg|thumb|800px|left|Graph of <math>f´(\eta)</math> for <math>k=0,33</math> ]]

|}

As can be seen, if we run the full program and see the vector <math>f2</math> the value of <math>\eta_0</math> for which <math>\ \vert f'(\eta)-1 \vert < 0,01</math> , if <math>\eta>\eta _0</math>, is <math>\eta_0\ge5,95</math> for <math> \eta \in \mbox{(0,20)}</math>.

=Horizontal velocity of the fluid=
Once we have numerically calculated the <math>f( \eta)</math> we proceed to calculate the horizontal component of the fluid velocity <math> u_1 </math>:

<math>\vec{u}=(u_1,u_2)=(\frac{\partial \psi}{\partial y},-\frac{\partial \psi}{\partial x})</math> Thus, as defined above:

<math>u_1=\frac{\partial \psi}{\partial y}=\frac{\partial }{\partial y}(\sqrt[]{\nu u_0 x} f(\eta))= \sqrt[]{\nu u_0 x} \frac{\partial f(\eta)}{\partial y}=\sqrt[]{\nu u_0 x} \frac{\partial f(\eta)}{\partial \eta} \frac{\partial \eta}{\partial y}=\sqrt[]{\nu u_0 x} \sqrt[]{\frac{u_0}{\nu x}} f’(\eta)=u_0 f’(\eta)</math>

To translate this result graphically we calculate <math>u_1(x_k,y)</math> with the Modified Euler method where <math>x_k=0.05,0.2,0.4,0.6,0.8</math> and <math>y \in \mbox{(0,3)}</math> with <math>h=0.01</math>.

{{matlab|codigo=
%We calculate u1 with the different values of xk(with modified Euler)
clear all
xk=[0.05,0.2,0.4,0.6,0.8];
nu=1; u0=2;
y0=0; yN=3; hy=0.01;
N=(yN-y0)/hy;
y=y0:hy:yN;
for m=1:5
%We define eta ('t') for each value of xk, each one
%corresponding to a row, and with 'y' in (0,3)
t(m,:)=sqrt(u0/(nu*xk(m)))*y;
h=sqrt(u0/(nu*xk(m)))*0.01;
f0=[0;0;0.33];
f=[f0(1);f0(2);f0(3)];
for n=1:N
A=[0 1 0;0 0 1;(-f0(3)/2) 0 0];
z=f0+h*A*f0;
f0=f0+(h/2)*(A*f0+A*z);
f(:,n+1)=[f0(1);f0(2);f0(3)];
end
Y(m,:)=f(2,:);
end
F=u0*Y;
hold on
plot(y,F(1,:),'k')
plot(y,F(2,:))
plot(y,F(3,:),'r')
plot(y,F(4,:),'m')
plot(y,F(5,:),'g')
legend('u_{1} for x_{k}=0.05','u_{1} for x_{k}=0.2','u_{1} for x_{k}=0.4','u_{1} for x_{k}=0.6','u_{1} for x_{k}=0.8','location','best')
xlabel('y')
ylabel('u_{1}')
hold off
}}

[[Archivo:Ap4a.png|thumb|800px|centre|In this picture are shown <math>5</math> graphs, each one corresponds to a different value of <math>k</math>. ]]

According to this graph we can appreciate that the fluid, when is moving along the <math>x</math> axis, has to achieve higher height to get limit velocity <math>u_0</math>, i.e., when the fluid moves, it must be getting over the plate to offset the perturbation that the plate causes to it ( to a higher value of <math>y</math> is greater the transition zone between zero velocity of the plate and the limit velocity <math>u_0</math> with which the fluid initially starts).

= Laminar boundary layer of the fluid =

According to the above findings, it can be deduced that there is for each value of <math>x</math> a limit value <math>y</math> from which the fluid velocity becomes constant speed again, with the same value that it had initially before reaching the area of the plate.
Obviously the value of the boundary layer will be related to the value <math>\eta_0</math>, calculated above, for which the function <math>f’(\eta)</math> will become almost constant. The relationship is expressed as follows:
<math>\eta = y \sqrt[]{ \frac{u_0}{\nu x}}</math>

If <math>\eta=\eta_0</math>, then

<center><math>\eta_0=y \sqrt[]{ \frac{2}{x}}; y=\frac{\eta_0\sqrt[]{x}}{\sqrt[]{2}}=g(x);</math></center>

Therefore we can interpret this function <math>g(x)</math> as the fluid boundary layer. Simple Matlab code is exposed for its representation:

{{matlab|codigo=
%plot function g(x)
x=[0:0.05:10];
eta0=5.95;
y=eta0*(x/2).^(1/2);
plot(x,y,'r')
xlabel('x')
ylabel('g(x)')
legend('g(x), interpreted as boundary layer','location','best')

}}

{|
|-
| [[Archivo: Figureg(x).jpg|thumb|500px|left|Graph of <math>g(x)</math> ]]
|}

In view of the graph, the findings are similar to previous ones. As the value <math>x</math> increases, you need a higher value of <math>y</math> for the fluid velocity stabilizes and becomes the speed that we had initially.

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Boundary layer in laminar fluids (Grupo 1B)

2014-03-07T09:05:38Z

Marino Rivera:

Boundary layer in laminar fluids (Grupo 1B)

2014-03-07T09:01:32Z

Marino Rivera:

Boundary layer in laminar fluids (Grupo 1B)

2014-03-06T23:06:48Z

Marino Rivera:

Boundary layer in laminar fluids (Grupo 1B)

2014-03-06T22:49:24Z

Marino Rivera:

Archivo:Ap4b.png

2014-02-26T23:01:12Z

Marino Rivera:

Archivo:Ap4a.png

2014-02-26T23:00:19Z

Marino Rivera: