PRACTICA6T_RES

Vista previa del material en texto

PRÁCTICA 6
semestre mayo-julio, 2017
1. Aplique el método de Taylor de orden tres a la ecuación y′ = −xy, con la condición inicial: y(0) = 1 para
estimar la solución de la ecuación diferencial alrededor de x = 0. Estimar la solución en x = 0.5
SOLUCIÓN:
y′ = −xy ⇒ y′(0, 1) = 0
y′′ = −y + xy′ = y(x− 1) ⇒ y′′(0, 1) = −1
y′′′ = y′(x− 1) + y = y(−x2 + x+ 1) ⇒ y′′′(0, 1) = 1
f(x) = f(x0) + f
′(x0)(x− 0) +
f ′′(x0)
2
(x− 0)2 + f
′′′(x0)
6
(x− 0)3 = 1− 1
2
x2 +
1
6
x3
f(0.5) = 1− 1
2
(0.5)2 +
1
6
(0.5)3 = 0.8958
2. Resuelva en forma anaĺıtica el problema de valores iniciales, en el intervalo [0, 1]. Grafique la solución.{
y′ = (1 + 2x)
√
y
y(0) = 1
donde y(0) = 1.
SOLUCIÓN:
y′ =
dy
dx
= (1 + 2x)
√
y ⇒ dy√
y
= (1 + 2x)dx ⇒ y−1/2dy = (1 + 2x)dx
−2√y = x+ x2 + C ⇒ −2
√
1 = C ⇒ C = ±2
y =
1
4
(x2 + x± 2)2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
1.5
2
2.5
3
3.5
4
 
Solución exacta
(a) Utilice el método de Euler con h = 0.5 y 0.25 para resolver el mismo problema. Grafique los resultados
en la misma gráfica para comparar en forma visual la exactitud de los dos tamaños de paso.
1
function f = fun1(t, y)
f = (1 + 2*t)*sqrt(y);
>> E = euler1(’fun1’, 0, 1, 1, 0.5)
E =
0 1.0000
0.5000 1.5000
1.0000 2.7247
>> E = euler1(’fun1’, 0, 1, 1, 0.25)
E =
0 1.0000
0.2500 1.2500
0.5000 1.6693
0.7500 2.3153
1.0000 3.2663
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
1.5
2
2.5
3
3.5
4
 
Solución exacta Solución exacta Euler con h = 0.5 Euler con h = 0.25
Solución utilizando método de Euler
(b) Emplee el método de Heun con h = 0.25 para resolver el mismo problema.
>> H = heun(’fun1’, 0, 1, 1, 0.5)
H =
0 1.0000
0.5000 1.2500
1.0000 1.8090
>> H = heun(’fun1’, 0, 1, 1, 0.25)
H =
0 1.0000
0.2500 1.1250
0.5000 1.3239
0.7500 1.6115
1.0000 2.0082
2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
1.5
2
2.5
3
3.5
4
 
Solución exacta Solución exacta Heun con h = 0.5 Heun con h = 0.25
Solución utilizando el método de Heun
(c) Emplee el método del punto medio con h = 0.25, para resolver el mismo problema.
>> R = puntomedio(’fun1’, 0, 1, 1, 0.5)
R =
0 1.0000
0.5000 1.8624
1.0000 3.8920
>> R = puntomedio(’fun1’, 0, 1, 1, 0.25)
R =
0 1.0000
0.2500 1.3346
0.5000 1.8836
0.7500 2.7277
1.0000 3.9710
3
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
1.5
2
2.5
3
3.5
4
 
Solución exacta Solución exacta Punto medio con h = 0.5 Punto medio con h = 0.25
Solución utilizando el método de punto medio
(d) Use el método de Runge-Kutta orden 2 y 4 con h = 0.25 para resolver el mismo problema.
>> R3 = rk3(’fun1’, 0, 1, 1, 0.25)
R3 =
0 1.0000
0.2500 1.3369
0.5000 1.8904
0.7500 2.7424
1.0000 3.9983
>> R4 = rk4(’fun1’, 0, 1, 1, 0.25)
R4 =
0 1.0000
0.2500 1.3369
0.5000 1.8906
0.7500 2.7431
1.0000 3.9998
4
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
1.5
2
2.5
3
3.5
4
 
Solución exacta Solución exacta RK−3 RK−4
Solución obtenida mediante RK-3 y RK-4
3. Un estanque se drena a través de un tubo como se observa en la figura. Con suposiciones simplificadoras, la
ecuación diferencial siguiente describe cómo cambia la profundidad con el tiempo:
dh
dt
= − πd
2
4A(h)
√
2g(h+ e)
donde h = profundidad (m), t = tiempo (s), d = diámetro del tubo (m), A(h) = área de la superficie del es-
tanque como función de la profundidad (m2), g = constante gravitacional (= 9.81 m/s2) y e = profundidad de
la salida del tubo por debajo del fondo del estanque (m). Con base en la tabla siguiente de área-profundidad,
resuelva esta ecuación diferencial para determinar cuánto tiempo tomaŕıa que el estanque se vaciara dado
que h(0) = 6 m, d = 0.25 m, e = 1 m.
h, m 6 5 4 3 2 1 0
A(h), ×104 m2 1.17 0.97 0.67 0.45 0.32 0.18 0
SOLUCIÓN:
f(x) = p1*x^5 + p2*x^4 + p3*x^3 + p4*x^2 + p5*x + p6
p1 = -10; p2 = 117.8; p3 = -388.6; p4 = 273.1; p5 = 1817; p6 = -1.299
dy/dt=-(pi*0.25^2)/(4*(-10*y^5+117.8*y^4-388.6*y^3+273.1*y^2+1817*y-1.299))*sqrt(2*9.81*(y+1))
donde con dy/dt hemos denotado dh/dt. Como el tiempo está dado en segundos aplicamos el método de
RK-4 con el siguiente comando:
rk4(’fun1’,0,68000,6,1)
donde fun1 es la función que almacena dy/dt. Se ha tomado el tiempo final de 68000 segundo porque la
función es racional y tiene discontinuidades en puntos cercanos al valor buscado. El primer valor del tiempo
5
obtenido para h = 0 se muestra en la parte de la tabla de los datos mostrada a continuación:
t y
6.7727e+ 04 1.3174e− 02
6.7728e+ 04 −8.9707e− 04
6.7729e+ 04 2.1752e− 02
Vemos que la solución comienza a fluctuar alrededor de 67727 y 67728 segundos. Graficando la solución
podemos observar este comportamiento.
0 1 2 3 4 5 6 7
x 10
4
−1
0
1
2
3
4
5
6
Por lo tanto, el estanque se vaciará aproximadamente después de 67727 segundos que equivale aproximada-
mente a 18.8 horas.
4. La ecuación diferencial básica de la curva elástica para una viga en voladizo (véase la figura) está dada por
EI
d2y
dx2
= −P (L− x)
donde E = módulo de elasticidad e I = momento de inercia. Resuelva para la deflexión de la viga con el
empleo de los métodos: (a) RK-4, (b) diferencias finitas centradas. Se aplican los valores siguientes de los
parámetros: E = 30 000 ksi, I = 800 in4, P = 1 kip, L = 10 ft. Compare sus resultados numéricos con la
solución anaĺıtica,
y = −PLx
2
2EI
+
Px3
6EI
Este problema puede resolverse de varias maneras. El más cómodo es convertir la EDO de segundo orden en
una EDO de primer orden. La ecuación permita hacer esto.
d2y
dx2
= − P
EI
(L− x) ⇒ dy
dx
= − P
EI
(
Lx− x
2
2
)
+ C
6
Según el problema las CI son y(0) = 0 y y′(0) = 0 (viga en voladizo). Por lo tanto, C = 0 y la ecuación por
resolver queda:
dy
dx
= − P
EI
(
Lx− x
2
2
)
a) RK-4
Conversión de unidades: P = 1 kip = 4448.22 N; L = 10 ft = 3.05 m; E = 30 000 ksi = 2.07× 1011 Pa;
I = 800 in4 = 6.48× 10−4 m4.
function f = fun2(x, ~)
E = 2.07e11; I = 6.48e-4; P = 4448.22; L = 3.05; f = -(P/(E*I))*(L*x - x^2/2);
rk4(’fun2’,0,3.05,0,0.305)
ans =
0 0
3.0500e-01 -4.5477e-06
6.1000e-01 -1.7563e-05
9.1500e-01 -3.8106e-05
1.2200e+00 -6.5235e-05
1.5250e+00 -9.8010e-05
1.8300e+00 -1.3549e-04
2.1350e+00 -1.7673e-04
2.4400e+00 -2.2080e-04
2.7450e+00 -2.6674e-04
3.0500e+00 -3.1363e-04
Por lo tanto, la deflexión en el extremo L = 3.05 m es aproximadamente 0.31363 mm.
0 0.5 1 1.5 2 2.5 3 3.5
−3
−2
−1
0
x 10
−4
x(m)
y
(m
)
RK−4 Solución exacta
b) Diferencias finitas
yi+1 = yi −Ah(Lxj − (1/2)x2j );
P = problema(0, 3.05, 0, 0.305)
P =
0 0
3.0500e-01 0
6.1000e-01 -8.9385e-06
9.1500e-01 -2.5875e-05
7
1.2200e+00 -4.9867e-05
1.5250e+00 -7.9976e-05
1.8300e+00 -1.1526e-04
2.1350e+00 -1.5478e-04
2.4400e+00 -1.9759e-04
2.7450e+00 -2.4275e-04
3.0500e+00 -2.8932e-04
0 0.5 1 1.5 2 2.5 3 3.5
−3
−2
−1
0
x 10
−4
x(m)
y
(m
)
Diferencias finitas, h = 0.305 Solución exacta
>> P = problema(0, 3.05, 0, 0.0305)
P =
0 0
3.0500e-02 0
6.1000e-02 -9.3619e-08
9.1500e-02 -2.7992e-07
1.2200e-01 -5.5795e-07
......................
2.9890e+00 -3.0186e-04
3.0195e+00 -3.0657e-04
3.0500e+00 -3.1127e-04
8
0 0.5 1 1.5 2 2.5 3 3.5
−3
−2
−1
0
x 10
−4
x(m)
y
(m
)
Diferencias finitas, h=0.0305
Solución exacta
El valor aproximado de la deflexión en el primer caso es 0.28932 mm y en el segundo caso 0.31127
mm. Notamos que se necesita un h mucho menor para que se obtiene un resultado similar del resultado
obtenido con RK-4.
5. La ecuación diferencial básica de la curva elástica para una viga con carga uniforme (véase la figura) está
dada por
EI
d2y
dx2
=
wLx
2
− wx
2
2
donde E = módulo de elasticidad e I = momento de inercia. Resuelva para la deflexión de la viga con los
métodos de (a) diferencias finitas (h = 1 ft), y b) disparo. Aplique los siguientes valores de parámetros: E =
30 000 ksi, I = 800 in4, w = 1 kip/in, L = 10 ft. Compare sus resultados numéricos con la solución anaĺıtica,
y =
wLx3
12EI
− wx
4
24EI
− wL3x
24EI
SOLUCIÓN:
a) Diferencias finitas
function f = fun1(x, y)
w = 15000; E = 2.07e+11; I = 6.48e-4; L = 3.00; f = w*(L*x - x.^2)/(2*E*I);
function y = p(t)
y = 0;
function y = q(t)
y = 0;
9
function y = r(t)
E = 2.07e+11; I = 6.48e-4; w = 15000; L = 3.00; y = w*(L*x - x^2)/(2*E*I);
function DDF = diffin(fun, a, b, n, alfa, beta)
h = (b - a)/(n + 1); T = a:h:b;
p = zeros(1,n); q = zeros(1,n); r = feval(fun,T(2:n+1));
Aa = 2 + (h^2)*q; Ab = -1 + h*p/2; Ac = -1 - h*p/2;
A = diag(Aa) + diag(Ab(1:n-1),1) + diag(Ac(2:n),-1);
d = -(h^2)*r; d(1) = d(1) - alfa*Ac(1); d(n) = d(n) - beta*Ab(n); d = d(:);
Y = A\d; Y = [alfa, Y’, beta];
DDF = [T’ Y’];
% PLOT
E = 2.07e+11; I = 6.48e-4; w = 15000; L = 3.00;
syms x
yexact=((w*L*x^3)/(12*E*I))-(((w*x^4)/(24*E*I)))-(((w*L^3*x)/(24*E*I)));
xexact = 0:0.3:L;
m = subs(yexact,x,xexact);
plot(DDF(:,1),DDF(:,2),’r’)
hold on, plot(xexact,m,’b’),
hold off
grid, xlabel(’Longitud x, (m)’), ylabel(’Deflexión y, (m)’)
>> diffin(’fun1’,0,3,9,0,0)
ans =
0 0
3.0000e-01 -3.7364e-05
6.0000e-01 -7.0652e-05
9.0000e-01 -9.6694e-05
1.2000e+00 -1.1322e-04
1.5000e+00 -1.1889e-04
1.8000e+00 -1.1322e-04
2.1000e+00 -9.6694e-05
2.4000e+00 -7.0652e-05
2.7000e+00 -3.7364e-05
3.0000e+00 0
0 0.5 1 1.5 2 2.5 3
−1.2
−1
−0.8
−0.6
−0.4
−0.2
0
x 10
−4
Longitud x, (m)
D
e
fl
e
x
ió
n
 y
