Ecuaciones diferenciales en
derivadas parciales
Pedro Corcuera
Dpto. Matemática Aplicada y
Ciencias de la Computación
Universidad de Cantabria
[email protected]
Objetivos
• Revisión de las ecuaciones diferenciales en
derivadas parciales y técnicas de solución
Ec. derivadas parciales
2
Diferenciación
Ec. derivadas parciales
3
Diferenciación de funciones continuas
• The mathematical definition of a derivative begins
with a difference approximation: ∆y = f (xi + ∆x) − f (xi )
∆x
∆x
and as ∆x is allowed to approach zero, the difference
becomes a derivative: dy = lim f (xi + ∆x) − f (xi )
dx
∆x→0
Ec. derivadas parciales
∆x
4
Fórmulas de diferenciación de gran
exactitud
• Taylor series expansion can be used to generate
high-accuracy formulas for derivatives by using linear
algebra to combine the expansion around several
points.
• Three categories for the formula include forward
finite-difference, backward finite-difference, and
centered finite-difference.
Ec. derivadas parciales
5
Aproximación por diferencia en adelanto
lim f ( x + Δx ) − f ( x )
f ′( x ) =
Δx → 0
Δx
For a finite ' Δx'
f ′( x ) ≈
f (x + ∆x ) − f ( x )
∆x
f(x)
x
x+Δx
Graphical Representation of forward difference approximation of first derivative.
Ec. derivadas parciales
6
Aproximación por diferencia en adelanto ejemplo
• The velocity of a rocket is given by


14 ×10 4
− 9.8t ,
ν (t ) = 2000 ln 

4
14 ×10 − 2100t 
0 ≤ t ≤ 30
• Where 'ν' is given in m/s and 't ' is given in seconds.
a)Use forward difference approximation of the first derivative
of ν(t ) to calculate the acceleration at t = 16s. Use a step
size of Δt = 2s.
b)Find the exact value of the acceleration of the rocket.
c) Calculate the absolute relative true error for part (b).
Ec. derivadas parciales
7
Aproximación por diferencia en adelanto ejemplo
• Solution
a)
a(ti ) ≈
ν (ti +1 ) −ν (ti )
a(16 ) ≈
Δt = 2
ti = 16
ti +1 = ti + ∆t = 16 + 2 = 18
∆t
ν (18) −ν (16)
2


14 × 10 4
ν (18) = 2000 ln 
− 9.8(18) = 453.02 m/s

4
14 × 10 − 2100(18)
Hence


14 ×10 4
ν (16 ) = 2000 ln 
 − 9.8(16 ) = 392.07 m/s
4
14 ×10 − 2100(16 ) 
ν (18) −ν (16)
453.02 − 392.07
≈ 30.474m/s 2
2
2
b) The exact value of a (16 ) can be calculated by differentiating


14 ×10 4
d
9
.
8
t
ν (t ) = 2000 ln 
−
as
a(t ) = [ν(t )]

4
14 ×10 − 2100t 
dt
a(16 ) ≈
≈
Ec. derivadas parciales
8
Aproximación por diferencia en adelanto ejemplo
Knowing that
d
[ln(t )] = 1 and
dt
t
1
d 1
=
−
dt  t 
t2
 14 × 10 4 − 2100t  d 

14 × 10 4


 − 9.8
a (t ) = 2000
4
4


14 × 10

 dt  14 × 10 − 2100t 
4


 14 × 10 4 − 2100t 
14
10
×
(− 2100 ) − 9.8

(− 1)
= 2000
2 
4
4

14 × 10


 (14 × 10 − 2100t ) 
− 4040 − 29.4t
− 200 + 3t
− 4040 − 29.4(16 )
a (16 ) =
= 29.674m/s 2
− 200 + 3(16)
=
c) The absolute relative true error is
29.674 − 30.474
True Value - Approximate Value
=
x100 = 2.6967 %
∈t =
⋅100
29.674
True Value
Ec. derivadas parciales
9
Efecto del tamaño de paso en el método
de diferencia dividida en adelanto
f ( x ) = 9e 4 x
Value of´ f ′(0.2) using forward difference method.
h
0.05
0.025
0.0125
0.00625
0.003125
0.001563
0.000781
0.000391
0.000195
9.77E-05
4.88E-05
f ' (0.2)
88.69336
84.26239
82.15626
81.12937
80.62231
80.37037
80.24479
80.18210
80.15078
80.13512
80.12730
εa %
Ea
-4.430976
-2.106121
-1.0269
-0.507052
-0.251944
-0.125579
-0.062691
-0.031321
-0.015654
-0.007826
5.258546
2.563555
1.265756
0.628923
0.313479
0.156494
0.078186
0.039078
0.019535
0.009767
Significant
digits
0
1
1
1
2
2
2
3
3
3
Ec. derivadas parciales
Et
εt %
-8.57389
-4.14291
-2.03679
-1.00989
-0.50284
-0.25090
-0.12532
-0.06263
-0.03130
-0.01565
-0.00782
10.70138
5.170918
2.542193
1.260482
0.627612
0.313152
0.156413
0.078166
0.039073
0.019534
0.009766
10
Efecto del tamaño de paso en el método
de diferencia dividida en adelanto
92
Number of times step size halved, n
f'(0.2)
88
0
0
84
2
4
6
8
10
12
-1
-2
Ea
80
-3
76
-4
0
1
2
3
4
5
6
7
8
9
10
11
12
-5
Num ber of tim es step size halved, n
4
Least number of significant digits
correct
6
5
|Ea| %
4
3
2
1
0
0
1
2
3
4
5
6
7
8
9
10
11
12
3
2
1
0
1
Num ber of tim es step size halved, n
Ec. derivadas parciales
2
3
4
5
6
7
8
9
10
11
Num ber of tim es step size halved, n
11
Efecto del tamaño de paso en el método
de diferencia dividida en adelanto
12
Num ber of tim es step size halved, n
10
0
Et
-3
2
4
6
8
10
12
8
|Et| %
0
6
4
-6
2
0
-9
1
2
3
4
5
6
7
8
9
10
11
Number of times step size halved, n
Ec. derivadas parciales
12
Aproximación por diferencia en atraso
• Sabemos que
f ′( x ) =
lim f ( x + Δx ) − f ( x )
Δx → 0
Δx
• For a finite ' Δx' ,
f ′( x ) ≈
f (x + ∆x ) − f ( x )
∆x
• If ' Δx' is chosen as a negative number,
f ′( x ) ≈
f ( x − ∆x ) − f ( x ) f (x ) − f (x − Δx )
=
Δx
− ∆x
Ec. derivadas parciales
13
Aproximación por diferencia en atraso
• This is a backward difference approximation as you
are taking a point backward from x. To find the value
of f ′(x ) at x = xi, we may choose another point ' Δx'
behind as x = xi −1 . This gives
f ′( xi ) ≈
f ( xi ) − f ( xi −1 )
∆x
f ( xi ) − f ( xi −1 )
=
xi − xi −1
where
Δx = xi − xi −1
f(x)
x
x-Δx
x
Graphical Representation of backward
difference approximation of first derivative.
Ec. derivadas parciales
14
Obtención de la adad a partir de las series
de Taylor
• Taylor’s theorem says that if you know the value of a
function f at a point xi and all its derivatives at that
point, provided the derivatives are continuous
between xi and xi +1 , then
f ( xi +1 ) = f ( xi ) + f ′( xi )( xi +1 − xi ) +
Substituting for convenience
f ′′( xi )
(xi +1 − xi )2 + 
2!
Δx = xi +1 − xi
f ′′( xi )
(Δx )2 + 
f (xi +1 ) = f (xi ) + f ′(xi )Δx +
2!
f ( xi +1 ) − f ( xi ) f ′′(xi )
′
(∆x ) + 
f ( xi ) =
−
∆x
2!
f ( xi +1 ) − f ( xi )
f ′( xi ) =
+ O(∆x )
∆x
Ec. derivadas parciales
15
Obtención de la adad a partir de las series
de Taylor
• The O(∆x ) term shows that the error in the approximation is of the order of ∆x . It is easy to derive from
Taylor series the formula for backward divided
difference approximation of the first derivative.
• As shown above, both forward and backward divided
difference approximation of the first derivative are
accurate on the order of O(∆x ) .
• Can we get better approximations? Yes, another
method is called the Central difference approximation of the first derivative.
Ec. derivadas parciales
16
Obtención de la adc a partir de las series
de Taylor
• From Taylor series
f ′′( xi )
f ′′′(xi )
2
(Δx ) +
(Δx )3 + 
2!
3!
′′′( )
f ′′( xi )
(Δx )2 − f xi (Δx )3 + 
f ( xi −1 ) = f ( xi ) − f ′( xi )Δx +
2!
3!
f ( xi +1 ) = f ( xi ) + f ′( xi )Δx +
(1)
(2)
Subtracting equation (2) from equation (1)
f ( xi +1 ) − f ( xi −1 ) = f ′( xi )(2∆x ) +
2 f ′′′( xi )
(∆x )3 + 
3!
f ′( xi ) =
f ( xi +1 ) − f ( xi −1 ) f ′′′( xi )
(∆x )2 + 
−
2∆x
3!
f ′( xi ) =
f (xi +1 ) − f (xi −1 )
2
+ O(∆x )
2∆x
Ec. derivadas parciales
17
Obtención de la adc a partir de las series
de Taylor
• Hence showing that we have obtained a more
2
accurate formula as the error is of the order of O(∆x )
f(x)
x
x-Δx
x
x+Δx
Graphical Representation of central difference approximation of first derivative
Ec. derivadas parciales
18
Aproximación por diferencia central ejemplo
• The velocity of a rocket is given by


