• Author / Uploaded
• Cristian Untitled

Citation preview

1

Resolutores de Ecuaciones Diferenciales Mathcad tiene una variedad de funciones para resolver numéricamente ecuaciones diferenciales parciales: Resolutor de Ecuaciones Diferenciales Ordinarias odesolve(x,b,step) rkfixed(y,x1,x2,npoints,D) Sistemas Alisados Bulstoer(y,x1,x2,npoints,D) Sistemas Stiff Stiffb(y,x1,x2,npoints,D,J) Stiffr(y,x1,x2,npoints,D,J) Sistemas lentamente variables Rkadapt(y,x1,x2,npoints,D) Encontrar el ultimo punto sobre el intervalo de integración bulstoer(y,x1,x2,acc,D,kmax,s) rkadapt(y,x1,x2,acc,D,kmax,s) stiffb(y,x1,x2,acc,D,J,kmax,s) stiffr(y,x1,x2,acc,D,J,kmax,s) Resolver problemas de valores límites dos-puntos bvalfit(v1,v2,x1,x2,xf,D,load1,load2,score) sbval(v,x1,x2,D,load,score) Resolver Ecuaciones Diferenciales a las derivadas parciales relax(a,b,c,d,e,f,u,rjac) multigrid(M,ncycle) Resolución de una sola ecuación diferencial ordinaria odesolve(x,b,[step]) Retorna una función de x la cual es una solución a la ecuación diferencial ordinaria (ODE), sujeta a a los constreñiminetos de Facultad de Ingeniería – Universidad de Mendoza

Dr. Ing. Jesús Rubén Azor Montoya

2 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS”

valor inicial o bordes provista en el solve block. La ODE debe ser lineal en sus derivadas más altas y el número de condiciones debe ser igual al orden de la ODE. . Argumentos: x es la variable de integración. x debe ser real. b es el punto terminal del intervalo de integración. b debe ser real. step (opcional) es el número de pasos usados internamente cuando se calcula la solución. Uso de la función odesolve: Pasos para usar la función odesolve para resolver una ecuación diferencial ordinaria: • Tipee la palabra Given para arrancar el solve block. •

Por debajo del Given, tipee la ecuación diferencial y sus constrñimientos usando boolean operators.

•

Tippee la función odesolve con la variable de integración, x, y el punto terminal b.

Notes: •

La ecuación diferencial puede ser escita usando los operadores derivada tales como d/dx y d2/dx2 o usando notación primada similar a y(x) e y'(x). (La combinación de teclas para el character prima es Ctrl+F7.)

•

Los constreñimientos estarán en la forma de y(a)=b o y'(a)=b. Mathcad no aceptará constrñimientos más complicados tales como y'(a)+y(a)=b.

•

El punto terminal b debe ser más grande que el valor inicial.

•

Por default, odesolve usa un runge-kutta de paso fijo de resolución. Para usar un método adaptivo, cliquee sobre odesolve con el botón derecho del mouse y eleija Adaptive desde el menú pop-up.

•

Para resolver sistemas de ecuaciones diferenciales o para resolver una que no es lineal en el término de derivada más alta, use rkfixed.

Ejemplo:

3

Los ejemplos de abajo demuestran cómo usar la función odesolve para resolver ecuaciones diferenciales ordinarias: Given

100 ⋅y''( x) + 10 ⋅y'( x) + 101 ⋅y( x)

y( 0)

0

y'( 0)

1  50 ⋅cos  ⋅x 4 

1

y := Odesolve( x , 150) 2 1.465 1 y( x) 0 − 0.698 −1

0

2

4

0

6

8

x

10 10

Given 4⋅

d2 2

f ( t) + f ( t)

t

f ( 0)

4

f ( 5)

13.5

dt

f := Odesolve( t , 5)

30 20 f ( t) 10

0

1

2

3

4

5

6

t

---------------------------------------------------------------------------------------------------Utilización de Matlab para resolución de Ecuaciones Diferenciales

DSOLVE Solución simbólica de ecuaciones diferenciales ordnarias DSOLVE('eqn1','eqn2', ...) acepta ecuaciones simbólicas representando ecuaciones diferenciales ordinarias y condiciones iniciales. Varias ecuaciones o Facultad de Ingeniería – Universidad de Mendoza

Dr. Ing. Jesús Rubén Azor Montoya

4 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS”

