1
Solución de la Ecuación del Calor en Una Dimensión en Estado Estacionario Mediante Diferencias Finitas F. I. Saldaña, W. J. Dután y W. N. Bernal, Estudiantes de Maestría, Universidad Estatal de Cuenca
Resumen—Se presenta el modelo de la ecuación de transferencia de calor en una dimensión y su resolución mediante diferencias finitas. El problema se analiza en estado estable y sujeto a condiciones de frontera mixtas, de Dirichlet y Neumann. Para completar el trabajo, se desarrollan los códigos en MATLAB para el esquema planteado y se presentan ejemplos de aplicación. Los códigos de MATLAB son claros y permiten al lector apreciar la implementación de los mismos. Palabras Clave—condiciones de frontera, EDP, ecuación del calor, diferencias finitas, MATLAB®.
I. NOMENCLATURA ( , ) es la derivada parcial de ( , ) respecto al tiempo. ( , ) es la segunda derivada parcial de ( , ) respecto a la variable x (una dimensión). II. INTRODUCCIÓN
E
L problema abordado en este documento es la resolución numérica de la ecuación del calor unidimensional en estado estable mediante diferencias finitas. La ecuación del calor es un modelo matemático (quizás el más sencillo) que trata de describir la evolución de la temperatura en un cuerpo sólido. Se establecen condiciones de frontera mixtas: de Dirichlet y de Neumann, y la presencia de una fuente de calor. En el modelado de sistemas mediante ecuaciones en derivadas parciales en muy pocas ocasiones pueden obtenerse soluciones analíticas a estos problemas. Aunque existen varios enfoques distintos para la resolución numérica de estos problemas, quizá el más sencillo y extendido en el medio ingenieril es el método de las diferencias finitas. Todos los problemas en derivadas parciales envuelven condiciones de frontera donde el valor de u o alguna derivada parcial de u es especificado en la frontera del dominio.
III. SOLUCIÓN DE LA ECUACIÓN DE CALOR Muchos problemas de valores de contorno surgen de soluciones en estado estable (estacionario) de problemas transitorios. En el presente caso, se modela el problema estacionario de conducción del calor por difusión en presencia de una fuente y sometido a condiciones de frontera, y se plantea la solución del mismo mediante el método de diferencias finitas. El método de diferencias finitas es una de varias técnicas para obtener la solución de la ecuación en derivadas parciales. A. La Ecuación de Calor Se considera el flujo de calor en una barra delgada constituida de algún material conductor del calor, sometida a una fuente de calor externa a lo largo de su longitud y condiciones de frontera en cada extremo. Si asumimos que las propiedades del material, la distribución inicial de temperatura, y la fuente varían únicamente con x, la distancia a lo largo de la barra, y una sección transversal uniforme de la barra, entonces la distribución de temperatura en cualquier tiempo variará únicamente con x y podemos modelar este sistema físico con una ecuación en derivadas parciales en una dimensión espacial. Puesto que la solución podría variar con el tiempo, u(x,t) denota la temperatura en el punto x al tiempo t, donde a
( , ) = ( κ ( ) ( , )) + ( , )
(1)
donde κ ( ) es el coeficiente de conducción del calor, el cual puede variar con x, y ( , ) es la fuente de calor (o sumidero, si < 0). u(x,0)=u0(x) distribución inicial de temperatura
u0(x)
u(a,t)=α(t)
u(b,t)=β(t) x
κ(x) L x=a
φ(x,t)
x=b
fuente (o sumidero) de calor
F. I. Saldaña y W. J. Dután, estudiantes de la Maestria en Planificación y Gestión Energética en la Universidad de Cuenca, actualmente prestan sus servicios en la Empresa Eléctrica Azogues.. W. N. Bernal, estudiante de la Maestria en Planificación y Gestión Energética en la Universidad de Cuenca, desarrolla proyectos de Ingeniería Industrial.
Fig. 1. Transferencia de calor en una barra
Se asume que la teoría básica de esta ecuación es familiar al lector.
2
B. Condiciones de Frontera Si el material es homogéneo, entonces κ ( ) ≡ κ es independiente de x y la ecuación del calor se reduce a ( , )= κ ( , )+ ( , ) (2) Al inicio del análisis se requieren condiciones iniciales ( , 0) = ( ) (3) y condiciones de frontera, por ejemplo la temperatura podría ser especificada en cada extremo ( , ) = ( ), ( , )= ( ) (4) Tales condiciones de frontera, donde el valor de la solución en sí es especificado, son llamadas condiciones de frontera de Dirichlet. Alternativamente, en uno o ambos extremos pueden especificarse los valores de los gradientes normales sobre el contorno, en cuyo caso = en este punto. Esta condición de frontera, la cual es una condición sobre la derivada de u más que sobre u en sí, es llamada una condición de frontera de Neumann. C. El Problema de Estado Estable En general, esperamos que la distribución de temperatura cambie con el tiempo. Sin embargo, si ( , ), ( ) y ( ) son independientes del tiempo, entonces cabría esperar que la solución eventualmente alcance una situación de estado estable ( ), la cual entonces permanece esencialmente sin cambio en tiempos posteriores. Típicamente, habrá un tiempo transitorio inicial, a medida que ( ) se aproxima a ( ) (a menos que ( ) ≡ ( )), pero, si estamos únicamente interesados en determinar la solución de estado estable en sí, entonces podemos establecer que = 0 y obtener una ecuación diferencial ordinaria en x y la solución u(x) − ( )= ( ) (5) donde introdujimos ( )=
( )
κ
las derivadas en la ecuación usando la expansión en series de Taylor truncada. La aproximación discreta resulta en un conjunto de ecuaciones algebraicas que son evaluadas (o solucionadas) para los valores de las incógnitas discretas. El dominio físico es discretizado por medio de la mallado del mismo, esto es, dividir el dominio en subintervalos. La malla es un conjunto de puntos donde la solución discreta será computada. Estos puntos son llamados nodos. El parámetro clave de la malla es ∆ , la distancia entre puntos adyacentes. La idea del método de diferencias finitas es reemplazar las derivadas continuas con las llamadas fórmulas de diferencias que envuelven solamente valores discretos asociados con posiciones en la malla. En el presente caso existen derivadas únicamente con respecto al espacio. En el límite cuando el espaciamiento de la malla tiende a cero, la solución numérica obtenida con cualquier esquema útil aproximará la solución a la ecuación diferencial original. Sin embargo, la tasa con la cual la solución numérica se aproxima a la solución real varía con el esquema. A continuación se ilustra el método de solución mediante diferencias finitas centrales como una secuencia de pasos. 1. Discretización del dominio (generación de la malla) En una dimensión todos los problemas toman lugar sobre un intervalo finito del eje x. Realizamos una partición uniforme del dominio [a,b] en − 1 subintervalos de longitud uniforme ∆ = y ! puntos de malla " = 0,∙∙∙ . De aquí " = + $∆ , etc.
(6)
Esta es una ecuación diferencial ordinaria de segundo orden que puede ser solucionada de manera analítica, mediante funciones de Green o integrales de Fourier. Además, se requieren dos condiciones de frontera para determinar una solución única. En nuestro caso se aplican las condiciones de frontera de Dirichlet y Neumann ( )= ( )= , (7) Se tiene, por lo tanto, un problema de frontera de dos puntos, puesto que una condición es establecida en cada uno de los extremos del intervalo donde la solución es deseada. La función φ(x) es especificada y deseamos determinar u(x) en el intervalo < < . Dicho problema es conocido como de Poisson en una dimensión, con condiciones de frontera mixtas (de Dirichlet y Neumann), y se puede solucionar explícitamente en ciertos casos, pero en la aplicación a sistemas físicos reales se emplean métodos numéricos en los que la ecuación en derivadas parciales es reemplazada con una aproximación discreta. D. Método de Diferencias Finitas El método de diferencias finitas es una de varias técnicas para obtener la solución de la ecuación de calor, que aproxima
Fig. 2. Discretización del dominio en una dimensión
2. Aproximación por diferencias finitas Para la primera derivada ( " ) = lim∆
( " ) = lim∆ ( " ) = lim∆ ( ") ≈
*( + ,∆ ) *( + ) ∆
→
→
→
*( + ) *( + ∆ ) ∆ *( +-. ) *( +/. ) 0∆
*+-. *+/. 0∆
hacia adelante hacia atrás
(8) (9)
central
(10)
diferencia central
(11)
( ") ≈
2+-. 2+/. 0∆
(12)
Note la diferencia entre ui (exacta) y Ui (aproximación), esto es, " = ( " ) ≈ 3" , donde Ui son los escalares incógnitas que estamos buscando. Para la segunda derivada ( " ) = lim∆
→
* +-. * +/. ∆
(13)
3
( " ) = lim∆
9 /9 9 /9 456∆7→8 +-. + 456∆7→8 + +/. ∆7 ∆7
→
∆
( ") ≈
*+/. 0*+ ,*+-.
(14) (15)
(∆ ):
sustituyendo esta expresión en (5), obtenemos −
2+/. 02+ ,2+-. (∆ ):
=
"
= ( " ), ∀$ = 1, … ,
(16)
La condición de frontera de Dirichlet (7) implica que 3 =
(17)
mientras que la condición de frontera de Neumann puede ser aproximada por diferencias hacia atrás − −
=
( )= ( ) κ =
= κ = lim∆
−
→
(18)
*( > ) *( >/. ) ∆
(19)
= 2> 2>/. κ ≈ ∆ ⇒
3 −3
!
=−
(20)
= κ ∆x
(21)
Esta diferencia hacia atrás produce una inconsistencia, puesto que introduce un error de A(∆ ) mientras que las diferencias centrales de la ecuación de gobierno introducen un error de A(∆ )0 ). Esto puede solucionarse tratando la discretización de la condición de borde de Neumann de una manera diferente. Primero, introducimos un punto de la malla ficticio +1 = + ∆ , con una temperatura asociada +1 . Esta temperatura no tiene significado físico alguno, dado que el punto +1 se encuentra fuera del dominio del problema. Obtenemos un sistema de N-1 ecuaciones para N incógnitas. La ecuación faltante la provee la condición de Neumann, pero aproximando mediante diferencias centrales: 3 +1 −3 −1 2∆
≈−
(22)
κ
Combinando (16), (17), y (22), obtenemos el sistema de ecuaciones algebraicas lineales: DE = DE C3 2 −1 ⋯ H−1 2 −1 G G 0 −1 2 ⋮ ⋱ ⋱ C=G G 0 ⋯ 0 G ⋮ ⋯ ⋮ G 0 ⋯ ⋯ F 0 ⋯ ⋯ 3! H 3 N G ! M G 30 M DE = G ⋮ M 3 G 3" M G ⋮ M G3 ! M F 3 L
⋯ 0 −1 ⋱ −1 ⋱ 0 ⋯
(23) ⋯ ⋯ 0 ⋱ 2 ⋮ 0 0
⋯ ⋯ ⋯ ⋱ −1 ⋱ −1 1
⋯ ⋯ ⋯ ⋮ 0 ⋱ 2 0
0 0N M 0M ⋮ M ⋯M 0M −1M −1L
(∆ )0 ! + 3 H 0 G (∆ ) 0 G (∆ )0 O G ⋮ DE = G (∆ )0 " G ⋮ G 0 G (∆ ) P G 2∆ F κ
N M M M M M M M M L
DE es un vector de donde A es la matriz de coeficientes, 3 incógnitas (valores de temperatura), y DE es un vector de valores conocidos. DE = DE 3. Solucionar el sistema lineal C3 Las derivadas en la ecuación diferencial han sido reemplazadas con aproximaciones de diferencias finitas en cada punto interior de la malla, reduciendo la ecuación diferencial en un set de ecuaciones algebraicas acopladas, las cuales pueden entonces ser solucionadas mediante un software. E. Implementación en MatLab El uso general de paquetes de software permite solucionar problemas de una manera más fácil, rápida y precisa. Los modelos matemáticos utilizados en la representación de fenómenos físicos, requieren de análisis especiales tanto en los diferentes métodos aplicados para solucionarlos como en los resultados que se obtienen. MATLAB proporciona gran funcionalidad en este campo e incluso incorpora un toolbox para EDP. Se implementó un script en MATLAB que resuelve la ecuación del calor en una dimensión en estado estable, con una fuente de calor, y con condiciones de frontera mixtas, mediante diferencias finitas. Las entradas del código de MATLAB requieren: • Los extremos de la barra, a y b, • El coeficiente de difusividad térmica κ [m2/s] • El número de puntos, incluyendo las condiciones de frontera, en los que se aproximará la solución. • La función fuente de calor ( ). • La condición de frontera de Dirichlet en el extremo izquierdo de la barra. • La condición de frontera de Neumann en el extremo derecho de la barra. Seguidamente se procede a la solución mediante diferencias centrales. Los resultados se presentan de manera tabulada y gráfica. El código del script implementado calorv.m permite, además la solución de la ecuación (5) con condiciones de Dirichlet en ambos extremos de la barra. Para este caso se debe ingresar calorv(1). En el caso de condiciones mixtas se utiliza la rutina calorv(2). F. Aplicación El código implementado en el script calorv.m se prueba mediante la aplicación al siguiente problema básico, y sus variantes: Se tiene una barra de un material cuya conductividad es Q = 0.5T0 /V, la longitud es L=1 m, y la generación de calor interno es ( ) = 10W ⁄TO . El extremo izquierdo está expuesto a una temperatura de 25 ⁰C, y el otro extremo a 5 ⁰C. Para solucionar el problema se considera que el extremo izquierdo se encuentra en el origen del eje x, con lo cual = 0 (sin embargo, a puede tomar otro valor). Como Y = 1, = 1.
4
Flujo de Calor en una Barra
Al especificar 5 puntos de aproximación se obtienen los resultados de la tabla I:
25
TABLA I SOLUCIÓN DE LA ECUACIÓN EN CINCO PUNTOS
Temperatura
x 0 0.2 0.4 0.6 0.8 1
20
Temperatura 25 22.6 19.4 15.4 10.6 5
15
10
5
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Longitud de la Barra Fig. 4. Temperatura para aproximación en diez puntos.
Flujo de Calor en una Barra 25
Si al problema anterior se consideran condiciones mixtas, es decir de Dirichlet en el extremo izquierdo y de Neumann en el derecho (está aislado ( = 0)) se obtiene:
20
Temperatura
0
Flujo de Calor en una Barra
15
34 33 10
5
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Longitud de la Barra
Fig. 3. Temperatura para aproximación en cinco puntos.
Con 10 puntos de evaluación se obtiene:
Temperatura
32 31 30 29 28 27 26
TABLA II SOLUCIÓN DE LA ECUACIÓN EN DIEZ PUNTOS
x 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Temperatura 25 23.9 22.6 21.1 19.4 17.5 15.4 13.1 10.6 7.9 5
25
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Longitud de la Barra Fig. 5. Temperatura para aproximación en diez puntos con condiciones: 25 de Dirichlet y 0 de Neumann. TABLA III SOLUCIÓN DE LA ECUACIÓN EN DIEZ PUNTOS CON FUENTE DE CALOR SENOIDAL
x 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Temperatura 25 26.7 28.2 29.5 30.6 31.5 32.2 32.7 33 33.1 33
Para el caso de una fuente variable ( ) = 2 0 VZ[( ), el extermo izquierdo a una temperatura de 25⁰C y el extremo derecho aislado ( = 0):
5
Flujo de Calor en una Barra
G. Análisis de Resultados El dominio se es 0
25.35
25.3
Temperatura
25.25
25.2
25.15
25.1
\]
25.05
25
flujo de calor en la ecuación del calor es proporcional a \^ , así por ejemplo en x=1 el extremo aislado da una condición de 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
\](!)
1
Longitud de la Barra Fig. 6. Temperatura para aproximación en cinco puntos con condiciones: 25 de Dirichlet y 0 de Neumann. Función forzante: ( ) = 2 0 VZ[( )
Flujo de Calor en una Barra 25.7
2 −1 H−1 2 G C = G 0 −1 0 G0 F0 0
25.6
Temperatura
25.5
25.4
25.2
25.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Longitud de la Barra Fig. 7. Temperatura para aproximación en cien puntos con condiciones: 25 de Dirichlet y 0 de Neumann. Función forzante: ( ) = 2 0 VZ[( )
Flujo de Calor en una Barra 60
55
Temperatura
50
45
40
35
30
25
0 0 −1 0 2 −1 −1 2 0 1
0 0 NM 0M −1M −1L
25.0013 H 0.0100 N G M = G 0.0325 M G 0.0735 M F 0 L
25.3
25
frontera de = 0. \^ La implementación del método en MATLAB depende de los nodos discretizados. La ejecución del script discretiza la longitud de de la barra, y encuentra los valores Ui que constituyen un conjunto de soluciones aproximadas. Para el primer caso de fuente variable y cinco nodos:
0
1
2
3
4
5
6
7
8
9
10
Longitud de la Barra Fig. 8. Temperatura para aproximación en diez puntos con condiciones: 25 de Dirichlet y 0 de Neumann. Longitud: 10. Función forzante ( ) =
2 VZ[( ).
25.0805 H25.1597N G M 3 = G25.2290M G25.2657M F25.2290L H. Conclusiones y desarrollo posterior El código implementado resuelve la ecuación (5) mediante diferencias finitas, la cual es una de varias técnicas para obtener su solución. En todas las soluciones numéricas, la ecuación en derivadas parciales es reemplazada con una aproximación discreta. En este contexto, la palabra “discreta” significa que la solución numérica es conocida solamente en un número finito de puntos en el dominio físico. El número de estos puntos puede ser seleccionado por el usuario del método numérico. De los resultados obtenidos para las diferentes corridas, en general se puede indicar que incrementar el número de puntos no solamente incrementa la resolución (detalle), sino que además incrementa la precisión de la solución numérica. Esto debido a que una fuente de error es el error de truncamiento, que puede ser reducido al incrementar el número de nodos y reducir el tamaño de los subintervalos. Con el fin de que los resultados no sean altamente
6
oscilatorios se deben elegir pasos razonablemente pequeños, de tal manera que el esquema sea estable. Es importante indicar que el problema analizado es estrictamente de condiciones de frontera y no requiere de condiciones iniciales. Además de las diferencias finitas existen otros métodos de solución, de los cuales el más importante es el de elementos finitos. Se plantea como desarrollo posterior la resolución de la ecuación de transferencia de calor generalizada, es decir considerando la evolución temporal (transitoria), y además usando el método de Crank-Nicholson.
IV. BIBLIOGRAFÍA [1] [2]
[3] [4] [5]
[6]
B. Hunt, R. Lipsman, J. Rosenberg, A Guide to MATLAB for Beginners and Experienced Users Systems, Cambridge PRESS, 2001. H. J. Lee, W. E. Schiesser, Ordinary and partial differential equation routines in C, C++, Fortran, Java, Maple, and MATLAB, USA, Chapman & Hall/CRC, 2004. S. C. Chapra, R. P. Canale, Numerical Methods for Engineers. (5ª ed.) Singapur, McGrawHill, 2006. G. D. Smith, Numerical Solution of Partial Differential Equations, Oxford, 3rd ed., 1985. W. J. Minkowycz, E. M. Sparrow and J. Y. Murthy, Handbook of Numerical Heat Transfer, Second Edition, 2006 John Wiley & Sons, Inc. http://www.mathworks.com/moler
V. ANEXO Se presenta el código del script calorv.m. %Programa que permite resolver la transferencia de calor en estado %estacionario en una dimensión por el método de diferencias finitas. % %Para utilizar el programa debe definir el tamaño de una barra, esto por %tratarse de una dimensión, tambien se debe establecer el material del cual %esta constituido, para lo cual se define el coeficinte de difusividad %térmica Kc. % %Como el método que se utiliza es el de diferencias finitas, se debe %considerar el número de puntos a estudiar o malla a construir. % %El problema permite estudiar el sistema cuando contiene una fuente %de calor intrna, que depende de la posición, por lo tanto este parámetro %tambien se debe ingresar. % %Al ser el problema una ecucación diferencial, se definen dos tipos de %condiciones de borde, la de Dirichlet para x(a)=alfa y la de Neumann %Ux=beta,|x=L, donde alfa es el valor inicial en el inicio de la barra y %beta es la condición de borde de la derivada en el otro extremo de la %barra (0 en el caso homogéneo). %
%Al finalizar los cálculos se presenta la grafica de la variación de %temperatura a lo largo de la barra aplicando las condiciones ingresadas. function calorv(tipo) %Menù para escoger los casos solicitados switch tipo case (1) c1; %Resuelve las ecuaciones diferenciales y devuelve un gràfica de X vs. Y case (2); c2; %Resuelve las ecuaciones diferenciales y devuelve un gràfica de X,Y y Z vs t otherwise error('Función no definida'); end
function c1 clear all clc global y s lamda disp('Programa para resolver la transferencia de calor estacionario en una'); disp('dimensión con condiciones de Dirichlet'); disp('por el método de diferencias finitas'); disp(' '); disp('Datos:'); a=input('Ingrese el inicio de la barra:'); L=input('Ingrese el final de la barra:'); %Es necesario conocer la %longitud de la barra para el càlculo en una dimensión k=input('Ingrese el coeficiente de difusividad térmica Kc:');%El coeficiente %de conductividad térmica es característico del material; n=input('Ingrese el numero de puntos en los que desea analizar el sistema:'); %Se deben incluir los bordes q=input('Ingrese la funciòn de calor q(x)=','s');%Es la funciòn que representa %la temperatura inicial D=input('Ingrese la condicion inicial para x(a):');%Se aplica la condición %de frontera de Dirichlet al inicio de la barra N=input('Ingrese la condicion inicial para x(L):');%Se aplica la %condición de frontera de Dirichlet al final de la barra y=linspace(a,L,n) %Crea un vector con incremento s hasta L que representa %los puntos de análisis del sistema. s=(y(1,2)-y(1,1)); %diferencial de desplazamiento que depende del nùmero %de puntos escogidos para el análisis lamda=s^2/k; %Ahora se va a formar el arrelgo AU=q %Creación del vector q, el cual contine las condiciones de Dirichlet for i=2:n-2 if i==2 x=y(1,i); b(i,1)=lamda.*eval(q)+D;%Condición de Dirichlet al inicio de la barra. else x=y(1,i); b(i,1)=lamda.*eval(q); end end x=y(1,n-1); b(n-1,1)=lamda.*eval(q)+N; %Condición de Dirichlet al final de la barra. b(1,:)=[]; b %Creaciòn de la matriz A, que contiene los coeficientes de la variables
7 %de temperatura en cada punto de análisis. for i=1:n-2 for j=1:n-2 if i==j A(i,j)=2; end end end for i=1:n-3 for j=1:n-3 if i==j j=j+1; A(i,j)=-1; end end end for i=2:n-2 for j=2:n-2 if i==j j=j-1; A(i,j)=-1; end end end
A % Vector tridiagonal que contiene los coeficientes de las ecuaciones de %temperatura T=A\b; U=[D;T;N] plot(y,U) %Gráfica de la Temperatura a lo largo de la barra grid on xlabel('Longitud de la Barra','FontSize',14);ylabel('Temperatura','FontSize',1 4);title('Flujo de Calor en una Barra','FontSize',16); function c2 clear all clc global y s lamda disp('Programa para resolver la transferencia de calor estacionario en una'); disp('dimensión, con condiciones mixtas (Dirichlet y Neumann),'); disp('por el método de diferencias finitas'); disp(' '); disp('Datos:'); a=input('Ingrese el inicio de la barra:'); L=input('Ingrese el final de la barra:'); %Es necesario conocer la %longitud de la barra para el càlculo en una dimensión k=input('Ingrese el coeficiente de difusividad térmica Kc:');%El coeficiente %de difusividad térmica es característico del material; n=input('Ingrese el numero de puntos en los que desea analizar el sistema:'); %Se deben incluir los bordes q=input('Ingrese la funciòn de calor q(x)=','s');%Es la funciòn que representa %la temperatura inicial D=input('Ingrese la condicion inicial para x(0):');%Se aplica la condición %de frontera de Dirichlet N=input('Ingrese la condicion de contorno para dx/dt|x=L:');%Se aplica la %condición de contorno de Neumann y=linspace(a,L,n) %Crea un vector con incremento s hasta L que representa %los puntos de análisis del sistema. s=(y(1,2)-y(1,1)); %diferencial de desplazamiento que depende del nùmero %de puntos escogidos para el análisis lamda=s^2/k; %Ahora se va a formar el arrelgo AU=q
%Creación del vector q, el cual contine las condiciones de Dirichlet y Neumann for i=2:n-1 if i==2 x=y(1,2); c(i,1)=lamda.*eval(q)+D; %Condición de Dirichlet. else x=y(1,i); c(i,1)=lamda.*eval(q); end end g=(2*s/k)*N; %Condición de Neumann. c(1,:)=[]; b=[c;g] % Vector q %Creaciòn de la matriz A, que contiene los coeficientes de la variables %de temperatura en cada punto de análisis. for i=1:n-2 for j=1:n-2 if i==j A(i,j)=2; end end end for i=1:n-2 for j=1:n-2 if i==j j=j+1; A(i,j)=-1; end end end for i=2:n-2 for j=2:n-2 if i==j j=j-1; A(i,j)=-1; end end end %Se completa el vector A con los elementos dbajo mostrados en razòn de que %se utilizó diferencias centrales para representar la derivada de primer %grado y al hacer esto se introdujo una nueva variable para la cual se %debe aplicar la condición de Neumann A(n-1,n-1)=-1; A(n-1,n-3)=1; A % Vector tridiagonal que contiene los coeficientes de las ecuaciones de %temperatura T=A\b U=[D;T] % Vector de resultados de Temperatura en la barra. plot(y,U) %Gráfica de la Temperatura vs. el desplazamieno grid on xlabel('Longitud de la Barra','FontSize',14);ylabel('Temperatura','FontSize',1 4);title('Flujo de Calor en una Barra','FontSize',16);