14 ×10 4
− 9.8t ,
ν (t ) = 2000 ln 

4
14 ×10 − 2100t 
0 ≤ t ≤ 30
• Where 'ν' is given in m/s and 't ' is given in seconds.
a)Use central divided difference approximation of the first
derivative of ν(t ) to calculate the acceleration at t = 16s .
Use a step size of Δt = 2s.
b)Calculate the absolute relative true error for part (a).
Ec. derivadas parciales
19
Aproximación por diferencia central ejemplo
• Solution
a)
a(ti ) ≈
ν (ti +1 ) −ν (ti −1 )
a(16 ) ≈
Δt = 2
ti = 16
ti −1 = ti − ∆t = 16 − 2 = 14
ti +1 = ti + ∆t = 16 + 2 = 18
2∆t
ν (18) −ν (14)
2⋅2


14 ×10 4
ν (18) = 2000 ln 
− 9.8(18) = 453.02 m/s

4
14 ×10 − 2100(18) 
Hence


14 ×10 4
ν (14) = 2000 ln 
− 9.8(14 ) = 334.24 m/s

4
14 ×10 − 2100 (14 ) 
ν (18) −ν (14)
453.02 − 334.24 ≈ 29.694 m/s 2
4
4
b) The absolute relative true error is knowing that the exact value at t = 16 s
is a(16 ) = 29.674 m/s 2
29.674 − 29.694
∈t =
⋅100 = 0.069157 %
29.674
a(16 ) ≈
≈
Ec. derivadas parciales
20
Efecto del tamaño de paso en el método
de diferencia dividida central
f ( x ) = 9e 4 x
Value of´ f ′(0.2) using forward difference method.
h
0.05
0.025
0.0125
0.00625
0.003125
0.001563
0.000781
0.000391
0.000195
9.77E-05
4.88E-05
f ' (0.2)
80.65467
80.25307
80.15286
80.12782
80.12156
80.12000
80.11960
80.11951
80.11948
80.11948
80.11947
εa %
Ea
-0.4016
-0.100212
-0.025041
-0.00626
-0.001565
-0.000391
-9.78E-05
-2.45E-05
-6.11E-06
-1.53E-06
0.500417
0.125026
0.031252
0.007813
0.001953
0.000488
0.000122
3.05E-05
7.63E-06
1.91E-06
Significant
digits
1
2
3
3
4
5
5
6
6
7
Ec. derivadas parciales
Et
εt %
-0.53520
-0.13360
-0.03339
-0.00835
-0.00209
-0.00052
-0.00013
-0.00003
-0.00001
0.00000
0.00000
0.668001
0.16675
0.041672
0.010417
0.002604
0.000651
0.000163
4.07E-05
1.02E-05
2.54E-06
6.36E-07
21
Efecto del tamaño de paso en el método
de diferencia dividida central
95
Num ber of steps involved, n
3
1
7
5
9
11
85
E(a)
f'(0.2)
0
75
1
3
5
7
9
11
-1
Least number of significant digits
correct
Num ber of tim es the step size is halved, n
0.6
0.5
|E(a)|,%
0.4
0.3
0.2
0.1
0
0
1
2
3
4
5
6
7
8
9
10
11
12
8
7
6
5
4
3
2
1
0
1
2
Num ber of steps involved, n
3
4
5
6
7
8
9
10
11
Num ber of steps involved, n
Ec. derivadas parciales
22
Efecto del tamaño de paso en el método
de diferencia dividida central
0.8
0.7
Num ber of steps involved, n
0.6
-0.1
E(t)
-0.2
-0.3
-0.4
-0.5
-0.6
0
2
4
6
8
10
12
|Et| %
0.0
0.5
0.4
0.3
0.2
0.1
0
1
2
3
4
5
6
7
8
9
10
11
Num ber of tim es step size halved,n
Ec. derivadas parciales
23
Comparación de métodos adad, adat, adc
• The results from the three difference approximations
are given in the following table.
Summary of a (16) using different divided
difference approximations
Type
of
Difference
Approximation
Forward
Backward
Central
a(16 )
(m / s )
2
30.475
28.915
29.695
Ec. derivadas parciales
∈t %
2.6967
2.5584
0.069157
24
Obtención del valor de la derivada con
una tolerancia específica
• In real life, one would not know the exact value of the
derivative – so how would one know how accurately they
have found the value of the derivative.
• A simple way would be to start with a step size and keep on
halving the step size until the absolute relative approximate
error is within a pre-specified tolerance.
• Take the example of finding v ′(t ) for
at


14 × 10 4
− 9.8t
ν (t ) = 2000 ln 

4
14 × 10 − 2100t 
t = 16 using the backward divided difference
Ec. derivadas parciales
scheme.
25
Obtención del valor de la derivada con
una tolerancia específica
• The values obtained using the backward difference
approximation method and the corresponding
absolute relative approximate errors are given in the
following table.
∆t
2
1
0.5
0.25
0.125
v ′(t )
28.915
29.289
29.480
29.577
29.625
∈a %
1.2792
0.64787
0.32604
0.16355
one can see that the absolute
relative approximate error
decreases as the step size is
reduced. At ∆t = 0.125 the absolute
relative approximate error is
0.16355%, meaning that at least 2
significant digits are correct in the
answer.
Ec. derivadas parciales
26
Fórmulas de diferencia finita en adelanto
Ec. derivadas parciales
27
Fórmulas de diferencia finita en atraso
Ec. derivadas parciales
28
Fórmulas de diferencia finita centrada
Ec. derivadas parciales
29
Sistemas de ecuaciones
lineales
Ec. derivadas parciales
30
Matrices
• A matrix consists of a rectangular array of elements
represented by a single symbol (example: [A]).
• An individual entry of a matrix is an element
(example: a23)
Ec. derivadas parciales
31
Matrices especiales
• Matrices where m=n are called square matrices.
• There are a number of special forms of square
matrices:
Symmetric
5 1 2


[A] = 1 3 7
2 7 8
Upper Triangular
a11

[A] = 

a12
a22
a13 

a23 
a33 
Diagonal
a11

a22
[A] = 




a33 
Identity
1



[A] =  1 

1
Lower Triangular
Banded
a11

a11



[A] = a21 a22
 [A] = a21
a31 a32 a33 



Ec. derivadas parciales
a12
a22
a32
a23
a33
a43



a34 

a44 
32
Multiplicación matricial
• The elements in the matrix [C] that results from
multiplying matrices [A] and [B] are calculated using:
n
c ij = ∑ aik bkj
k=1
Ec. derivadas parciales
33
Representación de Algebra Lineal
• Matrices provide a concise notation for representing
and solving simultaneous linear equations:
a11 x1 + a12 x 2 + a13 x 3 = b1
a21 x1 + a22 x 2 + a23 x 3 = b2
a31 x1 + a32 x 2 + a33 x 3 = b3
a11

a21
a31
a12
a22
a32
a13  x1   b1 
   
