Instituto Politécnico Nacional
Unidad Profesional Interdisciplinaria de Biotecnología Unidad de Aprendizaje. Métodos Numéricos
Tarea V
Interpolación de Lagrange Diferenciación Numrica
Profesora. Carmen Marín Albino Profesor. Jorge Luis Rosas Mendoa
Alumnos.
Abad Ju!re "scar Alfredo Alfredo Barr#n Rodrígue $afne %ictoria $orantes Ca&as Isaac Alí 'ern!nde 'ern!nde Marisol Moreno 'errera Ale(andro Noriega )aspar *dgar Leonardo
Grupo +MM,
No-iembre ,+ de ./,0
*1uipo 2
Pro!lema ". Considere la siguiente tabla para -apor de agua calentada a .// MPa3 3
a4 $etermine la entropía 5 para un -olumen específico - de /6,/2
m usando una kg
interpolaci#n de La)range cuadr!tica6 b4 *ncuentre el -olumen correspondiente a una entropía de 060
kJ usando una kg∙K
interpolaci#n de La)range in-ersa6 Tabla 1. Valores de vapor de agua calentada a 200 MPa.
( )
/6,/788
/6,,,++
/6,.9+/
( )
06+,+8
069+97
06800+
3
m # kg kJ $ kg∙K
$olución. *s necesario comprender 1ue el problema no presenta ma:or complicaci#n;
pues lo 1ue se est! solicitando es simplemente el c!lculo de -alores representati-os de una interpolaci#n de La)range atendiendo a los puntos 1ue se presentan para proponer una cur-a cuadr!tica6 Así; pues; se procede a realiar una gr!fica de los datos en Matlab de modo 1ue se pueda analiar el comportamiento de los mismos; la gr!fica
>>v= [0.10377 0.1144 0.12540]; >> S= [6.4147 6.5453 6.7664];
>> plot(v, S,'ro','markerfacecolor',''!,"r#$ o%, &lael('olme% e)pec*f#co (k"+m!'!,-lael('%trop*a (k/+k"!'!,le"e%$('ato) prope)to)'!,t#tle('%trop*a e% vapor $e a"a cale%ta$a a 200 a'!
Lo anterior puede significar 1ue esos puntos son factibles de encontrarse en una par!bola; de modo 1ue se puede seguir el método de interpolaci#n de La)range para encontrar un polinomio de interpolaci#n de segundo grado6 Así es como se plantea 1ue3 P2 ( x )= S 1 L1+ S 2 L2 + S3 L3
Para lo cual deben plantearse las ecuaciones 1ue describir!n a cada una de las rectas3 L1=
(v −v 2 )( v − v 3 ) ( v 1−v 2 )( v 1− v 3)
L2 =
(v −v 1 )( v −v 3 ) ( v 2−v 1 )( v 2− v 3)
L3=
( v −v 1 )( v −v 2 ) ( v 3−v 1 )( v 3 −v 2 )
Luego entonces; el c#digo en Matlab 1ue realiar! lo anteriormente descrito ser! el 1ue se muestra a continuaci#n; adem!s; los resultados se almacenan en un -ector para simplificar el c!lculo del polinomio de segundo grado3 )-m) (1!=e&pa%$((v(2!!(v(3!!+((v(1!v(2!!(v(1!v(3!!!!; (2!=e&pa%$((v(1!!(v(3!!+((v(2!v(1!!(v(2!v(3!!!!; (3!=e&pa%$((v(1!!(v(2!!+((v(3!v(1!!(v(3!v(2!!!!;
=bteniendo los siguientes -alores para cada una de las rectas; los -alores 1ue se presentan a continuaci#n son los obtenidos luego de aplicar un -alor de precisi#n aritmética de ocnicamente grandes n>meros debido a 1ue % se traba(# como una -ariable simb#lica3 L(1)= 4349.2084*V^2-1042.9402*V+62.3927 L(2)= -8552.1252*V^2+1959.8905*V-111.28674 L(3)= 4202.9168*V^2-916.95036*V+49.894036
?a 1ue los datos se encuentran organiados seg>n la correspondencia entre los puntos; es mu: f!cil realiar la determinaci#n del polinomio de segundo grado; puesto 1ue; si se realia un producto entre dos -ectores; se multiplicar! elemento a elemento : se obtendr! como resultado la suma de cada uno de esos productos de al siguiente manera3 2=)m(S.! P2=361.2583*V^2 - 66.529739*V + 9.4283848
Así es como el polinomio anterior ser! el 1ue defina a la par!bola 1ue tomar! a los tres puntos dados para la realiaci#n de la interpolaci#n de La)range; como se muestra su correcto a(uste en la siguiente gr!fica; cu:o c#digo de Matlab es3 >>ol$ o% >> e8plot(2,[m#%(v!,ma&(v!],'"'!,t#tle('alore) $e a"a cale%ta$a a 200 a'!,&lael('olme% e)pec*f#co (m+k"!'!,-lael('%trop*a
k/+k"!'!,le"e%$('ato) (<2(&!!'!
e&per#me%tale)','9%terpolac#:%
$e
ara%"e
Para culminar con este e(ercicio es necesario resol-er los incisos a : b; respecti-amente; lo cual es mu: simple puesto 1ue :a se cuenta con el polinomio de segundo grado; así 1ue3 3
a% Con un -olumen específico de /6,/2
m sobre el polinomio; se aplicar! el kg
siguiente c#digo de Matlab para la determinaci#n de la entropía3 5,@subsP.;/6,/24 $"&'.()'* !% Con una entropía de 060
kJ kg∙K
kJ kg∙K
se determinar! un -olumen específico por
medio de una interpolaci#n in-ersa; considerando 1ue se debe igualar tal ecuaci#n a cero; es decir3 V1=solve(361.2583*V^2 - 66.529739*V + 9.4283848-6.6) 3
V1=
0.11756727
m kg
3
0.066593865
m kg
3
*s e-idente 1ue el -alor 1ue debe ser considerado es el de /6,,8908.8
m puesto kg
1ue la par!bola en tal punto se encuentra en un estado ascendiente6 Pro!lema +. *n la siguiente tabla se muestra el -olumen específico de -apor s>per
caliente a -arias temperaturas a una presi#n absoluta de T ,-% # ,Lt+ 34g%
3,000
lb
¿2
6
+/0
+*1
+2(
(0'
("*
967,7
869272
262+.2
680/
,/697,,
$eterminar el -olumen específico a una temperatura de 89/ D6 5esolución
La f#rmula para con-ertir grados Da
° F −32 1.8
Para la determinaci#n del -olumen específico se con-ierten los grados Da
−32
750
1.8
=398.88 ° C
*s entonces 1ue; se realia el c#digo para la obtenci#n del Polinomio 1ue representa los puntos dados : posteriormente conocer el -olumen específico de -apor s>per caliente a 72622C6 clc; clear all; clo)e all; Se crea% lo) vectore) $e lo) $ato) $a$o) =[370 3?2 3@4 406 41?]; v=[5.@313 7.5?3? ?.?42? @.7@60 10.5311]; Se "raf#ca% lo) p%to) para ote%er % o)AeBo $e la relac#:% e%tre lo) $ato) plot(,v, 'o', 'markerfacecolor', 'r'! &lael('emperatra'!; -lael('volme% e)pec*f#co'!; t#tle('apor )Cper cal#e%te'!;
Se $eclara %a var#ale )#m:l#ca para po$er real#8ar lo) cDlclo) $e lo) a"ra%"#a%o) o) la"ra%"#a%o) )o% $ato) #%terpola$o) para tomar lo) meBore) co%B%to) $e $ato) - e%co%trar el mo$elo Ae meBor $e)cr#a ) comportam#e%to )-m) t (1!=e&pa%$((t(2!!(t(3!!(t(4!!(t(5!!!+(((1!(2!!((1! (3!!((1!(4!!((1!(5!!!; (2!=e&pa%$((t(1!!(t(3!!(t(4!!(t(5!!!+(((2!(1!!((2! (3!!((2!(4!!((2!(5!!!; (3!=e&pa%$((t(1!!(t(2!!(t(4!!(t(5!!!+(((3!(1!!((3! (2!!((3!(4!!((3!(5!!!; (4!=e&pa%$((t(1!!(t(2!!(t(3!!(t(5!!!+(((4!(1!!((4! (2!!((4!(3!!((4!(5!!!; (5!=e&pa%$((t(1!!(t(2!!(t(3!!(t(4!!!+(((5!(1!!((5! (2!!((5!(3!!((5!(4!!!; Se crea el ol#%om#o corre)po%$#e%te co% lo) a"ra%"#a%o) =vpa(v(1!(1!Ev(2!(2!Ev(3!(3!Ev(4!(4!Ev(5!(5!,4! ol$ o%; −22
P=7.025 x 10
4
−6
3
−3
2
t + 8.459 x 10 t − 11.06 x 10 t + 4.867 t −709.2
Se "raf#ca el mo$elo ote%#$o por #%terpolac#:% $e a"ra%"e Fae $e)tacar Ae el mo$elo ote%#$o e) m- parec#$o al comportam#e%to $e lo) $ato) rto) "raf#ca$o) a%ter#orme%te e8plot(,[m#%(!],[ma&(!]!; &lael('emperatra'!; -lael('volme% e)pec*f#co'!; t#tle('ol#%om#o $el vapor )Cper cal#e%te'!; "r#$ o%;
or Clt#mo, )e calcla el volme% e)pec*f#co $el vapor )Cper cal#e%te a 3@?.?? "ra$o) Fe%t*"ra$o), o)erva%$o la "rDf#ca )e e)t#ma Ae )erD % valor aproma$o a @.3 ve=))(,3@?.??!
*l -alor obtenido es de
ve = 9.2860
3
¿
: es un -alor mu: parecido al estimado; por lo
Kg
1ue; obser-ando la gr!fica; es un -alor 1ue corresponde al modelo obtenido por interpolaci#n de Lagrange6 Pro!lema '. Un a-i#n est! siendo rastreado por un radar : los datos se toman cada .
segundos en coordenadas polares θ , r 3 t (s)
200
202
204
206
208
210
θ
0.75
0.72
0.68
0.68
0.67
0.66
5120
5370
5560
5800
6030
6240
(rad) r (m)
*stime la -elocidad : aceleraci#n -ectoriales dadas por3 v =´r er + r θ ´eθ ⃗
⃗
⃗
a =( ´r − r θ´ ) er + ( r θ´ + 2 r´ ´θ ) eθ 2
⃗
⃗
⃗
$onde los puntos sobre r : θ significan el orden de la deri-ada; : son los -ectores unitarios6
er ⃗
:
eθ ⃗
5esolución
5e conoce el tama&o de paso al estar dado por el tiempo en 1ue se tomaron las coordenadas del a-i#n; por lo 1ue <@.; sin embargo; en el problema no se especifica x i : es entonces 1ue se decide realiar la deri-ada de la ra#n entre las coordenadas
por medio de diferencias entre los puntos proporcionados6 Al tener un punto antes de r,4; se comiena a obtener diferencias desde el siguiente punto r.4; perdiéndose entonces un punto por cada diferencia deri-ada4 realiada6 *s entonces predecible 1ue para el c!lculo de la segunda deri-ada por diferencia de puntos; el -ector resultante ser! de + elementos; siendo éste >ltimo -ector el indicador de 1ue todos los -ectores tendr!n 1ue ser de + elementos para la realiaci#n de productos entre -ectores para el c!lculo de las -elocidades : aceleraciones para cada tiempo6 5e deciden tomar los >ltimos + -alores para cada uno de los -ectores conocidos : calculados6 5iguiendo el planteamiento; se tiene el siguiente c#digo para resoluci#n al problema3 clc; clear all; clo)e all; Se crea% lo) vectore) co% lo) $ato) proporc#o%a$o) t=[200 202 204 206 20? 210]; teta=[0.75 0.72 0.70 0.6? 0.67 0.66]; r=[5120 5370 5560 5?00 6030 6240]; Se $eclara% lo) vectore) %#tar#o), el mot#vo $e lle%arlo co% 4 eleme%to) e) Ae al $er#var e% or$e% 1 - 2 la) formla) proporc#o%a$o), )e p#er$e % $ato por or$e%, e) e%to%ce) Ae al f#%al )e ce%ta co% % vector $e 4 - para real#8ar el pro$cto e%tre vectore) $ee% co%tar co% la m#)ma ca%t#$a$ $e eleme%to) er=[1 1 1 1]; eteta=[1 1 1 1]; Gl %o )aer la $#)ta%c#a $el p%to a e)t#mar, )e real#8a la $er#va$a co% %a relac#:% e%tre lo) p%to) $a$o) re)pecto al t#empo proporc#o%a$o, e) e%to%ce) Ae )e ot#e%e% $o) vectore), %o repre)e%ta%$o al ra$#o - el otro repre)e%ta%$o a teta $r(1!=(r(2!r(1!!+(t(2!t(1!!; $r(2!=(r(3!r(2!!+(t(3!t(2!!; $r(3!=(r(4!r(3!!+(t(4!t(3!!; $r(4!=(r(5!r(4!!+(t(5!t(4!!; $r(5!=(r(6!r(5!!+(t(6!t(5!!; $teta(1!=(teta(2!teta(1!!+(t(2!t(1!!; $teta(2!=(teta(3!teta(2!!+(t(3!t(2!!; $teta(3!=(teta(4!teta(3!!+(t(4!t(3!!; $teta(4!=(teta(5!teta(4!!+(t(5!t(4!!; $teta(5!=(teta(6!teta(5!!+(t(6!t(5!!; ara el cDlclo $e la )e"%$a $er#va$a, )e p#er$e % $ato mD), e) e%to%ce) Ae )e ce%ta co% 4 $ato) f#%ale) para po$er co%ocer la veloc#$a$ - la acelerac#:% para el av#:% e)t$#a$o $$r(1!=($r(2!$r(1!!+(t(2!t(1!!; $$r(2!=($r(3!$r(2!!+(t(3!t(2!!; $$r(3!=($r(4!$r(3!!+(t(4!t(3!!; $$r(4!=($r(5!$r(4!!+(t(5!t(4!!; $$teta(1!=($teta(2!$teta(1!!+(t(2!t(1!!; $$teta(2!=($teta(3!$teta(2!!+(t(3!t(2!!; $$teta(3!=($teta(4!$teta(3!!+(t(4!t(3!!; $$teta(4!=($teta(5!$teta(4!!+(t(5!t(4!!; Gl )er po)#le ):lo % pro$cto e%tre vectore) co% la m#)ma ca%t#$a$ $e eleme%to), )e toma% lo) Clt#mo) $ato) ote%#$o) ta%to $e lo) ra$#o) -a $a$o) como $e la) $er#va$a) $e or$e% 1 - 2 para r - teta $ator(1!=r(3!; $ator(2!=r(4!; $ator(3!=r(5!; $ator(4!=r(6!; $ato$r(1!=$r(2!;
$ato$r(2!=$r(3!; $ato$r(3!=$r(4!; $ato$r(4!=$r(5!; $ato$teta(1!=$teta(2!; $ato$teta(2!=$teta(3!; $ato$teta(3!=$teta(4!; $ato$teta(4!=$teta(5!; ) f#%alme%te Ae, )e ))t#t-e% lo) $ato) ote%#$o) e% la) ecac#o%e) $a$a) - )e co%oce la ma"%#t$ $e la veloc#$a$ - la acelerac#:% para lo) Clt#mo) 4 t#empo) $a$o) v=$ato$r.erE$ato$teta.eteta.$ator a=($$r$ator.$ato$r.($ato$teta.H2!!.erE ($ator.$$tetaE2.$ato$r.$$teta!.eteta
v =39.400062.0000 84.8500 73.8000 a =−53.4450 −57.1000− 4.1863−21.3800
Dinalmente; los -ectores obtenidos indican 1ue; para cada tiempo se tienen diferentes -elocidades : aceleraciones; mostr!ndose los resultados en la siguiente tabla3 Tiempo ,s% Velocidad ,m3s% Aceleración ,m3s 1%
10(
10'
10*
1"0
76+/// E976++9/
0.6//// E986,///
2+629// E+6,207
8762/// E.,672//
Pro!lema /. *l mercurio se dilata con la temperatura de forma 1ue ,7969/2 g de éste
elemento ocupa a diferentes temperaturas; los -ol>menes 1ue se indican en la siguiente tabla3 F C4
,0/
,2/
.//
../
.+/
.0/
.2/
7//
% cmG4
,/6.778 .
,/677/08 ,
,/6702,8
,/6+/9. 9
,/6++77 0
,/6+2..+ 9
,/69./22 +
,/69922 2
5abiendo 1ue el coeficiente de eHpansi#n c>bica α e n K − 1 ¿ se define como3 α =
( )
dV V dT 1
$etermine el coeficiente de eHpansi#n c>bico en cada punto : encuentre la cur-a 1ue me(or a(uste entre la temperatura : el coeficiente de eHpansi#n c>bica; con ella; estime el -alor del coeficiente de eHpansi#n c>bica para el mercurio a una temperatura de ,//C6 $olución. *s necesario obser-ar 1ue el problema busca encontrar la primera deri-ada del
-olumen con respecto de la temperatura; de modo 1ue lo necesario es realiar un an!lisis numérico con diferentes par!metros; pero tomando en cuenta 1ue la retícula es uniforme; por lo 1ue la resoluci#n del problema se simplifica completamente6 $e los datos se pueden rescatar los siguientes par!metros importantes3
*l anc
*n el primer punto se realia una diferencia finita
Así pues; se determina la primera deri-ada para cada uno de los puntos; atendiendo al
Luego entonces; se procede a realiar el método de diferencias finitas
V ( T i )=
V ( T i + 1 )−V ( T i ) h
ue ser! guardado en un -ector para an!lisis posterior de cada una de las deri-adas como3 dV(1)=(V(2)-V(1))/h;
dV(1)=0.00186495
Posteriormente; se realia un ciclo for para el an!lisis de cada una de las diferencias finitas centradas desde el segundo
V ( T i )=
V ( T i + 1 )−V ( T i−1 ) 2h
ue ser! guardado en un -ector para an!lisis posterior de cada una de las deri-adas como3 for #=2I7 $(#!=((#E1!(#1!!+(2! e%$
Por >ltimo; la >ltima de las deri-adas se dar! con el método de diferencias finitas
V ( T i )=
V ( T i ) −V (T i−1) h
ue ser! guardado en un -ector para an!lisis posterior de cada una de las deri-adas como3 dV(8)=(V(8)-V(7))/h;
dV(8)=0.0019502
Luego entonces; se tendr! un -ector 1ue almacenar! cada una de las deri-adas de primer grado obtenidas; el cual se muestra a continuaci#n3 T,- %
,0/ ,2/ .// ../ .+/ .0/ .2/ 7//
dV3dT
/6//,20+9 / /6//,28/,8 9 /6//,22,79 / /6//,27. 9 /6//,/2// / /6//,.78/ / /6//,+,/8 9 /6//,9/./ /
Una -e 1ue se tiene esta deri-ada; se puede proceder a calcular el coeficiente de eHpansi#n c>bico en cada punto; simplemente con el c#digo3 l!=(1./V).*dV
*sto dar! como resultado un nue-o -ector alfa; el cual albergar! los -alores de cada uno de los -ol>menes; el cual se muestra a continuaci#n3 T,- %
,0/ ,2/ .// ../ .+/ .0/ .2/ 7//
−1
α en 1.0 e − 03 ( K
)
/6,2,809/9/787 /6,2,/7,7,.0.2/9 /6,2,+9+.70/2,782 /6,2.//++27/..20 /6,2.028++990/88 /6,279,2+8,7,27 /6,2++87.7,,28. /6,2+082/7/,22
*l problema no termina a>n; pues si bien :a se tiene el coeficiente para cada uno de los -ol>menes dados originalmente; es necesario encontrar la relaci#n 1ue guarda la temperatura : el coeficiente de eHpansi#n c>bica; una primera aproHimaci#n a esta relaci#n se mostrar! en el momento en 1ue se grafi1uen ambos par!metros en un plano; lo 1ue se lle-ar! a cabo mediante el c#digo3 >> plot(,alfa,'ro','markerfacecolor',''!,"r#$ o%, t#tle('Kelac#:% e%tre temperatra - coef#c#e%te $e e&pa%)#:% cC#ca'!,&lael('emperatra (LF!'!,-lael('Foef#c#e%te $e e&pa%)#:% cC#ca (H1!'!
La gr!fica 1ue se muestra a continuaci#n es la resultante del c#digo :a redactado; una cuesti#n interesante en el an!lisis de esta gr!fica es la forma 1ue guarda; puesto 1ue es mu: similar a una sigmoidea; lo 1ue puede orillar a un comportamiento de una funci#n c>bica3
Para poder obtener la me(or recta de a(uste es necesario plantear un modelo matem!tico de grado tres 1ue se pueda a(ustar a los puntos en tal rango; por lo 1ue si3 y = a x
3
2
+ b x + x + d
*s el modelo planteado; deber! generarse una matri de +H+ para poder
∑ x ∑ x
N
∑ x ∑ x ∑ x
2
3
² 3 x 4 x
∑ ∑
∑ x ² ∑ x ³ ∑ x ∑ x 4
5
∑ x ³ ∑ x ∑ x ∑ x 4
5
6
d c
=
b a
∑ y ∑ yx ∑ yx ∑ yx
2
3
La matri tendr! la siguiente forma; con los c#digos de Matlab3 length(T)
sum(T)
sum(T.ˆ2)
sum(T.ˆ3)
d
sum(T)
sum(T.ˆ2)
sum(T.ˆ3)
sum(T.ˆ4)
c
sum(T.ˆ2)
sum(T.ˆ3)
sum(T.ˆ4)
sum(T.ˆ5)
b
sum(T.ˆ3)
sum(T.ˆ4)
sum(T.ˆ5)
sum(T.ˆ6)
a
=
sum( α ) sum( α .*T) sum( α .*T.ˆ2) sum( α .*T.ˆ3)
Así es como la matri tendr! los siguientes -alores numéricos en formato corto3 8 1840 440000 1082800 0
1840 440000 10828000
440000 10828000 108280000 2.!!82e"10 2.!!82e"10 !.2646e"12
d c b
0.0015 0.3365 80.5855
2.!!82e"10
!.2646e"12
a
1.!5e"04
1.31e"15
*ntonces; resol-iendo por comandos la matri; se tendr!3 M=#%v(G! " = 1.0e-03 * 0.215499298335864 -0.000483938465551 0.000002156270690 -0.000000002949196
Por lo 1ue el modelo de a(uste para la gr!fica 1ue se obser-# con anterioridad ser!3 3
2
y =0.0002154 x − 0.0000004839 x + 0.0000000021 56 x −0. 000000000002949196
Así es como se obser-a el modelo a(ustado con los puntos eHperimentales descritos anteriormente3 >> plot(,pol-val(3,!,'r'!,le"e%$('ato)N,'ol#%om#o aB)ta$o'!
Por >ltimo; se deber! definir el coeficiente de eHpansi#n c>bica del mercurio; considerando una temperatura de ,//C; recuérdese 1ue la funci#n 1ue me(or a(ust# los puntos eHperimentales fue3 3
2
α =0.0002154 T −0.0000004839 T + 0.000000002 156 T − 0.00000000000294919
Por lo 1ue; simplemente sustitu:endo en esa funci#n el -alor de la temperatura se podr! obtener un estimado lo suficientemente aceptable para el coeficiente de eHpansi#n c>bica3 α =0.0002 ( 10 0 )−0.00000048 ( 100 )+ 0.000000002 ( 100 )−0. 000000000002949 3
2
−1
α =1.857189622407978 e − 04 K
5er! el coeficiente de eHpansi#n c>bica a los ,//C en el mercurio6