condiciones iniciales pueden ser agrupadas juntas, separadas por comas, en un único argumento de entrada. Por omisión, la variable independiente es ' t '. La variable independente puede ser cambiada de 't' a alguna otra variable simbólica incluyendo esa variable como el último argumento. La letra 'D' denota derivada con respecto a la variable independiente, en este caso usualmente d/dt. Una 'D' seguida por un dígito denota derivación repetida; por ejemplo, D2 es d^2/dt^2. Cualesquiera caracteres siguiendo estos operadores de derivación son tomados como variables dependientes; por ejemplo, D3y denota al tercera derivada de y(t). Note que los nombes de las variables simbólicas no deberán contener la letra 'D'. Las condiciones iniciales son especificadas por ecuaciones tales como 'y(a)=b' o 'Dy(a) = b' donde y es una de las variables dependientes y a y b son constantes. Si el número de condiciones iniciales es menor que el número de variables dependientes, las soluciones resultantes obtendrán constantes arbitrarias, C1, C2, etc. Son posibles tres diferentes tipos de salidas. • Para una ecuación y una salida, es retornada la solución resultante, con soluciones múltiples para una ecuación no lineal en un vector simbólico. • Para varias ecuaciones e igual número de salidas, los resultados son ordenados en orden lexicográfico y asignados a las salidas. • Para varias ecuaciones y una única salida, se retorna una estructura conteniendo las soluciones. Si no se encuentra ninguna solución closed-form (explícita) , se intenta una solución implícita. Cuando se retorna una solución implícita, se da una advertencia. Si no se puede calcular una solución explícita o implícita, entonces se da una advertencia y se retorna el sym vacío. En algunos casos involucrando ecuaciones no-lineales, la salida será una ecuación diferencial de más bajo orden equivalente o una integral. Ejemplos: dsolve('Dx = -a*x') retorna ans = exp(-a*t)*C1 x = dsolve('Dx = -a*x','x(0) = 1','s') retorna x = exp(-a*s) y = dsolve('(Dy)^2 + y^2 = 1','y(0) = 0') retorna y= [ sin(t)] [ -sin(t)] S = dsolve('Df = f + g','Dg = -f + g','f(0) = 1','g(0) = 2') retorna una estructura S con campos S.f = exp(t)*cos(t)+2*exp(t)*sin(t)

5

S.g = -exp(t)*sin(t)+2*exp(t)*cos(t) Y = dsolve('Dy = y^2*(1-y)') Advertencia: No puede ser encontrada solución explícita; se retorna la solución implícita. Y= t+1/y-log(y)+log(-1+y)+C1=0 dsolve('Df = f + sin(t)', 'f(pi/2) = 0') dsolve('D2y = -a^2*y', 'y(0) = 1, Dy(pi/a) = 0') S = dsolve('Dx = y', 'Dy = -x', 'x(0)=0', 'y(0)=1') S = dsolve('Du=v, Dv=w, Dw=-u','u(0)=0, v(0)=0, w(0)=1') w = dsolve('D3w = -w','w(0)=1, Dw(0)=0, D2w(0)=0') y = dsolve('D2y = sin(y)'); pretty(y) Algunos comandos, tales como ode45 (un resolutor de ecuaciones diferenciales en forma numérica), requiere que su primer argumento sea una funcción — para ser preciso o bien una función inline, como en ode45(f, [0 2], 1). o una function handle, esto es, el nombre de una función built-in o una función M-file precedida por el símbolo especial @, como en ode45(@func, [0 2], 1)). Ecuaciones Diferenciales Ordinarias El objetivo de este laboratorio es aprender técnicas para la resolución numérica de problemas de valores iniciales (P.V.I.) para ecuaciones diferenciales ordinarias (E.D.O.) y sistemas de E.D.O. Matlab tiene varios comandos para la resolución numérica de P.V.I. para E.D.O.: Resolutores de ecuaciones diferenciales. Resolutores de problemas con valores iniciales para ODEs. (Si no tiene seguridad acerca del stiffness, intente primero ODE45, luego ODE15S.) [En matemática, una ecuación stiff es una ecuación diferencial en los que determinados

métodos numéricos para resolver la ecuación son numéricamente inestables, a menos que el tamaño de paso sea tomado extremadamente pequeño. Ha resultado difícil formular una definición precisa del stiff, pero la idea principal es que la ecuación incluye algunos términos que pueden conducir a una variación rápida en la solución.] ode45 – Resuelve ED non-stiff, por el medium order method. ode23 - Resuelve ED non-stiff, por el low order method. ode113 - Resuelve ED non-stiff,, por el variable order method. ode23t - Resuelve ED moderadamente non-stiff, y DAEs Index 1, trapezoidal rule. ode15s - Resuelve ED stiff y DAEs Index 1, variable order method. ode23s - Resuelve ED stiff, low order method. ode23tb - Resuelve ED stiff, low order method. Como se ve en esta lista, hay métodos para resolver E.D.O. stiff y no stiff. Además hay métodos de orden bajo, medio, alto y variable. Todos ellos tienen una sintaxis semejante. Por ejemplo, para resolver el P.V.I.

Facultad de Ingeniería – Universidad de Mendoza

Dr. Ing. Jesús Rubén Azor Montoya

6 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS”

en el intervalo [to, tf ] mediante el comando ode45 en su opción más sencilla, debe ejecutarse: [t,y]=ode45(’f’,[to tf],yo); donde: • f es el nombre de la función f(t, y) (típicamente definida mediante un programa function en un archivo f.m); • to y tf son los extremos del intervalo donde se desea conocer la solución; • yo es el valor de la solución en to (es decir el valor de la condición inicial y(to) = yo); • t devuelve los valores de la variable independiente t donde el método calcula el valor de la solución; • y devuelve los valores de la solución en cada uno de los puntos t. Estos comandos no requieren como dato un paso de integración h pues todos ellos determinan de manera automática en cada paso k, el tamaño del paso de integración hk necesario para mantener los errores por debajo de una tolerancia determinada. Los valores de t que entrega corresponden a los puntos tk = tk−1 + hk, k = 1, 2, . . . , en los que el comando necesitó calcular el valor de y(tk).