a23  x 2  = b2 
a33  x 3  b3 
[A]{x} = {b}
Ec. derivadas parciales
34
Solución con Matlab
• MATLAB provides two direct ways to solve systems of
linear algebraic equations [A]{x}={b}:
– Left-division
x = A\b
– Matrix inversion
x = inv(A)*b
• The matrix inverse is less efficient than left-division
and also only works for square, non-singular systems.
Ec. derivadas parciales
35
Matriz inversa
• Recall that if a matrix [A] is square, there is another
matrix [A]-1, called the inverse of [A], for which
[A][A]-1=[A]-1[A]=[I]
• The inverse can be computed in a column by column
fashion by generating solutions with unit vectors as
the right-hand-side constants:
1
[A]{x1} = 0
0
0
[A]{x 2 } = 1
0
[A]
−1
= [x1
x2
Ec. derivadas parciales
0
[A]{x 3} = 0
1
x3 ]
36
Matriz inversa y sistemas estímulo respuesta
• Recall that LU factorization can be used to efficiently
evaluate a system for multiple right-hand-side vectors
- thus, it is ideal for evaluating the multiple unit vectors
needed to compute the inverse.
• Many systems can be modeled as a linear
combination of equations, and thus written as a matrix
equation:
[Interactions]{response} = {stimuli}
• The system response can thus be found using the
matrix inverse.
Ec. derivadas parciales
37
Normas vectoriales y matriciales
• A norm is a real-valued function that provides a
measure of the size or “length” of multi-component
mathematical entities such as vectors and matrices.
• Vector norms and matrix norms may be computed
differently.
Ec. derivadas parciales
38
Normas vectoriales
• For a vector {X} of size n, the p-norm is:
n

p
= ∑ x i 

 i=1
