Problemas de frontera para
ecuaciones diferenciales
Cálculo Numérico
Práctica 3
Algoritmos de resolución

Algoritmos auxiliares
• Problema de valor inicial
• Sistemas lineales tridiagonales

Métodos de disparo
• Problema lineales y no lineales

Métodos de diferencias finitas
• Problema lineales y no lineales

Métodos variacionales
• Elementos lineales y cúbicos
Problema de valor inicial
y  ( t )  f ( t , y ( t ))
,
t [a , b]
,
y(t0 )  y0
siendo,
y ( t )  ( y 1 ( t ), y 2 ( t ), . . . , y m ( t ))
T
y  ( t )  ( y 1 ( t ), y 2 ( t ), . . . , y m ( t ))
T
f ( t , y )  ( f 1 ( t , y ), f 2 ( t , y ), . . . , f m ( t , y ))
y 0  ( y 01 , y 02 , . . . , y 0 m )
T
T
Algoritmo de MATLAB 5.2 para
el Problema de Valor Inicial

[X,Y] = ode23(‘fun’, [a, b], y0)
• fun.m: fichero que calcula las derivadas
function z=f(x,y)
% Ecuación no lineal de segundo orden
%
y" = (32 + 2x^3 -yy')/8
z=[y(2);(32 + 2*x^3 - y(1)*y(2))/8];

[T,Y] = ode23('fun',[1,3],[17,0])
Sistemas lineales tridiagonales
 a1

 c1







b1
a2
b2
c2
a3
b3



cn  2
a n 1
c n 1
  x1   d 1 

 

 x2   d 2 
 x   d 
 3    3 
     

 

bn 1 x n 1
d n 1

 

a n   x n   d n 
Algoritmo de Crout











% Eliminación
l(1)=a(1); y(1)=d(1)/l(1);
for i=2:n
u(i-1)=b(i-1)/l(i-1);
l(i)=a(i)-c(i-1)*u(i-1);
y(i)=(d(i)-c(i-1)*y(i-1))/l(i);
end
% Sustitución regresiva
x(n)=y(n);
for i=n-1:-1:1
x(i)=y(i)-u(i)*x(i+1);
Método de disparo no lineal