Si se desea conocer la solución para ciertos valores de t, puede alternativamente ejecutarse: [t,y]=ode45(’f’,tspan,yo); donde tspan es el vector de valores donde se desea conocer la solución. Por ejemplo, tspan=0:0.1:1. En ese caso, la salida t coincide con tspan e y contiene los valores de la solución en esos puntos. La tolerancia predeterminada de estos métodos es 10E−3, para el error relativo, y 10E−6, para el error absoluto. Si se desea calcular la solución con otras tolerancias, deben prefijarse las opciones elegidas mediante el comando odeset. Además, en la ejecución del comando para resolver la E.D.O., debe agregarse el parámetro adicional de opciones. La sintaxis para realizar esto es, por ejemplo: options=odeset(’RelTol’,1e-6,’AbsTol’,1.e-8); [t,y]=ode45(’f’,[to tf],yo,options); Si se ejecuta options=odeset(’RelTol’,1e-6,’AbsTol’,1.e-8) sin el “;” puede verse que hay otras opciones que pueden prefijarse, además de las tolerancias de los errores. Por ejemplo, si se desea resolver el P.V.I.

en el intervalo [0, 1.5], mediante el comando ode45 y visualizar la solución obtenida, debe crearse un fichero f.m como sigue:

7

function z=f(t,y) z=y; y ejecutarse: [t,y]=ode45(’f’,[0 1.5],1); plot(t,y,t,exp(t),'o') Así se obtiene la siguiente gráfica:

El siguiente ejemplo resuelve la misma ecuación en los puntos t=0:0.1:1.5, con error absoluto menor a 10−E6 y calcula los errores cometidos restando los valores calculados a los de la solución verdadera, que en este caso es y(t) = exp(t): options=odeset('AbsTol',1.e-6); tspan=0:.1:1.5; [t,y]=ode45('f',tspan,1,options); error=exp(t)-y error = 1.0e-006 * 0 -0.0003 -0.0248 -0.0448 -0.0076 -0.0415 -0.0694 -0.0200 -0.0669 -0.1056 -0.0402 -0.1048 -0.1586 -0.0721 -0.1612

Facultad de Ingeniería – Universidad de Mendoza

Dr. Ing. Jesús Rubén Azor Montoya

8 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS” -0.0989

La salida que se presenta indica que los errores son efectivamente menores en valor absoluto a 10E−6. La resolución de P.V.I. para sistemas de E.D.O. se realiza mediante los mismos comandos. En tal caso, f(t,y) debe ser una función a valores vectoriales (es decir un vector columna de funciones) e y un vector columna de variables de la misma dimensión. Además, la condición inicial yo también debe ser un vector columna de la misma dimensión. Por ejemplo, consideremos el P.V.I.

cuya solución exacta es:

Por lo tanto los puntos (x(t), y(t)) solución de este sistema de E.D.O, describen la circunferencia unitaria. Este sistema escrito vectorialmente resulta:

Para resolverlo debe crearse un fichero F.m como sigue:

function Z=F(t,Y) Z=[Y(2);-Y(1)]; Los siguientes comandos resuelven este P.V.I. en el intervalo [0, 2] y grafican la curva (x(t), y(t)), para 0 < t< 2, que se obtiene: [t,Y]=ode45('f',[0 2*pi],[1;0]); plot(Y(:,1),Y(:,2)); Así se obtiene la siguiente gráfica:

9

Supongamos que deseamos resolver y plotear la solución de la siguiente ecuación diferencial de segundo orden

eqn2='D2y+8*Dy+2*y=cos(x)'; inits2 = 'y(0)=0, Dy(0)=1'; y=dsolve(eqn2,inits2,'x') y = 1/65*cos(x)+8/65*sin(x)+(-1/130+53/1820*14ˆ(1/2))*exp((-4+14ˆ(1/2))*x) -1/1820*(53+14ˆ(1/2))*14ˆ(1/2)*exp(-(4+14ˆ(1/2))*x)

x=0:0.1:1; z=eval(vectorize(y)); plot(x,z)

Given

Facultad de Ingeniería – Universidad de Mendoza

Dr. Ing. Jesús Rubén Azor Montoya

10 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS”

1  50 ⋅cos  ⋅x 4 

100 ⋅y''( x) + 10 ⋅y'( x) + 101 ⋅y( x)

y( 0)

0

y'( 0)

1

y := Odesolve( x , 150) 2 1.465 1 y( x) 0 − 0.698 −1

0

2

0

4

6

8

x

10 10

Resolver con Matlab: eqn2='100*D2y+10*Dy+101*y=50*cos(0.25*x)'; inits2 = 'y(0)=0, Dy(0)=1'; y=dsolve(eqn2,inits2,'x') x=0:0.1:10; z=eval(vectorize(y)); plot(x,z)

