Ecuaciones de Segundo Orden

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) f ( 0) t 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) y'( 0) 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 <- define las condiciones iniciales     −y1 + 2 ⋅y0  y1 D ( t , y) :=  <- primera derivada <- segunda derivada Z := rkfixed( y , 0 , 0.5 , 400 , D) 0 0 2 1 3 1 1.25·10 -3 1.004 2.999 2.5·10 -3 1.007 2.998 3 3.75·10 -3 1.011 2.996 5·10 -3 1.015 2.995 2 4 5 6.25·10 -3 1.019 2.994 7.5·10 -3 1.022 2.993 7 8.75·10 -3 1.026 2.992 8 0.01 1.03 2.99 6 Z= 1 0 <- Evalua solución en 400 puntos entre 0 y 0.5 9 0.011 1.034 2.989 10 0.013 1.037 2.988 11 0.014 1.041 2.987 12 0.015 1.045 2.986 13 0.016 1.049 2.985 14 0.018 1.052 2.984 15 0.019 1.056 2.982 Facultad de Ingeniería – Universidad de Mendoza Dr. Ing. Jesús Rubén Azor Montoya 28 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS” Solución usando tamaños de pasos adaptivos Rkadapt(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 o un escalar. 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. D es una función vector-valued n-element conteniendo las primeras derivadas de las funciones desconocidas. Notes: • A diferencia de rkfixed la cual integra en pasos de igual tamaño para alcanzar una solución, Rkadapt examina cuán rápido una solución está cambiando y adapta su tamaño de paso acordemente. Rkadapt usará tamaños de paso no uniformes unternamente cuando resuelve la ecuación diferencial, pero retornará la solución en puntos igualmente espaciados. Sistemas Alisados (Smooth systems) Bulstoer(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 o un escalar. 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 Bulstoer. D es una función vector-valued n-element conteniendo las primeras derivadas de las funciones desconocidas. 29 Notas: • Bulstoer usa el método Bulirsch-Stoer el cual sera ligeramente más preciso que el método Runge-Kutta usado por rkfixed. • Use la función Bulstoer en vez de rkfixed cuando sepa que la solución es alisada (smooth) Ejemplo: Movimiento del Péndulo Simple Se analizará mediante la función Bulstoer, también para evaluar una integral para el período del péndulo. Considérese un péndulo simple en movimiento en un plano sin resistencia del aire ni fricción. Sea L la longitud de la varilla soporte, g la aceleración de la gravedad y  el ángulo que forma el alambre y la dirección de la gravedad. La varilla soporte es rígida y sin masa. La ecuación diferencial para el ángulo (t) es: L⋅ d 2 2 θ ( t) + g ⋅ sin( θ ( t) ) θ ( 0) 0 θ0 dt Valores de los parámetros de entrada: L := 3 ⋅m θ 0 := π 4 X1     D( t , X) := g  − ⋅sin( X0)   L  g = 9.807 2 s T := 6 Punto derecho extremo N := 100 Número de pasos de tiempo   θ 0   , 0 , T , N , D  0   P := Bulstoer m n := 0 , 1 .. N Facultad de Ingeniería – Universidad de Mendoza Dr. Ing. Jesús Rubén Azor Montoya 30 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS” π 1.384 4 (P〈1〉 )n (P〈2〉 )n − π 4 − 1.384 0 1 2 〈 0〉 P n ( ) 0 3 4 5 6 6 theta derivative of theta El periodo del péndulo (tiempo requerido para una oscilación completa) es una integral elíptica completa sel primer tipo:  θ 0   2  x := sin π 4⋅ ⌠2 L  ⋅ g   ⌡0 1 2 1−x dφ = 3.614 sec 2 ⋅ sin(φ ) Sistemas Stiff Stiffb(y, x1, x2, npoints, D, J) Stiffr(Y, x1, x2, npoints, D, J) Cada una de estas funciones 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 remanentes contienen 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 o un escalar. 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 by Stiffb or Stiffr. 31 D es una función vector-valued n-element conteniendo las primeras derivadas de las funciones desconocidas. J es una función que retorna la matriz n por (n+1) cuya primera columna contiene derivadas y cuyas filas remanentes forman la matriz Jacobiana para el sistema de ecuaciones diferenciales. Notes: • Stiffb usa el método Bulirsch-Stoer para sistemas stiff. Stiffr usa el método Rosenbrock. • Use uno de los dos resolutores de ecuación diferencial específicamente diseñados para sistemas stiff: Stiffb y Stiffr cuando resuelva un sistema stiff. Stiff Differential Equations A pesar de su éxito, los métodos de Runge-Kutta también tienen sus fallas. Una clase particular de ecuaciones diferenciales problemáticas son aquellas que exhiben stiffness. No hay una definición universalmente satisfactoria de stiffness. En la práctica se aprenderá a sospechar que una ODE (ecuación diferencial ordinaria) es stiff cuando un resolutor ODE falla para producir una solución eficiente. Las sospechas serán confirmadas si verdaderamente se obtienen mejores resultados cuando se usa un método diseñado para problemas stiff, tales como las rutinas de Matlab ode23s u ode15s. Ejemplo de ODE stiff Un ejemplo de modelo matemático, más que una aplicación práctica es: ( Y' ( t) −Y( t) − K ⋅ Y ( t) − e Y ( 0) 1+ δ −t ) K> > 1 aquí K es una constante grande, digamos 10 6 y δ es un número dado tamaño ordinario, digamos δ = 1. La solución se conoce que es: Y( t) − ( K+1) ⋅ t δ ⋅e −t +e En esta formula el primer término se desvanece rápidamente a medida que la ecuación diferencial fuerza la solución sobre la curva e -t . − ( K+1) ⋅ t Y ( t) := δ ⋅e −7 t := 10 Facultad de Ingeniería – Universidad de Mendoza −7 , 10 −7 + 10 .. 0.001 Dr. Ing. Jesús Rubén Azor Montoya 32 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS” 1 0.905 0.5 Y ( t) 0 0 2 .10 −7 1×10 6 4 .10 6 6 .10 t 6 8 .10 6 −5 10 Después esta solución cambia lentamente. − ( K+1) ⋅t Y ( t) := δ ⋅e 0.99 Y ( t) −t t := 0.01 , 0.02 .. 5 +e 1 0.5 −3 6.738×10 0 1 0.01 2 3 t 4 5 5 La región corta de cambio rápido es la capa frontera, la curva de variación lenta es conocida como la solución estado estacionario. Este comportamiento es típico de problemas stiff. Cuando se desea sólo el ultimo punto Si no se está interesado en ver todos los valores intermedios tomados por una solución, se puede ahorrar tiempo de procesamiento y reducir consumo de memoria usando una de las funciones de abajo. stiffb y stiffr para sistemas stiff. rkadapt para sistemas lentamente variables. bulstoer para sistemas alisados (smooth) A diferencia de sus equivalentes en mayúsculas, no se debe especificar el número de puntos en el resultado. En su lugar, debe especificar el nivel de exactitud y permitir a estas funciones generar una solución que tenga el número de elementos requeridos. Ejemplo: 33 Solve y'' y     −y   para x < 0 donde y(-1)=1 e y(1)=2 para x > 0    ( x < 0) ⋅y0 + ( x ≥ 0) ⋅( −y0)  y1 D( x , y) :=  v10 := 1 < - valor de intento para y(-1 v20 := 1 < - valor de intento para y(1)  1    v10  load1( x1 , v1) :=   2    v20  load2( x2 , v2) :=  score( xf , y) := y xf := 0 <- punto de dictontinuidad < - y(-1) < - valor de intento para y(-1 < - y(1) < - valor de intento para y(1) < - indica a Mathcad igualar las dos mitades de solución en x=xf S := bvalfitv1 ( , v2 , −1 , 1 , 0 , D , load1, load2, score) S = ( 0.092 −0.678 ) < - contiene [ y(-1) y(1) Uso de valores intermedios para derivar condiciones iniciales bvalfit(v1, v2, x1, x2, xf, D, load1, load2, score) Retorna un vector conteniendo aquellos valores iniciales dejados sin especificar en x1. Argumentos: v1 es un vector de intento para valores iniciales dejados sin especificar en x1. v2 es un vector que contiene intentos para valores iniciales dejados de especificar en x2. x1, x2 son puntos terminales del intervalo sobre la cual la solución a ecuaciones diferenciales serán evaluadas. xf es un punto entre x1 y x2 en los cuales la trayectorias de las soluciones comienzan en x1 y aquelllas que comienzan en x2 están constreñidas a ser igual. D es una función vector-valuada de n elementos conteniendo las primeras derivadas de las funciones desconocidads. load1 es una function vector-valuada cuyos n elementos corresponden a los valores de las n funciones incógnitas en x1. Facultad de Ingeniería – Universidad de Mendoza Dr. Ing. Jesús Rubén Azor Montoya 34 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS” Algunos de estos valores serán constantes especificadas por sus condiciones iniciales. Si un valor es desconocido se deberá usar el valor de intento correspondiente desde v1. load2 es análoga a load1 pero los valores tomados las n funciones incógnitas en x2. score es una función vector-valuada de n elementos usada para especificar cómo se desea que las soluciones igualen a xf. Usualmente se deseará definir score(xf, y) := y para hacer que las soluciones para todas las funciones incógnitas coincidan en xf. Notas: • bvalfit resuelve un problema de valores límites de dos puntos “disparando” desde los puntos terminales y equiparando trayectorias de la solución y sus derivadas en el punto intermedio. Este método se convierte en especialmente útil cuando la derivada tiene una discontinuidad en alguna parte del intervalo de integración. sbval (v1, v2, x1, x2, D, load, score) Retorna un vector conteniendo aquellos valores iniciales dejados sin especificar en x1. Argumentos: v1 es un vector de intento para valores iniciales dejados sin especificar en x1. v2 es un vector que contiene intentos para valores iniciales dejados de especificar en x2. x1, x2 son puntos terminales del intervalo sobre la cual la solución a ecuaciones diferenciales serán evaluadas. D es una function vector-valuada de n elementos conteniendo las primeras derivadas de las funciones desconocidads. load es una función vector-valuada cuyos n elementos corresponden a los valores de las n funciones incógnitas en x1. Algunos de estos valores serán constantes especificadas por sus condiciones iniciales. Si un valor es desconocido se deberá usar el valor de intento correspondiente desde v1. score es una function vector-valuada que tiene el mismo número de que v. Cada elemento es la diferencia entre una condición inicial en x2, como originalmente se especificó, y el estimado correspondiente desde la solución. El vector store mide cuán ajustadamente una solución propuesta se adapta a las condiciones iniciales en x2. un valor de 0 para cualquier elemento indica una adaptación perfecta. Notas: • sbval resuelve un problema de valores límites de dos puntos “disparando” desde los puntos terminales y equiparando trayectorias de la solución y sus derivadas en el punto intermedio. 35 • sbval computa los valores iniciales que la solución debe tener para ajustar al valor inicial que se especificó. Se deben entonces tomar los valores iniciales retornados por sbval y resolver el problema de valor inicial resultante. • Cuando se conozca la solución en x1 y en x2, y se tengan n valores conocidos, se deberá usar sbval para evaluar los valores iniciales perdidos en x1. Una vez que se tienen estos valores iniciales perdidods, se tendrá un problema con valores iniciales en vez de un problema de valores límites de dos puntos. Ecuación de Poisson relax(a, b, c, d, e, f, u, rjac) Retorna una matriz cuadrada en la cual: (1) una locación de elemento en la matriz corresponde a su locación dentro de la región cuadrada, y (2) su valor Aproxima el valor de la solución en este punto. Argumentos: a, b, c, d, e son matrices cuadradas todas del mismo tamaño conteniendo los coeficientes de la ecuación diferencial. f es una matriz cuadrada conteniendo el término fuente en cada punto dentro del cuadrado. u es una matriz cuadrada conteniendo valores límites junto con los flancos de la región y los intentos iniciales para la solución dentro de la región. rjac es un radio spectral de la iteración Jacobi. Esto controla la convergencia del algoritmo de relajación. Puede oscilar desde 0 a 1 pero su valor óptimo depende de los detalles del problema. Notas: • Esta function usa el método de relajación para converger a la solución. • Se deberá usa la función relax si se conoce el valor tomado por la function desconocida u(x, y) sobre el total de los cuatro lados de una región cuadrada. • Si la condición de contorno es cero sobre el total de los cuatro lados del cuadrado dominio de integración, use en su lugar la función multigrid. Facultad de Ingeniería – Universidad de Mendoza Dr. Ing. Jesús Rubén Azor Montoya 36 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS” Uso de la función relax M R R := 1 R := 32 3⋅ MR , R := −2 , 4 4 M R R := 2 8 i := 0 .. 32 ai , k := 1 , 2 k := 0 .. 32 b := a vi , k := 0 d := a c := a e := −4 ⋅a f := −M S1 := relax( a , b , c , d , e , f , v , 0.95) S1 Ecuación de Poisson con condiciones de contorno cero multigrid(M, ncycle) Retorna una matriz cuadrada en la cual: (1) una locación de elemento en la matriz corresponde a su locación dentro de la región cuadrada, y (2) su valor aproxima el valor de la solución en este punto. Argumentos: M es una matriz cuadrada de 1+2^n filas cuyos elementos corresponden al término fuente en el correspondiente punto en el dominio cuadrado. 37 ncycle es el número de ciclos en cada nivel de la iteración multigrid. Un valor de 2 dará generalmente una Buena aproximación de la solución. Notas: Si u(x, y) es cero sobre los cuatro lados del cuadrado, use la función multigrid. Esta function frecuentemente resolverá el problema más rápido que relax. • Uso de multigrid: R := 32 MR , R := 0 M R R := 1 3⋅ M R R := 2 M R 3 ⋅R := −2 4 , , 4 4 8 4 , 2 S := multigridM ( , 2) S Facultad de Ingeniería – Universidad de Mendoza Dr. Ing. Jesús Rubén Azor Montoya 38 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS” $$$$$$$$$$$$$$$$$$$$$$$$ MATLAB contiene dos funciones para calcular soluciones numérica de ecuaciones diferenciales ordinarias; "ode23" y "ode45". A continuación se describen los argumentos de los comandos MATLAB ODE. [x,y] = ode23('función',a,b,inicial) Esta instrucción regresa un conjunto de coordenadas "x" y "y" que representan a la función y=f(x), los valores se calculan a través de métodos Runge-Kuta de segundo y tercer orden. El nombre "función", define una función que representa a una ecuación diferencial ordinaria, ODE23 proporciona los valores de la ecuación diferencial y'=g(x,y). Los valores "a" y "b" especifican los extremos del intervalo en el cual se desea evaluar a la función y=f(x). El valor inicial y = f(a) especifica el valor de la función en el extremo izquierdo del intervalo [a,b]. [x,y] = ode45('función',a,b,inicial) Esta instrucción regresa un conjunto de coordenadas "x" y "y" que representan a la función y=f(x), los valores se calculan a través de métodos Runge-Kuta de cuarto y quinto orden. El nombre "función", define una función que representa a una ecuación diferencial ordinaria, ODE45 proporciona los valores de la ecuación diferencial y'=g(x,y). Los valores "a" y "b" especifican los extremos del intervalo en el cual se desea evaluar a la función y=f(x). El valor inicial y = f(a) especifica el valor de la función en el extremo izquierdo del intervalo [a,b]. Las instrucciones "ODE23" y "ODE45" contienen dos parámetros adicionales. Se usa un quinto parámetro para especificar una tolerancia relacionada con el tamaño del paso; las tolerancias por omisión son 0.001 para ODE23 y 0.000001 para ODE45. Existe un sexto parámetro que sirve para solicitar que la función exhiba resultados intermedios, es decir, que realice rastreo; el valor por omisión "0" indica que no se desean rastrear los resultados. Como ilustración de la función ODE de MATLAB, se presentan los pasos para calcular soluciones numéricas de ecuaciones diferenciales, las siguientes instrucciones MATLAB definen las funciones requeridas para evaluar la ecuación diferencial deseada. Se graban las siguientes instrucciones con su editor ASCII favorito, en lo particular yo uso el "editeur", si Ud. desea usarlo también, esta disponible en la siguiente URL; http://proton.ucting.udg.mx/shareware/editeur/editeur.zip 39 function dy = g1(x,y) % % g1 % esta función evalúa una ODE % ecuación diferencial de primer grado % dy = 3*x.^2; El siguiente paso consiste en grabar este archivo como "g1.m", sobre algún subdirectorio de trabajo valido para el MATLAB, las siguientes instrucciones resuelven g(x,y) dentro del intervalo [2,4] con condición inicial 0.5 para y=f(2). % Determinar la Solución de la EDO % % dy = 3*x.^2; % [t,y] = ode23('g1',[2,4],0.5); plot(t,y,'o'),... title('Solución de la Ecuación dy = 3*x.^2'),... xlabel('Tiempo'),ylabel('y = f(t)'),grid Sobre el subdirectorio de trabajo valido se graba este archivo como "mat1.m" y se escribe mat1, generándose la siguiente solución gráfica. Facultad de Ingeniería – Universidad de Mendoza Dr. Ing. Jesús Rubén Azor Montoya 40 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS” Mathcad.pdf Las dos funciones más fáciles de usar, Odesolve y Pdesolve, son usadas dentro de Solve Blocks. Hay muchas otras funciones que ayudan a relover sistemas de ecuaciones diferenciales, incluyendo Adams, AdamsBDF, BDF, Bulstoes, bvalfit, Jacob, multigrid, numol, Radau, relax, Rkadapt, rkfixed, sbval, statespace, Stiffb y Stiffr. La función Odesolve es usada en un Solve Block para resolver una única ecuación diferencial (ED) o un sistema de EDs. Retorna la solución como una función de la variable independiente. Hace esto salvando soluciones en un número específico de puntos (npoints) igualmente espaciados en el intervalo de la solución, y luego interpolando entre estos puntos usando la función lspline. El Solve Block pone el valor inicial o constreñimientos de límites. La ODE debe ser lineal en su término derivado más alto, y el número de condiciones iniciales o en bordes debe ser igual al orden (u órdenes) de la(s) ODE(s). La función tiene la forma Odesolve([vector],x,b,[npoints]). Los argumentos son como sigue: • Vector es usado sólo para sistemas de ODEs y es un vector de nombres de función (con ningún nombre de variable incluido) como ellas aparecen dentro del Solve Block. • x es el nombre de la variable de integración. • b es el punto final del intervalo de la solución. (El punto inicial del intervalo de la solución es especificado por las condiciones iniciales). 41 • npoints es opcional y es el número entero de puntos equiespaciados usados para interpolar la función solución. El valor por default es de 1000. Si npoints es incrementado, luego la solución interpolada es más exacta. El valor default es usualmente adecuado, pero si se desea resolver sobre un gran intervalo, se debe poner a npoints a un valor mayor a 1000. El incremento de npoints lleva ala incremento en el tiempo de cálculo. ----- dsolve - Symbolic solution of ordinary differential equations Syntax dsolve('eq1','eq2',...,'cond1','cond2',...,'v') dsolve(...,'IgnoreAnalyticConstraints',value) Descripción dsolve('eq1','eq2',...,'cond1','cond2',...,'v') resuelve simbólicamente las EDs ordinarias eq1, eq2,... usando v como la variable independiente. Aquí cond1,cond2,... especifica condiciones iniciales o en bordes o ambas. También se usa la siguiente sintaxis: dsolve('eq1, eq2',...,'cond1,cond2',...,'v'). La variable independiente por default es t. The letter D denotes differentiation with respect to the independent variable. The primary default is d/dx. The letter D followed by a digit denotes repeated differentiation. For example, D2 is d2/dx2. Any character immediately following a differentiation operator is a dependent variable. For example, D3y denotes the third derivative of y(x) or y(t). You can specify initial and boundary conditions by equations like y(a) = b or Dy(a) = b, where y is a dependent variable and a and b are constants. If the number of the specified initial conditions is less than the number of dependent variables, the resulting solutions contain the arbitrary constants C1, C2,.... Ejemplo: >> dsolve('D2y=-3*y','y(0)=1') Facultad de Ingeniería – Universidad de Mendoza Dr. Ing. Jesús Rubén Azor Montoya 42 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS” ans = C1*sin(3^(1/2)*t)+cos(3^(1/2)*t) You can input each equation or a condition as a separate symbolic equation. The dsolve command accepts up to 12 input arguments. dsolve puede producir los siguientes tres tipos de salida: * Para una ecuación y una salida, dsolve retorna la solución resultante con soluciones múltiples para una ecuación no-lineal en un vector simbólico. * Para varias ecuaciones y un número igual de salidas, dsolve ordena los resultados alfabéticamente y los asigna a las salidas. * Para varias ecuaciones y una única salida, dsolve retorna una estructura conteniendo las soluciones. Si dsolve no puede encontrar una solución closed-form (explícita), intenta buscar una solución implícita. Cuando dsolve retorna una solución implícita, indica una advertencia. Si dsolve no puede encontrar una solución ni implícita ni explícita, da una advertencia y retorna el sym vacío. En tal caso, se puede encontrar una solución numérica, usando las funciones MATLAB ode23 u ode45. En algunos casos involucrando ecuaciones no-lineales, la salida es una ecuación diferencial equivalente de más bajo orden o una integral. dsolve(...,'IgnoreAnalyticConstraints',value) acepta los siguientes valores: * value = 'all' aplica simplificaciones puramente algebraicas a las expresiones sobre ambos lados de las ecuaciones. Estas simplificaciones pueden no ser generalmente válidas. El valor default de esta opción es all. * value = 'none' resuleve EDs ordinarias sin suposiciones adicionales. Los resultados obtenidos con esta opción son correctos para todos los valores de los argumentos. Nota: Por default, el solver no garantiza la corrección y completitud de los resultados. Si no se pone la opción IgnoreAnalyticConstraints a none,siempre verifique los resultados retornados por el comando dsolve. Examples Solving Ordinary Differential Equations Symbolically dsolve('Dx = -a*x') ans = C2/exp(a*t) Specifying the Dependent Variable 43 The following differential equation presents f as a dependent variable: dsolve('Df = f + sin(t)') ans = C4*exp(t) - sin(t)/2 - cos(t)/2 Specifying the Independent Variable dsolve('(Dy)^2 + y^2 = 1','s') ans = 1 -1 cosh(C7 + s*i) cosh(C11 - s*i) Setting Initial and Boundary Conditions dsolve('Dy = a*y', 'y(0) = b') ans = b*exp(a*t) dsolve('D2y = -a^2*y', 'y(0) = 1', 'Dy(pi/a) = 0') ans = (1/exp(a*t*i))/2 + exp(a*t*i)/2 Solving a System of Differential Equations z = dsolve('Dx = y', 'Dy = -x') z = x: [1x1 sym] y: [1x1 sym] Enter z.x and z.y to see the results: z.x ans = C20*cos(t) + C19*sin(t) z.y ans = C19*cos(t) - C20*sin(t) Using the IgnoreAnalyticConstraints Option By default, the solver applies the set of purely algebraic simplifications that are not correct in general, but that can result in simple and practical solutions: Facultad de Ingeniería – Universidad de Mendoza Dr. Ing. Jesús Rubén Azor Montoya 44 “HERRAMIENTAS COMPUTACIONALES EN CIENCIAS EXACTAS” y = dsolve('Dy=1+y^2','y(0)=1') y = tan(pi/4 + t) To obtain complete and generally correct solutions, set the value of the option IgnoreAnalyticConstraints to none: y = dsolve('Dy=1+y^2','y(0)=1',... 'IgnoreAnalyticConstraints','none') y = piecewise([C29 in Z_, tan(pi/4 + t + pi*C29)]) The algebraic simplifications also allow you to obtain solutions for the equations that the solver cannot compute when it uses strict mathematical rules: dsolve('Dv=19.6*v-0.00196*v^2','v(0)=1') ans = 10000/(exp(log(9999) - (98*t)/5) + 1) versus dsolve('Dv=19.6*v-0.00196*v^2','v(0)=1',... 'IgnoreAnalyticConstraints','none') Warning: Possibly spurious solutions [solvelib::checkSolutions] ans = piecewise([C38 in Z_, 1/(exp(log(9999) - (98*t)/5 +... 2*pi*C38*i)/10000 + 1/10000)]) Diagnostics If dsolve cannot find an analytic solution for an equation, it prints the warning: Warning: Explicit solution could not be found. and returns an empty sym object. See Also