Resolver el problema de contorno
y''  f (x , y, y')
y (a )  

y(b)  
iterando las soluciones de los PVI
y''  f (x , y, y')
y (a )  

x  a , b 
x  a , b 
y ' (a )  t k
eligiendo los parámetros t = tk para que
lim y ( b , t k )  y ( b )  
k
Disparo con la secante

tk: ángulo de tiro,
21
20
19
y'(a) = tk,
y(b,tk)
18
17
k=1,2,...
16
15
14
13

Iterar los tk según
la fórmula
t k 1
12
1
1.5
y(b, t k )  
 tk 
y ( b , t k )  y ( b , t k 1 )
t k  t k 1
2
2.5
3
3.5
4
4.5
5
5.5
6
k
Algoritmo de disparo secante


Entradas: f, a, b, , , tol, maxiter
Proceso
•
•
•
•

Estimar t0 y t1
Disparar con y(a) = , y’(a) = t0 para hallar y(b,t0)
Disparar con y(a) = , y’(a) = t1 para hallar y(b,t1)
Mientras |y(b,tk)  )| > tol y k < maxiter
Hallar tk+1 por la fórmula de la secante
Disparar con y’(a)=tk+1 para hallar y(b,tk+1)
Salida: y
Ejemplo
y''
1
8
32  2 x
y (1)  17
3
 yy '
y (3) 

x  1,3 
43
3

Solución exacta:
17
y(x )  x
2

16
16
x

15
Criterio de parada:
14
13
y ( 3, t k ) 
43
3
 10
5
12
11
1
1.2
1.4
1.6
1.8
2
2.2
2.4
2.6
2.8
3
Diferencias finitas: caso no lineal

Problema de contorno
y"  f ( x , y , y' ), x  a , b ;

y( a )   , y( b )  
Nodos de discretización
a = x0 < x1 < ... < xn < xn+1 = b

Aproximaciones en los nodos
y i  y ( xi )
y" i  y" ( x i ) 
y' i  y' ( x i ) 
y i 1  2 y i  y i 1
h
2
y i 1  y i 1
2h
Discretización del problema no
lineal
Para

Derivada
primera  zi
i  1 ,2 ,...., N
y i 1  2 y i  y i 1
h
2
y i 1  y i 1 

 f  xi , yi ,
0
2h


Derivada
segunda
y0  
y N 1  
Condiciones de contorno
Sistema no lineal
2 y 1  y 2  h f  x1 , y 1 , z 1     0


2
 y1  2 y 2  y 3  h f  x 2 , y 2 , z 2   0





2
 y N  2  2 y N 1  y N  h f  x N 1 , y N 1 , z N 1   0

2
 y N  1  2 y N  h f  x N , y N , z N     0 
2

Lo resolvemos por Newton
Jacobiano

Diagonal: para i = 1, 2, ..., n
a i  a i ,i  2  h f y ( x i , y i , z i )
2

Superdiagonal: para i = 1, 2, ..., n1
bi  a i ,i  1   1  ( h 2 ) f z ( x i , y i , z i )

Superdiagonal: para i = 1, 2, ..., n1
c i  a i  1 ,i   1  ( h 2 ) f z ( x i  1 , y i  1 , z i  1 )
Término independiente
    2 y 1  y 2  h 2 f  x1 , y 1 , z 1 

  y1  2 y 2  y 3  h 2 f  x 2 , y 2 , z 2 

d  




2
 y

2
y



h
f xN , yN ,zN
N 1
N

diff(y,2)







 
Ejemplo
y' ' 
1
8
32  2 x
3
 yy '

x  1 , 3 
function [f,fy,fz] = fun(x,y,z)
% y" = (32 + 2x^3 -yy')/8
% Valor de y"
f = (32 + 2*x.^3 -y.*z)/8;
% Parcial respecto a y
fy = -z/8;
% Parcial respecto a y'
fz = -y/8;
MétodosVariacionales:
Rayleigh-Ritz

TEOREMA: Bajo ciertas condiciones para las
funciones p(x), q(x) y f(x), y(x) es la solución
del problema de frontera

d 
dy 
p
(
x
)

  q( x) y  f ( x)
dx 
dx 
0  x 1
(1)
y ( 0 )  y (1)  0
si y sólo si y(x) es la única función que
minimiza la integral
I (u ) 
 p ( x ) u ' ( x ) 
1
0
2

 q ( x ) u ( x )   2 f ( x ) u ( x ) dx
2
(2)
Funciones base

La integral I se minimiza en el subespacio
generado por las funciones base
1 ,  2 ,...,  n 
Las funciones base
son linealmente independientes
verifican las condiciones de contorno


 i ( 0 )   i (1)  0
i  1, 2 ,..., n
Solución aproximada

Hallaremos una aproximación a la solución
y(x) de (1),
n
 ( x) 
 c
i
i
( x)
i 1


eligiendo los coeficientes para que
minimicen la integral I(f).
Se obtiene el sistema lineal A c = b donde
   p ( x )
1
a ij
bj 
0

1
0
'
i

( x ) ( x )  q ( x ) i ( x ) j ( x ) dx
f ( x ) j ( x ) dx
'
j
Funciones base: polinomios
lineales a trozos

Dada una partición de [0,1]
0  x 0  x 1  x 2  ...  x n  x n 1  1
donde
h i  x i 1  x i
0

x  x i 1

 h i 1
i ( x)  
 x i 1  x

hi

0
i  0 ,1,..., n
0  x  x i 1
x i 1  x  x i
i  1, 2 ,..., n
x i  x  x i 1
x i 1  x  1
Coeficientes del sistema
a ij 
  p ( x ) ( x )
a ii 
  dx    dx
1
0
xi
x i 1
'
i
x i 1
'
j

( x )  q ( x ) i ( x ) j ( x ) dx

xi
2
 1

( x  xi 1 )
  
p( x ) 
q ( x )  dx 
2
2
x i 1
hi  1
 hi  1

xi

x i 1
xi
2
 1

( xi 1  x )
q ( x )  dx
 2 p( x ) 
2
hi
 hi

Coeficientes del sistema
a i ,i  1 

x i 1
xi
 1

( x  x i )( x i  1  x )
q ( x )  dx
  2 p( x ) 
2
hi
 hi

Términos independientes
bi 

xi
( x  xi 1 )
x i 1
hi  1
f ( x ) dx 

x i 1
xi
( x  xi )
hi
f ( x ) dx
Algoritmo de elementos finitos
lineales


Entrada: problema, partición del intervalo
Proceso:
• En cada subintervalo [xi, xi+1] hallar las integrales
que aparecen el las fórmulas de los coeficientes
del sistema.
• Combinar adecuadamente las integrales
calculadas para obtener los coeficientes.
• Resolver el sistema lineal.

Salida: aproximación lineal a trozos de la solución.
Integrales a evaluar: i = 0, ..., n
ph i 
ql i 

x i 1
( x  xi )
xi
hi
rl i 

xi
hi
2
p ( x ) dx
( x  xi )
qr i 
q ( x ) dx
2
hi
xi
1
2
qh i 
x i 1

x i 1

x i 1

x i 1
xi
( x  x i )( x i  1  x )
xi
f ( x ) dx
hi
2
rr i 

hi
2
2
q ( x ) dx
q ( x ) dx
x i 1
xi
( xi 1  x )
( xi 1  x )
hi
f ( x ) dx
Coeficientes del sistema
Matriz del sistema

Diagonal: para i = 1, 2, ..., n
ai = phi1 + phi + qli1 + qri

Sub y superdiagonal: para i = 1, 2, ..., n1
bi = ci = qhi  phi
Términos independientes: para i = 1, 2, ..., n
di = rli1 + rri
 y ' ' p y  2p sen( p x )
2
Ejemplo


2
x  0 ,1, y ( 0 )  y (1)  0
Solución exacta: y(x)=sen (px)
Tomamos h=0.1, xi = 0.1 i
i=0,1,…,9,10
c 9  0 . 31028667
c 4  0 . 95496419
c 8  0 . 59020033
c 3  0 . 81234106
c 7  0 . 81234106
c 2  0 . 59020033
c 6  0 . 95496419
c1  0 . 31028667
c 5  1 . 00410877
FIN
Descargar

Problemas de frontera para ecuaciones diferenciales