Resolución de Sistemas de ED con Matlab Supóngase que deseamos resolver y graficar las soluciones del sistema de tres ED ordinarias

11

Primero, para encontrar una solución general, procedemos como en el caso de una sóla ED, excepto que cada ecuación es abrazada por su par de comillas (simples): [x,y,z]=dsolve(’Dx=x+2*y-z’,’Dy=x+z’,’Dz=4*x-4*y+5*z’) x = -C1*exp(3*t)-C2*exp(t)-2*C3*exp(2*t) y = C1*exp(3*t)+C2*exp(t)+C3*exp(2*t) z = 4*C1*exp(3*t)+2*C2*exp(t)+4*C3*exp(2*t) Tenga en cuenta que ya no se ha especificado ninguna variable independiente, MATLAB utiliza por defecto t.. Para resolver un problema de valores iniciales, simplemente se define un conjunto de valores iniciales y se añaden al final del comando dsolve (). Supongamos que tenemos x (0) = 1, y (0) = 2, y z (0) = 3. Luego: inits=’x(0)=1,y(0)=2,z(0)=3’; [x,y,z]=dsolve(’Dx=x+2*y-z’,’Dy=x+z’,’Dz=4*x-4*y+5*z’,inits) x = -5/2*exp(3*t)-5/2*exp(t)+6*exp(2*t) y = 5/2*exp(3*t)+5/2*exp(t)-3*exp(2*t) z = 10*exp(3*t)+5*exp(t)-12*exp(2*t) Finalmente, graficando esta solución t=linspace(0,.5,25); xx=eval(vectorize(x)); yy=eval(vectorize(y)); zz=eval(vectorize(z)); plot(t, xx, t, yy, t, zz)

Facultad de Ingeniería – Universidad de Mendoza

Dr. Ing. Jesús Rubén Azor Montoya

12 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS”

Búsqueda de soluciones numéricas MATLAB tiene una serie de herramientas para resolver numéricamente las ecuaciones diferenciales ordinarias. Nos centraremos en los dos principales, el incorporado en las funciones ode23 y ode45, que implementan versiones de Runge-Kutta 2do/3er-orden y 4to/5to-orden, respectivamente. Ejemplo. Aproximar numéricamente la solución de el ED de primer orden

sobre el intervalo x  [0, .5]. Para cualquier ED en la forma y′ = f(x, y), comenzamos definiendo la función f(x, y). Para ecuaciones únicas, podemos definir f(x, y) como una función inline. Aquí, f=inline('x*y^2+y') f = Inline function: f(x,y) = x*yˆ2+y El uso básico para el resolutor ode45 de Matlab es ode45(function,domain,initial condition). Esto es, usamos [x,y]=ode45(f,[0 .5],1) y MATLAB retorna dos vectores columna, el primero con valores de x y el segundo con valores de y. Ya que x e y son vectores con los correspondientes componentes, podemos graficar los valores con plot(x,y)

13

Elección de la partición. En la aproximación de esta solución, el algoritmo ode45 ha seleccionado una partición determinada del intervalo [0, 0.5], y MATLAB ha devuelto un valor de y en cada punto de esta partición. A menudo es el caso en la práctica en que nos gustaría especificar la partición de valores en los que MATLAB devuelve una aproximación. Por ejemplo, sólo puede ser que desee para aproximar y(0.1), y(0,2), ..., y(0,5). Podemos especificar esto al introducir el vector de valores [0, 0.1, 0.2, 0.3, 0.4, 0.5], como el dominio en el ode45. Es decir, que utilizamos xvalues=0:.1:.5; [x,y]=ode45(f,xvalues,1) x = 0 0.1000 0.2000 0.3000 0.4000 0.5000 y = 1.0000 1.1111 1.2500 1.4286 1.6667 2.0000 Opciones. Hay varias opciones disponibles para el resolutor ode45, dando al usuario un control limitado sobre el algoritmo. Dos opciones importantes son la tolerancia relativa y absoluta, respectivamente RelTol y AbsTol. En cada paso del algoritmo ode45, un error se aproxima a ese paso. Si yk es la aproximación de y(xk) en el paso k, y ek es el error aproximado en este paso, a continuación, MATLAB elige su partición para asegurar

Facultad de Ingeniería – Universidad de Mendoza

Dr. Ing. Jesús Rubén Azor Montoya

14 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS”