1/ p
X
p
• Important examples of vector p-norms include:
p = 1 : sum of the absolute values
n
X 1 = ∑ xi
i =1
p = 2 : Euclidian norm (length)
p = ∞ : maximum − magnitude
X
X
Ec. derivadas parciales
2
= Xe=
∞
= max xi
n
2
x
∑ i
i =1
1≤i ≤ n
39
Normas matriciales
• Common matrix norms for a matrix [A] include:
n
column - sum norm
A 1 = max ∑ aij
1≤ j≤n
Frobenius norm
Af =
i=1
n n
2
a
∑∑ ij
i=1 j=1
n
row - sum norm
A ∞ = max ∑ aij
1≤i≤n
spectral norm (2 norm)
j=1
1/2
A 2 = (µ max )
• Note - µmax is the largest eigenvalue of [A]T[A].
Ec. derivadas parciales
40
Número de condición de una matriz
• The matrix condition number Cond[A] is obtained by
calculating Cond[A]=||A||·||A-1||
• In can be shown that: ∆X ≤ Cond[A] ∆A
X
A
• The relative error of the norm of the computed
solution can be as large as the relative error of the
norm of the coefficients of [A] multiplied by the
condition number.
• If the coefficients of [A] are known to t digit precision,
the solution [X] may be valid to only
t-log10(Cond[A]) digits.
Ec. derivadas parciales
41
Comandos Matlab
• MATLAB has built-in functions to compute both norms
and condition numbers:
– norm(X,p)
• Compute the p norm of vector X, where p can be any
number, inf, or ‘fro’ (for the Euclidean norm)
– norm(A,p)
• Compute a norm of matrix A, where p can be 1, 2, inf, or
‘fro’ (for the Frobenius norm)
– cond(X,p) or cond(A,p)
• Calculate the condition number of vector X or matrix A using
the norm specified by p.
Ec. derivadas parciales
42
Método iterativo: Gauss - Seidel
• The Gauss-Seidel method is the most commonly used iterative
method for solving linear algebraic equations [A]{x}={b}.
• The method solves each equation in a system for a particular
variable, and then uses that value in later equations to solve later
variables. For a 3x3 system with nonzero elements along the
diagonal, for example, the jth iteration values are found from the j1th iteration using:
j−1
j−1
b1 − a12 x2 − a13 x3
x =
a11
j
1
b2 − a21 x1j − a23 x3j−1
x =
a22
j
2
b3 − a31 x1j − a32 x2j
x =
a33
j
3
Ec. derivadas parciales
43
Iteración de Jacobi
•
The Jacobi iteration is similar to the Gauss-Seidel
method, except the j-1th information is used to
update all variables in the jth iteration:
a) Gauss-Seidel
b) Jacobi
Ec. derivadas parciales
44
Convergencia
• The convergence of an iterative method can be
calculated by determining the relative percent change
of each element in {x}. For example, for the ith
element in the jth iteration,
ε a ,i
xij − xij −1
=
×100%
j
xi
• The method is ended when all elements have
converged to a set tolerance.
Ec. derivadas parciales
45
Dominancia diagonal
• The Gauss-Seidel method may diverge, but if the
system is diagonally dominant, it will definitely
converge.
• Diagonal dominance means:
n
aii > ∑ aij
j=1
j≠i
Ec. derivadas parciales
46
Relajación
• To enhance convergence, an iterative program can
introduce relaxation where the value at a particular
iteration is made up of a combination of the old value
and the newly calculated value:
xinew = λxinew + (1 − λ )xiold
where λ is a weighting factor that is assigned a value
between 0 and 2.
– 0<λ<1: underrelaxation
– λ=1: no relaxation
– 1<λ≤2: overrelaxation
Ec. derivadas parciales
47
Ecuaciones diferenciales en
derivadas parciales
Ec. derivadas parciales
48
Ecuaciones en derivadas parciales
• Ordinary Differential Equations have only one
independent variable
dy
3 + 5 y 2 = 3e − x , y (0) = 5
dx
• Partial Differential Equations have more than one
independent variable
∂ 2u ∂ 2u
3 2 + 2 = x2 + y2
∂x
∂y
• subject to certain conditions: where u is the
dependent variable, and x and y are the independent
variables.
Ec. derivadas parciales
49
Clasificación de EDPs de 2 orden
∂ 2u
∂ 2u
∂ 2u
A 2 +B
+C 2 + D = 0
∂x∂y
∂x
∂y
• where A, B, and C are functions of x and y , and
.D is a function of
∂u ∂u
x, y, u and , .
∂x ∂y
• can be:
Elliptic if B2 – 4AC < 0
Parabolic if B2 – 4AC = 0
Hyperbolic if B2 – 4AC > 0
Ec. derivadas parciales
50
Ejemplos de EDPs de 2 orden
• Elliptic
• Parabolic
A = 1, B = 0, C = 1
∂ 2T ∂ 2T
+ 2 =0
2
∂x
∂y
Laplace equation
A = k , B = 0, C = 0
∂T
∂ T
=k 2
∂t
∂x
2
• Hyperbolic
∂2 y 1 ∂2 y
= 2 2
2
∂x
c ∂t
A = 1, B = 0, C = −
Heat equation
1
c2
Wave equation
Ec. derivadas parciales
51
Ejemplo físico de una PDE elíptica
y
Tt
W
Tr
Tl
x
Tb
L
• Schematic diagram of a plate with specified temperature boundary
conditions
∂ 2T ∂ 2T
• The Laplace equation governs the temperature:
+ 2 =0
2
∂x
Ec. derivadas parciales
∂y
52
Discretizando la PDE elíptica
y
Tt
(0, n)
L
∆x =
m
∆x
∆y
Tl
W
∆y =
n
Tr
∆y
(i − 1, j )
(i, j )
x
(0,0)
Tb
(i, j + 1)
∆x
(i, j )
(i + 1, j )
(i, j − 1)
(m,0)
∂ 2T
T ( x + ∆x, y ) − 2T ( x, y ) + T ( x − ∆x, y )
( x, y ) ≅
2
∂x
(∆x )2
∂ 2T
T ( x, y + ∆y ) − 2T ( x, y ) + T ( x, y − ∆y )
≅
x
y
(
,
)
∂y 2
(∆y )2
Ec. derivadas parciales
53
Discretizando la PDE elíptica
y
Tt
(0, n)
∆x
∆y
Tl
(i, j + 1)
∆x
∆y
Tr
(i − 1, j )
(i, j )
(i, j − 1)
x
(0,0)
Tb
(i + 1, j )
(i, j )
(m,0)
∂ 2T
T ( x + ∆x, y ) − 2T ( x, y ) + T ( x − ∆x, y )
(
x
,
y
)
≅
∂x 2
(∆x )2
∂ 2T
T ( x, y + ∆y ) − 2T ( x, y ) + T ( x, y − ∆y )
(
x
,
y
)
≅
∂y 2
(∆y )2
Ec. derivadas parciales
∂ 2T
∂x 2
∂ 2T
∂y 2
≅
Ti +1, j − 2Ti , j + Ti −1, j
i, j
≅
i, j
(∆x )2
Ti , j +1 − 2Ti , j + Ti , j −1
(∆y )2
54
Discretizando la PDE elíptica
∂ 2T
∂ 2T
+
=0
2
2
∂x
∂y
• Substituting these approximations into the Laplace equation yields:
Ti +1, j − 2Ti , j + Ti −1, j Ti , j +1 − 2Ti , j + Ti , j −1
+
=0
2
2
(∆x )
(∆y )
• if, ∆x = ∆y
• the Laplace equation can be rewritten as
Ti +1, j + Ti −1, j + Ti , j +1 + Ti , j −1 − 4Ti , j = 0
(Eq. 1)
• there are several numerical methods that can be used to solve the
problem:
Direct Method
Gauss-Seidel Method
Lieberman Method
Ec. derivadas parciales
55
Ejemplo 1: Método directo
• Consider a plate 2.4 m × 3.0 m that is subjected to the boundary
conditions shown below. Find the temperature at the interior nodes
using a square grid with a length of 0.6 m by using the direct method.
y
300 °C
W
75 °C
3.0 m
100 °C
50 °C
x
2.4 m
L
Ec. derivadas parciales
56
Ejemplo 1: Método directo
• We discretize the plate by taking, ∆x = ∆y = 0.6m
L
W
= 4 n= =5
m=
∆y
∆x
• The nodal temperatures at the boundary nodes are given by:
y
T 0,5
T0, 4
75 °C
300 °C
T1,5
T2,5
T3,5
T4,5
T1, 4
T2, 4
T3, 4
T4, 4
T1,3
T2,3
T3,3
T4,3
T1, 2
T2, 2
T3, 2
T4, 2
T0,1
T1,1
T2,1
T3,1
T4,1
T0, 0
T1, 0
T2, 0
T3, 0
T4, 0
T0,3
T0, 2
T0, j = 75, j = 1,2,3,4
100 °C
T4, j = 100, j = 1,2,3,4
Ti , 0 = 50, i = 1,2,3
Ti ,5 = 300, i = 1,2,3
x
50 °C
Ec. derivadas parciales
57
Ejemplo 1: Método directo
• the equation for the temperature at the node (2,3)
y
T 0,5
T0, 4
T0,3
T0, 2
T0,1
• i=2 and j=3
T0, 0
T1,5
T2,5
T3,5
T4,5
T1, 4
T2, 4
T3, 4
T4, 4
T1,3
T2,3
T3,3
T4,3
T1, 2
T2, 2
T3, 2
T4, 2
T1,1
T2,1
T3,1
T4,1
x
T1, 0
T2, 0
T3, 0
T4, 0
Ti +1, j + Ti −1, j + Ti , j +1 + Ti , j −1 − 4Ti , j = 0
T3,3 + T1,3 + T2, 4 + T2, 2 − 4T2,3 = 0
T1,3 + T2, 2 − 4T2,3 + T2, 4 + T3,3 = 0
Ec. derivadas parciales
58
Ejemplo 1: Método directo
• We can develop similar equations for every interior node leaving us with
an equal number of equations and unknowns.
• For this problem the number of equations generated is 12
Ec. derivadas parciales
59
Ejemplo 1: Método directo
• The corner nodal temperature of T0,5 , T4,5 , T4,0 , T0,0 are not needed
• To get the temperature at the interior nodes we have to write Equation 1
for all the combinations of i and j, i = 1,...., m − 1; j = 1,...., n − 1
i=1 and j=1
i=1 and j=2
i=1 and j=3
i=1 and j=4
i=2 and j=1
i=2 and j=2
i=2 and j=3
i=2 and j=4
i=3 and j=1
i=3 and j=2
i=3 and j=3
i=3 and j=4
− 4T1,1 + T1, 2 + T2,1 = −125
T1,1 − 4T1, 2 + T1,3 + T2, 2 = −75
T1, 2 − 4T1,3 + T1, 4 + T2,3 = −75
T1,3 − 4T1, 4 + T2, 4 = −375
T1,1 − 4T2,1 + T2, 2 + T3,1 = −50
T1, 2 + T2,1 − 4T2, 2 + T2,3 + T3, 2 = 0
T1,3 + T2, 2 − 4T2,3 + T2, 4 + T3,3 = 0
T1, 4 + T2,3 − 4T2, 4 + T3, 4 = −300
T2,1 − 4T3,1 + T3, 2 = −150
T2, 2 + T3,1 − 4T3, 2 + T3,3 = −100
T2,3 + T3, 2 − 4T3,3 + T3, 4 = −100
T2, 4 + T3,3 − 4T3, 4 = −400
Ec. derivadas parciales
60
Ejemplo 1: Método directo
• We can use Excel and matrix operations to solve the linear equations
system
T1,1 T1,2 T1,3 T1,4 T2,1 T2,2 T2,3 T2,4 T3,1 T3,2 T3,3 T3,4
RHE
-4
1
0
0
1
0
0
0
0
0
0
0
-125
T1,1
74.8719
1
-4
1
0
0
1
0
0
0
0
0
0
-75
T1,2
95.8959
0
1
-4
1
0
0
1
0
0
0
0
0
-75
T1,3
127.8036
0
0
1
-4
1
0
0
1
0
0
0
0
-375
T1,4
196.9288
1
0
0
0
-4
1
0
0
1
0
0
0
-50
T2,1
78.5917
0
1
0
0
1
-4
1
0
0
1
0
0
0
T2,2
105.9082
75.0
95.9 105.9 105.8 100.0
0
0
1
0
0
1
-4
1
0
0
1
0
0
T2,3
143.3896
0
0
0
1
0
0
1
-4
0
0
0
1
-300
T2,4
206.3200
75.0
74.9
78.6
83.6 100.0
0
0
0
0
1
0
0
0
-4
1
0
0
-150
T3,1
83.5868
50.0
50.0
50.0
0
0
0
0
0
1
0
0
1
-4
1
0
-100
T3,2
105.7554
0
0
0
0
0
0
1
0
0
1
-4
1
-100
T3,3
133.5267
0
0
0
0
0
0
0
1
0
0
1
-4
-400
T3,4
184.9617
Ec. derivadas parciales
300.0 300.0 300.0
75.0 196.9 206.3 185.0 100.0
75.0 127.8 143.4 133.5 100.0
61
Método Gauss-Seidel
• Recall the discretized equation
Ti +1, j + Ti −1, j + Ti , j +1 + Ti , j −1 − 4Ti , j = 0
• This can be rewritten as
Ti , j =
Ti +1, j + Ti −1, j + Ti , j +1 + Ti , j −1
4
• For the Gauss-Seidel Method, this equation is solved
iteratively for all interior nodes until a pre-specified
tolerance is met.
Ec. derivadas parciales
62
Ejemplo 2: Método Gauss-Seidel
• Consider a plate 2.4 m × 3.0 m that is subjected to the boundary
conditions shown below. Find the temperature at the interior nodes
using a square grid with a length of 0.6 m using the Gauss-Siedel
method. Assume the initial temperature at all interior nodes to be 0 °C .
y
W
300 °C
75 °C
3.0 m
100 °C
50 °C
x
2.4 m
L
63
Ejemplo 2: Método Gauss-Seidel
• Discretizing the plate by taking, ∆x = ∆y = 0.6m
L
W
= 4 n= =5
m=
∆y
∆x
• The nodal temperatures at the boundary nodes are given by:
y
T 0,5
T0, 4
75 °C
300 °C
T1,5
T2,5
T3,5
T4,5
T1, 4
T2, 4
T3, 4
T4, 4
T1,3
T2,3
T3,3
T4,3
T1, 2
T2, 2
T3, 2
T4, 2
T0,1
T1,1
T2,1
T3,1
T4,1
T0, 0
T1, 0
T2, 0
T3, 0
T4, 0
T0,3
T0, 2
T0, j = 75, j = 1,2,3,4
100 °C
T4, j = 100, j = 1,2,3,4
Ti , 0 = 50, i = 1,2,3
Ti ,5 = 300, i = 1,2,3
x
50 °C
Ec. derivadas parciales
64
Ejemplo 2: Método Gauss-Seidel
• Now we can begin to solve for the temperature at each interior node
using
T +T +T +T
Ti , j =
i +1, j
i −1, j
i , j +1
4
i , j −1
, i = 1,2,3,4; j = 1,2,3,4,5
• Assume all internal nodes to have an initial temperature of zero.
• Iteration 1:
i=1 and j=1 T1,1 = 31.25 º C
i=2 and j=3 T2,3 = 9.27735 º C
i=1 and j=2 T1, 2 = 26.5625 º C
i=2 and j=4 T2, 4 = 102.344 º C
i=1 and j=3 T1,3 = 25.3906 º C
i=3 and j=1 T3,1 = 42.5781 º C
i=1 and j=4 T1, 4 = 100.098 º C
i=3 and j=2 T3, 2 = 38.5742 º C
i=2 and j=1 T2,1 = 20.3125 º C
i=3 and j=3 T3,3 = 36.9629 º C
i=2 and j=2 T2, 2 = 11.7188 º C
i=3 and j=4 T3, 4 = 134.827 º C
Ec. derivadas parciales
65
Ejemplo 2: Método Gauss-Seidel
• Iteration 2:
we take the temperatures from iteration 1 and calculate the
present
previous
approximated error. ε = Ti, j − Ti, j ×100
a i, j
Ti ,present
j
i=1, j=1 T1,1 = 42.9688 º C ε a 1,1 = 27.27%
i=2, j=3 T2,3 = 56.4881 º C ε a 2,3 = 83.58%
i=1, j=2 T1, 2 = 38.7596 º C ε a 1, 2 = 31.49% i=2, j=4 T2, 4 = 156.150 º C ε a 2, 4 = 34.46%
i=1, j=3 T1,3 = 55.7862 º C ε a 1,3 = 54.49% i=3, j=1 T3,1 = 56.3477 º C ε a 3,1 = 24.44%
i=1, j=4 T1, 4 = 133.283 º C ε a 1, 4 = 24.90% i=3, j=2 T3, 2 = 56.0425 º C ε a 3, 2 = 31.70%
i=2, j=1 T2,1 = 36.8164 º C ε a 2,1 = 44.83% i=3, j=3 T3,3 = 86.8394 º C ε a 3,3 = 57.44%
i=2, j=2 T2, 2 = 30.8594 º C ε a = 62.03% i=3, j=4 T3, 4 = 160.747 º C ε a 3, 4 = 16.12%
2, 2
Ec. derivadas parciales
66
Ejemplo 2: Método Gauss-Seidel
Node
Temperature Distribution in the Plate (°C)
Number of Iterations
1
2
10
T1,1
31.2500
42.9688
73.0239
T1, 2
T1,3
26.5625
38.7695
91.9585
25.3906
55.7861
119.0976
T1, 4
T2,1
T2, 2
T2,3
T2, 4
100.0977
133.2825
172.9755
20.3125
36.8164
76.6127
11.7188
30.8594
102.1577
9.2773
56.4880
137.3802
102.3438
156.1493
198.1055
T3,1
T3, 2
T3,3
T3, 4
42.5781
56.3477
82.4837
38.5742
56.0425
103.7757
36.9629
86.8393
130.8056
134.8267
160.7471
182.2278
Ec. derivadas parciales
67
Ejemplo 2: Método Gauss-Seidel in Excel
• The numerical solution of Laplace equation at a point is the average of
four neighbors
Ti +1, j + Ti −1, j + Ti , j +1 + Ti , j −1
Ti , j =
• Example for cell S8:
=(S7+S9+R8+T8)/4
4
• Enter the boundary conditions in the appropriate cells.
• Copy and paste to cover the cells where values of the potential are to
be calculated. This calculation contains a "circular reference“.
Ec. derivadas parciales
68
Ejemplo 2: Método Gauss-Seidel in Excel
• To allow circular references and enable iterations:
File → Options → Formulas
On the "Calculations options" form select "Enable iterative calculation"
We can increase the Maximum Iterations (100 is the deafult) and
reduce the Maximum Change (0.001 is the default). Iterations will stop
when the maximum iteration is reached or the change is less than the
maximum change.
• F9 to recalculate.
Ec. derivadas parciales
69
Ejemplo 2: Método Gauss-Seidel in Excel
• Color cell based on value
To achieve the cell color based on value: Inicio → Estilos →
Formato condicional → Escalas de color → Más reglas
We can chose a 3 color scale with blue for mimimum, white or gray for
midpoint and red for maximum.
Ec. derivadas parciales
70
Ejemplo 2: Método Gauss-Seidel in Excel
• Plotting the results
Normally we use the chart type Surface or Contour.
Ec. derivadas parciales
71
Ejemplo 3: Uso de Solver
Finite Difference Solution
x\y
8
7
6
5
4
3
2
1
0
100
100
100
100
100
100
100
0
x\y
sum =
0
33.90
56.89
71.80
81.62
88.31
93.11
96.80
100
2
0
26.95
48.82
65.03
76.65
84.97
91.07
95.85
100
3
0
25.07
46.39
62.86
75.00
83.83
90.37
95.52
100
4
0
26.95
48.82
65.03
76.65
84.97
91.07
95.85
100
5
0
33.90
56.89
71.80
81.62
88.31
93.11
96.80
100
6
0
51.74
73.07
83.64
89.71
93.56
96.23
98.26
100
7
100
100
100
100
100
100
100
8
=(-4*D5+D4+D6+C5+E5)^2
Residuals-squared
8
7
6
5
4
3
2
1
0
0
51.74
73.07
83.64
89.71
93.56
96.23
98.26
100
1
0
2.7745E-13
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
1
2
3
4
5
6
7
8
=SUMA(D17:J23)
Ec. derivadas parciales
72
Método de Lieberman
• Recall the equation used in the Gauss-Siedel Method,
Ti , j =
Ti +1, j + Ti −1, j + Ti , j +1 + Ti , j −1
4
• Because the Gauss-Siedel Method is guaranteed to converge, we can
accelerate the process by using overrelaxation. In this case,
old
Ti ,relaxed
= λTi ,new
j + (1 − λ )Ti , j
j
• The λ is known as the “overrelaxation parameter" and is in the range
0 < λ < 2.
Ec. derivadas parciales
73
Condiciones de contorno alternativas
• In the past examples, the boundary conditions on the plate had a
specified temperature on each edge.
What if the conditions are different ? For example, what if one of the
edges of the plate is insulated.
• In this case, the boundary condition would be the derivative of the
temperature. Because if the right edge of the plate is insulated, then
the temperatures on the right edge nodes also become unknowns.
y
300 °C
75 °C
Insulated
3.0 m
50 °C
2.4 m
Ec. derivadas parciales
x
74
Condiciones de contorno alternativas
• The finite difference equation in this case for the right edge for the
nodes (m, j ) for j = 2,3,..n − 1
Tm +1, j + Tm −1, j + Tm , j −1 + Tm , j +1 − 4Tm , j = 0
• However the node (m + 1, j ) is not inside the plate. The derivative
boundary condition needs to be used to account for these additional
unknown nodal temperatures on the right edge. This is done by
approximating the derivative at the edge node (m, j ) as
∂T
∂x
≅
m, j
y
Tm +1, j − Tm −1, j
300 °C
2(∆x)
75 °C
Insulated
3.0 m
50 °C
2.4 m
Ec. derivadas parciales
x
75
Condiciones de contorno alternativas
• Rearranging this approximation gives us,
Tm +1, j = Tm −1, j + 2(∆x)
∂T
∂x
m, j
• We can then substitute this into the original equation gives us,
2Tm −1, j + 2(∆x)
∂T
∂x
+ Tm , j −1 + Tm , j +1 − 4Tm , j = 0
m, j
• Recall that is the edge is insulated then,
∂T
∂x
• Substituting this again yields,
=0
m, j
2Tm −1, j + Tm , j −1 + Tm , j +1 − 4Tm , j = 0
Ec. derivadas parciales
76
Ecuaciones en derivadas parciales
parabólicas
• The general form for a second order linear PDE with
two independent variables and one dependent
variable2 is
2
2
∂u
∂u
∂u
∂u
∂u
+ Fu + G = 0
+C 2 + D + E
A 2 +B
∂y
∂x
∂y
∂x∂y
∂x
• The criteria for an equation of this type to be
considered parabolic: B 2 − 4 AC = 0
• Examine the
heat-conduction
equation
given
by
2
∂ T ∂T
k k = thermal conductivity of rod material,
where α =
α 2 =
ρ = density of rod material,
ρ
C
∂t
∂x
C = specific heat of the rod material.
where A = α , B = 0, C = 0, D = 0, E = −1, F = 0, G = 0
thus we can classify this equation as parabolic.
Ec. derivadas parciales
77
Ejemplo de una EDP parabólica
• Consider the flow of heat within a metal rod of length
L, one end of which is held at a known high
temperature, the other end at a lower temperature.
– Heat will flow from the hot end to the cooler end.
– We'll assume that the rod is perfectly insulated, so that
heat loss through the sides can be neglected.
• We want to calculate the temperature along the length
of the rod as a function of time.
Ec. derivadas parciales
78
Discretización de una EDP Parabólica
• For a rod of length L divided into n + 1 nodes ∆x = L
n
• The time is similarly broken into time steps of ∆t
j
T
• Hence i corresponds to the temperature at node i ,that is,
x = (i )(∆x ) and time t = ( j )(∆t )
x
∆x
∆x
i
i −1
i +1
Schematic diagram showing interior nodes
Ec. derivadas parciales
79
Solución EDP Parabólica:
Método explícito
• If we define
∆x =
L
n
we can then write the finite central divided
difference approximation of the left hand side at a general interior node
j
j
j
2
( i ) as ∂ T ≅ Ti +1 − 2Ti + Ti −1 where ( j ) is the node number along
∂x 2
(∆x )2
the time.
• The time derivative on the right hand side is approximated by the
forward divided difference method as, ∂T
Ti j +1 − Ti j
≅
∂t i , j
∆t
i, j
x
∆x
∆x
i
i −1
i +1
Schematic diagram showing interior nodes
Ec. derivadas parciales
80
Solución EDP Parabólica:
Método explícito
• Substituting these approximations into the governing equation yields
α
Ti +j1 − 2Ti j + Ti −j1
(∆x )2
Ti j +1 − Ti j
=
∆t
• Solving for the temp at the time node
Ti
j +1
j +1
gives
∆t
j
j
j
= Ti + α
T
−
2
+
T
T
i +1
i
i −1
(∆x) 2
j
(
• choosing, λ = α ∆t
2
(∆x)
(
)
j +1
j
j
j
j
• we can write the equation as, Ti = Ti + λ Ti +1 − 2Ti + Ti −1
)
• we can be solved explicitly: for each internal location node of the rod for
time node j + 1 in terms of the temperature at time node j. If we know
the temperature at node j = 0 , and the boundary temperatures, we can
find the temperature at the next time step. We continue the process until
we reach the time at which we are interested in finding the temperature.
Ec. derivadas parciales
81
Ejemplo 1 EDP Parabólica:
Método explícito
• Consider a steel rod that is subjected to a temperature of 100°C on the
left end and 25°C on the right end. If the rod is of length 0.05m ,use the
explicit method to find the temperature distribution in the rod from t = 0
and t = 9 seconds. Use ∆x = 0.01m , ∆t = 3s .
• Given: k = 54
W
kg
J
, ρ = 7800 3 , C = 490
kg − K
m−K
m
• The initial temperature of the rod is 20°C.
i=0
1
2
3
T =100 ° C
4
5
T = 25 °C
0.01m
Ec. derivadas parciales
82
Ejemplo 1 EDP Parabólica:
Método explícito
• Number of time steps =
t final − tinitial
∆t
=
9−0
=3
3
k
54
=
= 1.4129 ×10 −5 m 2 / s
ρC 7800 × 490
Then, λ = α ∆t 2 = 1.4129 ×10−5 3 2 = 0.4239
(∆x )
(0.01)
• Recall, α =
•
• Boundary Conditions T0 j = 100°C 
 for all j = 0,1,2,3