, 
(m
)
Diferencias finitas Solución exacta
b) Disparo
10
Datos: E = 2.07e+ 11; I = 6.48e− 4; w = 15000; L = 3.00; A = w/(24EI).
Problema 1: Problema 2:
y′′ = A(Lx− x2) y′′ = 0
y(0) = 0, y′(0) = 0 y(0) = 0, y′(0) = 1
Sistema 1 SIstema 2:
y′ = u1(x) y
′ = v1(x)
u′1 = A(Lx− x2) v′1 = 0
y(0) = 0, u1(0) = 0 y(0) = 0, v1(0) = 1
Ecuaciones del metodo de Euler
u1,i+1 = u1,i + hu2,i v1,i+1 = v1,i + hv2,i
u2,i+1 = u2,i + hA(Lxi − x2i ) v2,i+1 = v2,i
x0 = 0 u1,0 = 0, u2,0 = 0 x0 = 0 v1,0 = 0, v2,0 = 1
Finalmente, m = −u1,b
v1,b
y yi = u1,i +m · v1,i. Los resultados se muestran en la tabla siguiente:
x u1 u2 v1 v2 m yapr yex
0 0 0 0.3 1 -0.000113225 0 0
0.3 0 1.3587E-05 0.6 1 -0.000113225 -6.79348E-05 -3.7024e-05
0.6 4.07609E-06 3.77415E-05 0.9 1 -0.000113225 -9.78261E-05 -7.0048e-05
0.9 1.53986E-05 6.94444E-05 1.2 1 -0.000113225 -0.000120471 -9.5901e-05
1.2 3.62319E-05 0.000105676 1.5 1 -0.000113225 -0.000133605 -1.1232e-04
1.5 6.79348E-05 0.000143418 1.8 1 -0.000113225 -0.00013587 -1.1794e-04
1.8 0.00011096 0.00017965 2.1 1 -0.000113225 -0.000126812 -1.1232e-04
2.1 0.000164855 0.000211353 2.4 1 -0.000113225 -0.000106884 -9.5901e-05
2.4 0.000228261 0.000235507 2.7 1 -0.000113225 -7.74457E-05 -7.0048e-05
2.7 0.000298913 0.000249094 3.0 1 -0.000113225 -4.07609E-05 -3.7024e-05
3.0 0.000373641 0.000249094 3.3 1 -0.000113225 0 0
0 0.5 1 1.5 2 2.5 3
−1.4
−1.2
−1
−0.8
−0.6
−0.4
−0.2
0
x 10
−4
Longitud x, m
D
e
fl
e
x
ió
n
 y
, 
m
 
Solución exacta Solución aproximada
11