ek ≤ max (yk * RelTol , AbsTol), donde los valores por defecto son RelTol = 0,001 y AbsTol = 0.000001. Como un ejemplo de cuándo puede ser que desee cambiar estos valores, observamos que si yk llega a ser grande, entonces al error ek se le permitirá crecer bastante. En este caso, aumentar el valor de RelTol. Para la ecuación y' = xy2 + y, con y (0) = 1, los valores de y llegan a ser muy grandes cuando x se acerca a 1. De hecho, con las tolerancias de error por defecto, nos encontramos con que el comando [x, y] = ode45 (f, [0,1], 1); conduce a un mensaje de error, causado por el hecho de que los valores de y son cada vez más grandes a medida que x se acerca a 1. (Note en la parte superior del vector de la columna para y que se multiplica por 1014.) Con el fin de solucionar este problema, seleccione un valor menor para RelTol. options=odeset ('RelTol', 1e-10); [x, y] = ode45 (f, [0,1], 1,options); max (y) ans = 2.425060345544448e 07 Además de emplear el comando option, se ha calculado el valor máximo de y(x) para mostrar que o mostrar que sí es bastante grande, aunque no tan grandes como se sugiere en los últimos cálculos. Ecuaciones de primer orden con M-files Alternativamente, se puede resolver la misma ODE definiendo primero a f(x, y) como un M-file firstode.m. function yprime = firstode(x,y); % FIRSTODE: Computes yprime = x*yˆ2+y yprime = x*yˆ2 + y; En este caso, sólo requerimos un cambio en el commando ode45: debemos usar un puntero @ para indicar el M-file. Esto es, usamos los siguientes comandos. xspan = [0,.5]; y0 = 1; [x,y]=ode23(@firstode,xspan,y0); plot(x,y)

15

Sistemas de EDs Resolver un sistema de EDs en MATLAB es muy similar a resolver una única ecuación, aunque un sistema de EDs no puede ser definido como una función inline debemos definirlo como un M-file. Ejemplo. Resolver el sistema de ecuaciones de Lorenz

[Las ecuaciones de Lorenz tienen algunas propiedades de las ecuaciones en derivados atmosféricos. Las soluciones de las ecuaciones de Lorenz han servido como ejemplo de comportamiento caótico.] Donde a los efectos de este ejemplo, vamos a tomar σ = 10, β = 8 / 3, y ρ = 28, así como x (0) = - 8, y (0) = 8, y, z (0) = 27. El M-file que contiene las ecuaciones de Lorenz aparece a continuación. function xprime = lorenz1(t,x); %LORENZ: Computes the derivatives involved in solving the %Lorenz equations. sig=10; beta=8/3; rho=28; xprime=[-sig*x(1) + sig*x(2); rho*x(1) - x(2) - x(1)*x(3); -beta*x(3) + x(1)*x(2)]; Observe que x se almacena como x(1), y como x(2), y z como x(3). Además, xprime es un vector columna, como se desprende de la coma después de la primera aparición de x(2). Si en la ventana de comandos, escribimos

Facultad de Ingeniería – Universidad de Mendoza

Dr. Ing. Jesús Rubén Azor Montoya

16 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS”

x0=[-8 8 27]; tspan=[0,20]; [t,x]=ode45(@lorenz1,tspan,x0); Aunque no se dan aquí, la salida de este último comando consiste en una columna de las tiempos seguido por una matriz con tres columnas, la primera de las cuales se corresponde con los valores de x en los tiempos asociados, y lo mismo para la segunda y tercera columna para y y z. La matriz se ha denotado x en la declaración llamante ode45, y en general cualquier coordenada de la matriz se puede especificar como x(m,n), donde m denota la fila y n denota la columna. En lo que se estará más interesado es en referencia a las columnas de x, las cuales se corresponden con valores de los componentes del sistema. En este sentido, podemos denotar todas las filas o todas las columnas por un colon (:). Por ejemplo, x(:,1) se refiere a todas las filas de la primera columna de la matriz x, es decir, se refiere a todos los valores de nuestro componente original x. Con esta información, podemos fácilmente trazar el atractor extraño de Lorenz, que es un gráfico de z en función de x: plot(x(:,1),x(:,3))

Por supuesto, también podemos trazar cada componente de la solución en función de t, y una forma útil de hacer esto es apilar los resultados. Podemos crear la siguiente figura con: subplot(3,1,1) plot(t,x(:,1)) subplot(3,1,2) plot(t,x(:,2)) subplot(3,1,3) plot(t,x(:,3))

17

Pasando parámetros Al analizar el sistema de ecuaciones diferenciales, a menudo se quiere experimentar con diferentes valores de los parámetros. Por ejemplo, en el estudio de las ecuaciones de Lorenz se podría considerar el comportamiento en función de los valores de σ , β y ρ. Por supuesto, una forma de cambiar esto es manualmente volviendo a abrir el M-file lorenz1.m cada vez que se quiere probar con nuevos valores, pero no sólo es una forma lenta de hacerlo, sino que es difícil de manejar para automatizar. Lo que podemos hacer en cambio es pasar valores de los parámetros directamente a nuestro M-file a través de la instrucción de llamada ode45. Para ver cómo funciona esto, lo primero es alterar lorenz1.m en lorenz2.m, el último de los cuales acepta un vector de parámetros que denotamos con p. funtion XPRIME = lorenz2 (t, x, p); % LORENZ: Calcula los derivados necesarios para resolver el % ecuaciones de Lorenz. sig = p (1), beta = p (2); rho = p (3); XPRIME = * [-señal x (1) + * señal x (2); * rho x (1) - x (2) - x (1) * x (3)-beta * x (3) + x (1) * x (2)]; Ahora puede enviar valores de los parámetros con ode45. > p> = [10 08/03 28]; >> [t, x] = ode45 (@ lorenz1, tspan, x0, [], p); Podemos enviar ahora los valores de parámetros con ode45. p=[10 8/3 28]; [t,x]=ode45(@lorenz2,tspan,x0,[],p); Ecuaciones de Segundo Orden El primer paso en resolver una ED ordinaria de segundo orden (o más) en MATLAB es escribir la ecuación como un sistema de primer orden. Como ejemplo se retomará uno anterior. Tomando y1(x) = y(x) e y2(x) = y′(x), tenemos el sistema Facultad de Ingeniería – Universidad de Mendoza