T5 = 25°C 
• All internal nodes are at 20°C for t = 0 sec: Ti 0 = 20°C , for all i = 1,2,3,4
j
T00 = 100°C
T10 = 20°C 

T20 = 20°C 
 Interior nodes
T30 = 20°C 
T40 = 20°C 
We can now calculate the temperature at each node
explicitly using the equation formulated earlier,
(
Ti j +1 = Ti j + λ Ti +j1 − 2Ti j + Ti −j1
)
T50 = 25°C
Ec. derivadas parciales
83
Ejemplo 1 EDP Parabólica:
Método explícito
• Nodal temperatures vs. Time
t = 0 sec j = 0
T00 = 100°C
T10 = 20°C 

T20 = 20°C 
 Interior nodes
T30 = 20°C 
T40 = 20°C 
T50 = 25°C
t = 6 sec j = 2
t = 9 sec j = 3
T01 = 100°C − Boundary Condition
T02 = 100°C
T03 = 100°C
T11 = 53.912°C
T12 = 59.073°C
T13 = 65.953°C
T21 = 20°C
T22 = 34.375°C
T23 = 39.132°C
T31 = 20°C
T32 = 20.889°C
T33 = 27.266°C
T41 = 22.120°C
T42 = 22.442°C
T43 = 22.872°C
t = 3 sec j = 1
T51 = 25°C − Boundary Condition T52 = 25°C
Ec. derivadas parciales
T53 = 25°C
84
Ejemplo 2 EDP Parabólica:
Método explícito
50
Time-dependent Temperature Distribution in a Brass Rod
(Temperature values in bold are constant)
10
0,09
(hcap)
thermal conductivity of brass, cal/sec/cm/deg
0,26
(k)
density of brass, g/cm3
8,4
(rho)
Coefficient e in general PDE, =k/(hcap*rho)
0,33(e)
∆x
∆t
f=e*Dt/(Dx^2)
1
(Dx)
1
(Dt)
0,33
(f)
time t (sec)
Distance x (cm)
0
1
2
3
4
5
6
7
8
9
10
0
100
0
0
0
0
0
0
0
0
0
0
1
100
32,9
0,0
0,0
0,0
0,0
0,0
0,0
0,0
0,0
0
2
100
44,2
10,8
0,0
0,0
0,0
0,0
0,0
0,0
0,0
0
3
100
51,6
18,2
3,6
0,0
0,0
0,0
0,0
0,0
0,0
0
4
100
56,5
24,4
7,2
1,2
0,0
0,0
0,0
0,0
0,0
0
5
100
60,3
29,3
10,9
2,8
0,4
0,0
0,0
0,0
0,0
0
6
100
63,2
33,4
14,3
4,7
1,0
0,1
0,0
0,0
0,0
0
7
100
65,5
36,9
17,4
6,6
1,9
0,4
0,0
0,0
0,0
0
8
100
67,5
39,9
20,3
8,6
3,0
0,8
0,1
0,0
0,0
0
9
100
69,1
42,5
22,9
10,6
4,1
1,3
0,3
0,1
0,0
0
10
100
70,5
44,8
25,3
12,5
5,3
1,9
0,5
0,1
0,0
0
11
100
71,8
46,9
27,5
14,4
6,6
2,6
0,9
0,2
0,0
0
12
100
72,9
48,7
29,6
16,1
7,8
3,3
1,2
0,4
0,1
0
13
100
73,8
50,4
31,4
17,8
9,1
4,1
1,6
0,6
0,2
0
14
100
74,7
51,9
33,2
19,4
10,3
4,9
2,1
0,8
0,2
0
15
100
75,5
53,2
34,8
21,0
11,5
5,8
2,6
1,0
0,3
0
16
100
76,2
54,5
36,3
22,4
12,7
6,6
3,1
1,3
0,5
0
17
100
76,9
55,7
37,7
23,8
13,9
7,5
3,7
1,6
0,6
0
18
100
77,5
56,8
39,1
25,1
15,1
8,4
4,3
2,0
0,7
0
19
100
78,1
57,8
40,3
26,4
16,2
9,2
4,9
2,3
0,9
0
20
100
78,6
58,7
41,5
27,6
17,2
10,1
5,5
2,7
1,1
0
21
100
79,1
59,6
42,6
28,8
18,3
10,9
6,1
3,1
1,2
0
22
100
79,6
60,4
43,6
29,9
19,3
11,7
6,7
3,5
1,4
0
23
100
80,0
61,2
44,6
30,9
20,3
12,6
7,3
3,8
1,6
0
24
100
80,4
61,9
45,6
31,9
21,2
13,4
7,9
4,2
1,8
0
25
100
80,8
62,6
46,5
32,9
22,2
14,2
8,5
4,6
2,0
0
40
Temperature ºC
heat capacity of brass, cal/g/deg
30
Temperature at x=5cm.
20
10
0
0
20
40
60
Time sec.
80
100
Temperature distribution along
the length of the rod
120
100
Temperature ºC
length, cm
80
60
t = 20 sec.
40
t = 50 sec.
20
t = 100 sec.
0
0
Ec. derivadas parciales
2
4
6
8
Loc ation on rod cm.
10
85
Solución EDP Parabólica:
Método implícito
• Using the explicit method, we were able to find the temperature at each
node, one equation at a time.
• However, the temperature at a specific node was only dependent on the
temperature of the neighboring nodes from the previous time step. This
is contrary to what we expect from the physical problem.
• The implicit method allows us to solve this and other problems by
developing a system of simultaneous linear equations for the
temperature at all interior nodes at a particular time.
• The second derivative is approximated by the CDD and the first
derivative by the BDD scheme at time level j+1 at node ( i ) as
∂ 2T ∂T
α 2 =
∂x
∂t
Ti +j1+1 − 2Ti j +1 + Ti −j1+1
∂ 2T
≈
2
∂x i , j +1
(∆x )2
Ti j +1 − Ti j
∂T
≈
∂t i , j +1
∆t
Ec. derivadas parciales
86
Solución EDP Parabólica:
Método implícito
• Substituting these approximations into the heat conduction equation
yields
Ti +j1+1 − 2Ti j +1 + Ti −j1+1 Ti j +1 − Ti j
α
=
2
∆t
(∆x )
∂ 2T ∂T
α 2 =
∂t
∂x
• Rearranging yields
− λTi −j1+1 + (1 + 2λ )Ti j +1 − λTi +j1+1 = Ti j
given that
λ =α
∆t
(∆x )2
• The rearranged equation can be written for every node during each
time step. These equations can then be solved as a simultaneous
system of linear equations to find the nodal temperatures at a particular
time.
Ec. derivadas parciales
87
Ejemplo 2 EDP Parabólica:
Método implícito
• Consider a steel rod that is subjected to a temperature of 100°C on the
left end and 25°C on the right end. If the rod is of length 0.05m ,use the
implicit method to find the temperature distribution in the rod from t = 0
and t = 9 seconds. Use ∆x = 0.01m , ∆t = 3s .
• Given: k = 54
W
kg
J
, ρ = 7800 3 , C = 490
kg − K
m−K
m
• The initial temperature of the rod is 20°C.
i=0
1
2
3
T =100 ° C
4
5
T = 25 °C
0.01m
Ec. derivadas parciales
88
Ejemplo 2 EDP Parabólica:
Método implícito
• Number of time steps =
t final − tinitial
∆t
=
9−0
=3
3
k
54
=
= 1.4129 ×10 −5 m 2 / s
ρC 7800 × 490
Then, λ = α ∆t 2 = 1.4129 ×10−5 3 2 = 0.4239
(∆x )
(0.01)
• Recall, α =
•
• Boundary Conditions T0 j = 100°C 
 for all j = 0,1,2,3