Dr. Ing. Jesús Rubén Azor Montoya

18 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS”

Método visto arriba: Codificado en Matlab (usando dsolve, solución simbólica): eqn2 = 'D2y + 8*Dy + 2*y = cos(x)'; inits2 = 'y(0)=0, Dy(0)=1'; y=dsolve(eqn2,inits2,'x') y = exp((-4+14^(1/2))*x)*(53/1820*14^(1/2)-1/130)+exp(-(4+14^(1/2))*x)*(53/1820*14^(1/2)-1/130)+1/65*cos(x)+8/65*sin(x)

x=0:0.01:1; z = eval(vectorize(y)); plot(x,z)

Método sistema de EDs (usando ode45, solución numérica) Primero, se construye el M-file basado en:

function xprime = ee(x,y); % % xprime=[y(2);-8*y(2)-2*y(1)+cos(x)];

Se ejecuta:

19

x0=[0 1]; % condiciones iniciales tspan=[0,1]; % intervalo [x,y]=ode45(@ee,tspan,x0); plot(x,y(:,1))

Transformadas de Laplace Una de las más útiles transformadas en matemática es la transformada de Laplace. MATLAB tiene rutinas built-in para calcular la Transformada de Laplace como su inversa. Por ejemplo, para computar la transformada de f(t)=t2, simplemente tipee syms t; laplace(t^2) ans = 2/s^3 Para invertir, digamos, F(s) = 1/(1+s), tipee syms s; ilaplace(1/(1+s)) ans = exp(-t) Problemas de contorno Por diversas razones de mérito discutible mayoría de los cursos de introducción a ecuaciones diferenciales ordinarias se centran principalmente en problemas de valores iniciales (IVP). Otra clase de EDs que surgen a menudo en las aplicaciones son los problemas de contorno (BVPs). Consideremos, por ejemplo, la ecuación diferencial y''- 3y '+ 2y = 0 y (0) = 0 y (1) = 10, Facultad de Ingeniería – Universidad de Mendoza

Dr. Ing. Jesús Rubén Azor Montoya

20 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS”

donde nuestras condiciones y(0) = 0 e y(1) = 10 se especifican en el límite del intervalo de interés x ε [0, 1]. (Aunque nuestra solución normalmente se extiende más allá de este intervalo, la mayoría de los escenarios comunes en los problemas de valores en la frontera es el caso en el que sólo estamos interesados en los valores de la variable independiente entre los puntos extremos especificados.) El primer paso en la solución de este tipo de ecuación es escribirlo como un sistema de primer orden con y1 = y e y2 = y', por lo cual tenemos

Grabamos este sistema en el M-file bvpexample.m. function yprime = bvpexample(t,y) %BVPEXAMPLE: Differential equation for boundary value %problem example. yprime=[y(2); -2*y(1)+3*y(2)]; Luego, escribimos las condiciones de contorno como el M-file bc.m, lo cual registra los residuos de contorno. function res=bc(y0,y1) %BC: Evaluates the residue of the boundary condition res=[y0(1);y1(1)-10]; Por residuos, nos referimos a la parte izquierda de la condición de frontera, una vez que ha sido puesta en 0. En este caso, la segunda condición de contorno es y(1) = 10, de modo que su residuo es y(1) - 10, que se registra en el segundo componente del vector que devuelve bc.m. Las variables y0 e y1 representan la solución en x=0 y en x=1, respectivamente, mientras que el 1 en el paréntesis indica el primer componente del vector. En el caso de que la segunda condición de contorno era y '(1)=10, reemplazaría a y1(1) - 10 con y1(2) - 10. Ahora estamos en condiciones de comenzar a resolver el problema de contorno. En el siguiente código, en primer lugar se especifica una cuadrícula de valores de x para resolver en una estimación inicial del vector que se daría para un problema de valor inicial [y (0), y '(0)]. (Por supuesto, y(0) conocido, pero y'(0) debe ser una conjetura. En términos generales, MATLAB va a resolver una familia de problemas de valores iniciales, buscando aquel para el cual las condiciones de contorno se cumplen.) Resolvemos el problema del valor límite con el resolutor built-in de MATLAB bvp4c. sol=bvpinit(linspace(0,1,25),[0 1]); sol=bvp4c(@bvpexample,@bc,sol); sol.x ans = Columns 1 through 9 0 0.0417 0.0833 0.1250 0.1667 0.2083 0.2500 0.2917 0.3333 Columns 10 through 18 0.3750 0.4167 0.4583 0.5000 0.5417 0.5833 0.6250 0.6667 0.7083 Columns 19 through 25

21 0.7500 0.7917 0.8333 0.8750 0.9167 0.9583 1.0000

sol.y ans = Columns 1 through 9 0 0.0950 0.2022 0.3230 0.4587 0.6108 0.7808 0.9706 1.1821 2.1410 2.4220 2.7315 3.0721 3.4467 3.8584 4.3106 4.8072 5.3521 Columns 10 through 18 1.4173 1.6787 1.9686 2.2899 2.6455 3.0386 3.4728 3.9521 4.4805 5.9497 6.6050 7.3230 8.1096 8.9710 9.9138 10.9455 12.0742 13.3084 Columns 19 through 25 5.0627 5.7037 6.4090 7.1845 8.0367 8.9726 9.9999 14.6578 16.1327 17.7443 19.5049 21.4277 23.5274 25.8196

Observamos que en este caso, MATLAB devuelve la solución como una estructura cuyo primer componente de la estructura sol.x simplemente contiene sólo los valores de x especificados. El segundo es sol.y, que es una matriz que contiene como primera fila los valores de y(x) en los puntos de la grilla x especificada, y como segunda fila los valores correspondientes de y '(x). Métodos Numéricos A pesar de que se pueden resolver EDOs en MATLAB sin ningún conocimiento de los métodos numéricos que emplea, a menudo es útil para comprender los principios básicos subyacentes. En esta sección se utilizará el teorema de Taylor para obtener métodos para aproximar la solución de una ecuación diferencial. Método de Euler Consideremos la ecuación diferencial general de primer orden

y supongamos que queremos resolver esta ecuación en el intervalo de valores de x [x0, xn]. Nuestro objetivo aspira a aproximar el valor de la solución y(x) en cada uno de los valores de x en una partición P = [x0, x1, x2, ..., xn]. Puesto que y(x0) está dado, el primer valor que necesitamos estimar es y(x1). Por el Teorema de Taylor, podemos escribir

donde c ε (x0, x0). Observando desde nuestra ecuación que y′(x0) = f(x0, y(x0)), tenemos

Si nuestra partición P tiene pequeños subintervalos, entonces x1 − x0 será pequeño y se puede considerar la cantidad pequeña

como un término error. Esto es, tenemos Facultad de Ingeniería – Universidad de Mendoza

Dr. Ing. Jesús Rubén Azor Montoya

22 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS”

(6.2) Ahora podemos calcular y(x2) de una manera similar usando el teorema de Taylor para escribir

Una vez más, tenemos desde nuestra ecuación que y'(x1) = f(x1, y(x1)), y así

Si eliminamos el término

como un error, entonces tenemos

donde el valor y(x1) requerido aquí puede ser aproximado por el valor de (6.2). Más generalmente, para cualquier k = 1, 2, ..., n - 1 se puede aproximar y(xk+1) de la relación

donde y(xk) se conocerán a partir del cálculo anterior. Al igual que con los métodos numéricos de la integración, es habitual en la práctica de tomar la partición en subintervalos de igual anchura,

(En el estudio de métodos numéricos para ecuaciones diferenciales, esta cantidad es a menudo denotada con h) En este caso, tenemos la relación general

Si decimos que los valores de y0, y1, ..., yn denotan nuestras aproximaciones para y en los puntos x0, x1, ..., xn (es decir, y0 = y(x0), y1 ≈ y(x1), etc), entonces podemos aproximar y(x) sobre la partición de P calculando iterativamente (6.3) Ejemplo: Utilice el método de Euler (6.3) con n = 10 para resolver la ecuación diferencial

en el intervalo [0, 1]. Llevaremos a cabo las primeras iteraciones en detalle, y a continuación vamos a escribir un archivo de MATLAB M-file para aplicar el

23

método en su totalidad. En primer lugar, el valor inicial y(0) = π nos da los valores x0 = 0 e y0=0. Si nuestra partición se compone de subintervalos de igual anchura, entonces x1 =∆x = 1/10 = 0.1, y de acuerdo con (6.3)

Ahora tenemos el punto (x1, y1) = (.1,π), y podemos utilizar esto y (6.3) para calcular

Ahora tenemos (x2, y2) = (0.2, 3.1725), y podemos utilizar esto para calcular

Más generalmente, podemos usar el M-file euler1.m function [xvalues, yvalues] = euler1(f,x0,xn,y0,n) %EULER: MATLAB function M-file that solve the %ODE y’=f, y(x0)=y0 on [x0,y0] using a partition %with n equally spaced subintervals dx = (xn-x0)/n; x(1) = x0; y(1) = y0; for k=1:n x(k+1)=x(k) + dx; y(k+1)= y(k) + f(x(k),y(k))*dx; end xvalues = x'; yvalues = y'; Podemos implementar este archivo con el siguiente código, que crea la figura 6.1. f=inline('sin(x*y)'); [x,y]=euler1(f,0,1,pi,10); plot(x,y)

Facultad de Ingeniería – Universidad de Mendoza

Dr. Ing. Jesús Rubén Azor Montoya

24 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS”

Con ∆x=0.01 [x,y]=euler1(f,0,1,pi,100); plot(x,y)

Resolutores ODE Avanzados Además de los resolutores ODE ode23 y ode45, ambos basados en el método de Runge-Kutta, MATLAB tiene otros resolutores adicionales, que aparecen a continuación junto con las sugerencias para su uso que da el MATLAB-help. • Resolutores Multipaso - ode113. Si se utilizan tolerancias estrictas de error o la resolución de un archivo ODE computacionalmente intensivo. • Problemas stiff (véase más adelante) - ode15s. Si ode45 es lento porque el problema es stiff. - ode23s. Si se está utilizando tolerancias de error crudas para resolver sistemas stiff y matriz de masa es constante. - ode23t. Si el problema sólo es moderadamente stiff y necesitas una solución sin amortiguación numérica. - ode23tb. Si se utilizan tolerancias de error crudas para resolver sistemas stiff. EDs stiffrígido Con EDs stiff nos referimos a una ED para la cual los errores numéricos crecen dramáticamente en el tiempo. Por ejemplo, considere la ED ordinaria y '= -100y + 100 t + 1,

y(0) = 1.

25

Ya que la variable dependiente, y, en la ecuación está multiplicada por 100, pequeños errores en nuestra aproximación tenderán a magnificarse. En general, debemos tomar considerablemente pasos más pequeños en tiempo para resolver EDs stiff, y esto puede alargar dramáticamente el tiempo de resolución. A menudo, las soluciones se puede calcular de manera más eficiente usando uno de los resolutores diseñado para problemas stiff. Ecuaciones diferenciales en derivadas parciales (PDE) en el espacio unidimensional Para las PDE de valores inicial-borde con tiempo t y una variable espacial x, MATLAB tiene un resolutor built-in llamado pdepe. Ecuaciones únicas Ejemplo. Supongamos, por ejemplo, que queremos resolver la ecuación del calor

MATLAB specifica tal PDE parabólica en la forma

Las formas standard de una ecuación diferencial El paso más importante en “preparación” de una ecuación diferencial para los resolutores de Mathcad es adquirirla en la forma Standard que entiendan los resolutotes. El proceso es tomar la ecuación diferencial y liberarse de cualquier derivada de orden más alto que aparezca, dejando sólo primeras derivadas. Su ecuación diferencial es entonces un sistema de ecuaciones diferenciales de primer orden. Por ejemplo, la ecuación diferencial: d2 2

dx

y( x) + 3 ⋅

d y( x) − 7 ⋅y( x) dx

4 ⋅x

contiene una segunda derivada la cual puede ser escrita como una primera derivada:

Facultad de Ingeniería – Universidad de Mendoza

Dr. Ing. Jesús Rubén Azor Montoya

26 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS”

d2 2

y( x)

dx

 d d  y( x)  dx  dx 

define dos funciones como: y0 ( x)

y( x)

y1 ( x)

d y0 ( x) dx

Ahora la ecuación diferencial contiene dos funciones y es esencialmente un sistema de dos ecuaciones diferenciales. d y0 ( x) dx d dx

y1 ( x)

y1 ( x) 4 ⋅x + 7 ⋅y0 ( x) − 3 ⋅y1 ( x)

El próximo paso es capturar esta información en una función vector único valuado D para usar con los resolutores de Ecuación Diferencial de Mathcad:



   4 ⋅x + 7 ⋅y0 − 3 ⋅y1 

DY( x , y) := 

y1

Sistemas de ecuaciones diferenciales ordinarias rkfixed(y, x1, x2, npoints, D) Retorna una matriz en la cual (1) la primera columna contiene los puntos en los cuales es evaluada la solución y (2) las columnas remanente contiene los valores correspondientes de la solución y sus primeras n-1 derivadas. Argumentos:
y debe ser o bien un vector de n valores inicales u un valor inicial único.
x1, x2 son puntos extremos del intervalo sobre el cual la solución a las ecuaciones diferenciales será evaluada. Los valores iniciales en y son los valores en x1.
npoints es el número de puntos más allá del punto inicial para el cual la solución es aproximada. Esto controla el número de filas (1 + npoints) en la matriz retornada por rkfixed.

27


D es una función vector-valued n-element conteniendo las primeras derivadas de las funciones desconocidas. Notas: •

rkfixed usa el método Runge-Kutta method de cuarto orden para resolver una ecuación diferencial de primer orden.

•

Se puede usar rkfixed para resolver una ecuación diferencial así también como un sistema de ecuaciones diferenciales.

Uso de la función rkfixed Para sistemas de ecuaciones diferenciales o por una que no es lineal en el término derivada de más alto orden, use rkfixed. El ejemplo de abajo muestra cómo usar rkfixed para evaluar la solución de una ecuación diferencial de segundo orden a puntos igualmente espaciados. Solve y( 0)

y' + 2 ⋅y

y'' 1

y'( 0)

1 y :=   3

3