T5 = 25°C 
• All internal nodes are at 20°C for t = 0 sec: Ti 0 = 20°C , for all i = 1,2,3,4
j
T00 = 100°C
T10 = 20°C 

T20 = 20°C 
 Interior nodes
T30 = 20°C 
T40 = 20°C 
T50 = 25°C
We can now form the system of equations for the
first time step by writing the approximated heat
conduction equation for each node
− λTi −j1+1 + (1 + 2λ )Ti j +1 − λTi +j1+1 = Ti j
Ec. derivadas parciales
89
Ejemplo 2 EDP Parabólica:
Método implícito
• Nodal temperatures when t = 3 sec
• For the first time step we can write four such equations with four
unknowns, expressing them in matrix form yields
0
0
 1.8478 − 0.4239
 T11  62.390
− 0.4239 1.8478 − 0.4239
  1   20 
0


 T2  = 
1

0
− 0.4239 1.8478 − 0.4239 T3   20 


 1  
0
0
−
0
.
4239
1
.
8478

 T4  30.598
• The above coefficient matrix is tri-diagonal, so special algorithms
(e.g.Thomas’ algorithm) can be used to solve. The solution is given by
T11   39.451
 1 

T2  =  24.792
T31   21.438
 1 

T4  21.477 
T01   100 
 1 

T1   39.451
T21   24.792
 1 = 

T3   21.438
T 1  21.477 
 41  

T5   25 
Ec. derivadas parciales
90
Ejemplo 2 EDP Parabólica:
Método implícito
• Nodal temperatures when: t = 3 sec
T01   100 
 1 

39
.
451
T
1
  

1
T2  24.792
 1 = 

21
.
438
T
 3 

1
T  21.477 
 41  

25
T
 5  

t = 6 sec
T02   100 
 2 

51
.
326
T
1
  

2
T2  30.669 
 2 = 

23
.
876
T
 3 

2
T  22.836
 42  

25
T


 5  
Ec. derivadas parciales
t = 9 sec
T03   100 
 3 

.
043
59
T
1
  

T23  36.292 
 3 = 

T3  26.809
T 3   24.243
 43  

25
T


 5  
91
Solución EDP Parabólica:
Método Crank-Nicolson
2
• Using the implicit method our approximation of ∂ T was of accuracy
O(∆x) 2, while our approximation of ∂T was of
∂t
∂x 2
O(∆t )
accuracy.
• One can achieve similar orders of accuracy by approximating the
second derivative, on the left hand side of the heat equation, at the
midpoint of the time step. Doing so yields
∂ 2T
∂x 2
≈
i, j
α  Ti +j1 − 2Ti j + Ti −j1

2
(∆x )2
Ti +j1+1 − 2Ti j +1 + Ti −j1+1 
+

(∆x )2

• The first derivative, on the right hand side of the heat equation, is
approximated using the forward divided difference method at time level
j +1,
j +1
j
∂T
∂t
≈
i, j
Ti
− Ti
∆t
Ec. derivadas parciales
92
Solución EDP Parabólica:
Método Crank-Nicolson
• Substituting these approximations into the governing equation for heat
conductance yields
α  Ti +j1 − 2Ti j + Ti −j1

2
(∆x )2
Ti +j1+1 − 2Ti j +1 + Ti −j1+1  Ti j +1 − Ti j
+
=
2
∆t
(∆x )

• giving
• where
− λTi −j1+1 + 2(1 + λ )Ti j +1 − λTi +j1+1 = λTi −j1 + 2(1 − λ )Ti j + λTi +j1
λ =α
∆t
(∆x )2
• Having rewritten the equation in this form allows us to discretize the
physical problem. We then solve a system of simultaneous linear
equations to find the temperature at every node at any point in time.
Ec. derivadas parciales
93
Ejemplo 3 EDP Parabólica:
Método Crank-Nicolson
• Consider a steel rod that is subjected to a temperature of 100°C on the
left end and 25°C on the right end. If the rod is of length 0.05m ,use the
Crank-Nicolson method to find the temperature distribution in the rod
from t = 0 and t = 9 seconds. Use ∆x = 0.01m, ∆t = 3s .
• Given: k = 54
W
kg
J
, ρ = 7800 3 , C = 490
kg − K
m−K
m
• The initial temperature of the rod is 20°C.
i=0
1
2
3
T =100 ° C
4
5
T = 25 °C
0.01m
Ec. derivadas parciales
94
Ejemplo 3 EDP Parabólica:
Método Crank-Nicolson
• Number of time steps =
t final − tinitial
∆t
=
9−0
=3
3
k
54
=
= 1.4129 ×10 −5 m 2 / s
ρC 7800 × 490
Then, λ = α ∆t 2 = 1.4129 ×10−5 3 2 = 0.4239
(∆x )
(0.01)
• Recall, α =
•
• Boundary Conditions T0 j = 100°C 
 for all j = 0,1,2,3
T5 = 25°C 
• All internal nodes are at 20°C for t = 0 sec: Ti 0 = 20°C , for all i = 1,2,3,4
j
T00 = 100°C
We can now form the system of equations for the
T10 = 20°C 

first time step by writing the approximated heat
T20 = 20°C 
 Interior nodes conduction equation for each node
T30 = 20°C 
T40 = 20°C 
− λTi −j1+1 + 2(1 + λ )Ti j +1 − λTi +j1+1 = λTi −j1 + 2(1 − λ )Ti j + λTi +j1
T50 = 25°C
Ec. derivadas parciales
95
Ejemplo 3 EDP Parabólica:
Método Crank-Nicolson
• Nodal temperatures when t = 3 sec
• For the first time step we can write four such equations with four
unknowns, expressing them in matrix form yields
0
0  T11  116.30 
 2.8478 − 0.4239
  1  40.000
− 0.4239 2.8478 − 0.4239
0

 T2  = 

1
 0
− 0.4239 2.8478 − 0.4239 T3  40.000

 1  

0
0
−
0
.
4239
2
.
8478
52
.
718
T


  4  
• The above coefficient matrix is tri-diagonal, so special algorithms
(e.g.Thomas’ algorithm) can be used to solve. The solution is given by
T11   44.372
 1 

23
.
746
T
2
 =

T31  20.797 
 1 

T4  21.607 
T01   100 
 1 

T1   44.372
T21   23.746
 1 = 

20
.
797
T
 3 

1
T  21.607 
 41  

T5   25 
Ec. derivadas parciales
96
Ejemplo 3 EDP Parabólica:
Método Crank-Nicolson
• Nodal temperatures when: t = 3 sec
T01   100 
 1 

44
.
372
T
1
  

1
T2   23.746
 1 = 

20
.
797
T
 3 

1
T  21.607 
 41  

25
T


 5  

t = 6 sec
t = 9 sec
T02   100 
 2 

55
.
883
T
1
  

2
T2  31.075
 2 = 

23
.
174
T
 3 

2
T  22.730
 42  

25
T


 5  
T03   100 
 3 

62
.
604
T
1
  

3
T2  37.613
 3 = 

26
.
562
T
 3 

3
T  24.042
 43  

T5   25 
Ec. derivadas parciales
97
Comparación de métodos:
temperaturas en 9 seg.
• The table below allows you to compare the results from all three
methods discussed in juxtaposition with the analytical solution.
Node
Explicit
Implicit
T13
T23
T33
T43
65.953
39.132
59.043
36.292
27.266
22.872
26.809
24.243
Ec. derivadas parciales
CrankNicolson
62.604
37.613
26.562
24.042
Analytical
62.510
37.084
25.844
23.610
98
Ecuaciones en derivadas parciales
hiperbólicas
• The general form for a second order linear PDE with
two independent variables and one dependent
variable2 is
2
2
∂u
∂u
∂u
∂u
∂u
+ Fu + G = 0
+C 2 + D + E
A 2 +B
∂y
∂x
∂y
∂x∂y
∂x
• The criteria for an equation of this type to be
considered hyperbolic: B 2 − 4 AC > 0
• The wave equation (oscillatory systems) given by
∂2 y
∂2 y
Tg
k
=
where
k
=
w
∂x 2
∂t 2
where A = 1, B = 0, C = −1
T = tension,
g = gravitational constant,
w = weight/unit =W/L ,W=weight, L=length
thus we can classify this equation as hyperbolic.
Ec. derivadas parciales
99
Ejemplo de una EDP hiperbólica
• A string of certain length and weight is under a fixed
tension. Initially the mid-point of the string is displaced
some distance from its equilibrium position and
released.
• We want to calculate the displacement as a function
of time at fixed intervals along the length of the string.
Ec. derivadas parciales
100
Solución EDP hiperbólica:
Método explícito
• Once again, we can solve the problem by replacing derivatives by finite
differences.
j +1
j
j −1
j
j
 j

Ti
− 2Ti + Ti
(∆t )2
• which, when rearranged, yields
=
Tg Ti +1 − 2Ti + Ti −1


w 
(∆x )2

2
 Tg (∆t )2  j
Tg (∆t )
j
j
j −1
T
(Ti +1 + Ti −1 ) − Ti + 21 −
Ti =
2
2  i
w (∆x )
w (∆x ) 

• If we set Tg (∆t )2 w(∆x )2 = 1 , the above equation is simplified to
j +1
Ti j +1 = Ti +j1 + Ti −j1 − Ti j −1
• When employing the simplified equation, the value of ∆t is determined
by the expression ∆t = ∆x .To begin the calculations (value at t1),
Tg / w
it is required values of the function at t0 = 0 and also a value at t-1. We
can get a value for the function at t-1 by making use of the fact that the
0
0
T
+
T
1
i
+
1
i
−
1 for the first row.
function is periodic. We can use Ti =
2
Ec. derivadas parciales
101
Ejemplo 1 EDP hiperbólica
• A string 50 cm long and weighing 0.5 g is under a tension of 33 kg.
Initially the mid-point of the string is displaced 0.5 cm from its
equilibrium position and released. We want to calculate the
displacement as a function of time at 5 cm intervals along the length of
the string, using equation
Ti j +1 = Ti +j1 + Ti −j1 − Ti j −1
• From equation ∆t =
∆x
the ∆t must be 8.8 x 10-5 seconds.
Tg / w
Ec. derivadas parciales
102
Ejemplo 1 EDP hiperbólica
The Wave Equation:Vibration of a String
time t (sec)
gravitational constant, cm/sec2
Dx
Dt
0
8,8E-05
1,8E-04
2,6E-04
3,5E-04
4,4E-04
5,3E-04
6,2E-04
7,0E-04
7,9E-04
8,8E-04
9,7E-04
1,1E-03
1,1E-03
1,2E-03
1,3E-03
1,4E-03
1,5E-03
1,6E-03
1,7E-03
1,8E-03
1,8E-03
1,9E-03
2,0E-03
2,1E-03
2,2E-03
2,3E-03
2,4E-03
2,5E-03
2,5E-03
2,6E-03
2,7E-03
2,8E-03
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
5
0,1
0,1
0,1
0,1
0,1
0,0
-0,1
-0,1
-0,1
-0,1
-0,1
-0,1
-0,1
-0,1
-0,1
0,0
0,1
0,1
0,1
0,1
0,1
0,1
0,1
0,1
0,1
0,0
-0,1
-0,1
-0,1
-0,1
-0,1
-0,1
-0,1
10
0,2
0,2
0,2
0,2
0,1
0,0
-0,1
-0,2
-0,2
-0,2
-0,2
-0,2
-0,2
-0,2
-0,1
0,0
0,1
0,2
0,2
0,2
0,2
0,2
0,2
0,2
0,1
0,0
-0,1
-0,2
-0,2
-0,2
-0,2
-0,2
-0,2
15
0,3
0,3
0,3
0,2
0,1
0,0
-0,1
-0,2
-0,3
-0,3
-0,3
-0,3
-0,3
-0,2
-0,1
0,0
0,1
0,2
0,3
0,3
0,3
0,3
0,3
0,2
0,1
0,0
-0,1
-0,2
-0,3
-0,3
-0,3
-0,3
-0,3
Distance x (cm)
20
25
30
0,4
0,5
0,4
0,4
0,4
0,4
0,3
0,3
0,3
0,2
0,2
0,2
0,1
0,1
0,1
0,0
0,0
0,0
-0,1
-0,1
-0,1
-0,2
-0,2
-0,2
-0,3
-0,3
-0,3
-0,4
-0,4
-0,4
-0,4
-0,5
-0,4
-0,4
-0,4
-0,4
-0,3
-0,3
-0,3
-0,2
-0,2
-0,2
-0,1
-0,1
-0,1
0,0
0,0
0,0
0,1
0,1
0,1
0,2
0,2
0,2
0,3
0,3
0,3
0,4
0,4
0,4
0,4
0,5
0,4
0,4
0,4
0,4
0,3
0,3
0,3
0,2
0,2
0,2
0,1
0,1
0,1
0,0
0,0
0,0
-0,1
-0,1
-0,1
-0,2
-0,2
-0,2
-0,3
-0,3
-0,3
-0,4
-0,4
-0,4
-0,4
-0,5
-0,4
-0,4
-0,4
-0,4
-0,3
-0,3
-0,3
50
33000
0,5
0,01
(L)
(T)
(Wt)
(w)
980
5
8,79E-05
(g)
(Dx)
(Dt)
35
0,3
0,3
0,3
0,2
0,1
0,0
-0,1
-0,2
-0,3
-0,3
-0,3
-0,3
-0,3
-0,2
-0,1
0,0
0,1
0,2
0,3
0,3
0,3
0,3
0,3
0,2
0,1
0,0
-0,1
-0,2
-0,3
-0,3
-0,3
-0,3
-0,3
40
0,2
0,2
0,2
0,2
0,1
0,0
-0,1
-0,2
-0,2
-0,2
-0,2
-0,2
-0,2
-0,2
-0,1
0,0
0,1
0,2
0,2
0,2
0,2
0,2
0,2
0,2
0,1
0,0
-0,1
-0,2
-0,2
-0,2
-0,2
-0,2
-0,2
45
0,1
0,1
0,1
0,1
0,1
0,0
-0,1
-0,1
-0,1
-0,1
-0,1
-0,1
-0,1
-0,1
-0,1
0,0
0,1
0,1
0,1
0,1
0,1
0,1
0,1
0,1
0,1
0,0
-0,1
-0,1
-0,1
-0,1
-0,1
-0,1
-0,1
50
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Displacement evolution at 25 cm
0.6
0.4
Displacement (cm.)
length, cm
tension, g
weight,g
weight per unit length, g/cm
0.2
0
-0.2
0
0.0005
0.001
0.0015
0.002
0.0025
0.003
-0.4
-0.6
Ec. derivadas parciales
Time (sec.)
103

Summary of the model specified

Meisterschaft NLB / Championnat LNB