Ecuaciones Diferenciales Con MATLAB Cesar Perez Lopez

GMatl#Iforks.

ECUAGIONES DIFERENGIALES CON MATLAB

Ejemplos y ejercicios resueltos

CESAR PÉREZ LÓPEZ

INDICE MilTRODUCCIÓN PRÁCTICA A MATLAB

MT,TCTONES DIFERENCIALES DE PRIMER ORDEN. ECUACIONES m[,{cTAS, EN VARIABLES SEPARADAS, HOMOGENEAS Y LINEALES........25

1l EcuACIoNES DIFERENCIALES DE PRIMER oRDEN ]JECUACIONESENVARIABLESSEPARADAS..................

]j

............... EXACTAS...

ECUACIONES DIFERENCIALES HOMOGÉNEAS

:' 4 ECUACIONES DIFERENCIALES

1,i ECUACIONES DIFERENCIALES LINEAIES

]I

...............25 ...............28 ........... 30 ............ 33 .................... 35

TACIONES DIFERENCIALES DE ORDEN SUPEROR. TRANSFORMADA DE ...................37

I,ñPLACEYTIPOSESPECIALESDEECUACIONES...... 3.I ECUACIONES ORDINARIAS DE ORDEN

SUPERIOR.

.....,.,..,...,...37

Jl

ECUACIONES LINEALES DE ORDEN SUPERIORHOMOGÉNEAS EN ..................39 COEFICIENTES CONSTANTES.

]S ECUACIONES NO HOMOGÉNEAS CON COEFICIENTES CONSTANTES, \-ARTACIÓN DE

PARÁMETROS...........

....,.,41

J.+ ECUACIONES NO HOMOGÉNEAS CON COEFICIENTES VARIABLES.

CAUCHY-EULER.... 3*i TRANSFORMADA DE LAPLACE .................. -1.6 POLINOMTOS ORTOGONALES ECL]ACIONES DE

,.,....44 .............45 ............. 48

3.6.1 Polinomios de Chebychev de primera y segunda especie.....................'........... 48

Legendre................... 3.6.3 Polinomios asociados de Legendre 3.6.4 Polinomios de Hermite................... 3.6.2 Polinomios de

........................49 ..........................49 ..........................49

4

CON ECUACIONES DIFERENCIALES

MÑ 50

1Íiilll"-:i:::T'"r;::i:::'""": ""'1""""""'

so

""""""""""""" 50 """"""""""""' 50

Jacobi 3.6'7 Polinomios de Gegenbauer de 3.6.8 Polinomios DE AIRY Y 3.7 FLTNCIONES

""""""' BESSEL

.."...".

51

ECUACIONESDIFERENCIALESponITÉroDoSAPRoXIMADOS...................53

"""""53 4.lEcuAcIoNEs"o*Ñ"o*o*'**á*AI-lNo'un{EArEsYNoI.lN{FAtrs' """"""' 53 rr¡Étooosepnoxnr¡anos TAYLoR """"""" LAS sERIES DE """""' 56 4.2 Er- MÉroDo DE 4.3

Er-MÉroDoDERtrNGE-Kurre""""

DIFERENcIALEs srsrEMAs DE EcuAcroNEs TTAS"""""""""' DIFERENCIAS FTN

Y EcuAcIoi::.::................ se ........59

5.lSISTEMAS,"u**io*II.IEALES;OMoGENEASCoNCoEFICIENTES coN CoEFICIENTES CONSTA\TTES

OE*oáo*"^*r*ono*oo^t* CONSTANTES 53EcuAcloNEsÑ'*u**'FINnAs"""'

s2 STSTEMAS

........60

""""""""'

A LAs EcuAcroNrs EN MATLAB' APLIcAcroNEs NuuÉnTco o cÁl,cut """

MATLAB Y LA PRoGRAMactÓN 6.2 EotronoE rExro""""

6.1

6.3

"""""

61

65

"""""' 65 """""'67

ScuPrs"

"""""""'73 6.4Ft-lNctoNesyM-rtcrrenos.FLINCTI0N,EVALYFEVAL...'.....................68 6.5

Verueelss

LocALES Y

GLoBALEs

""""""""""'75

6.6TtposDEDAros 6.7

Cor¡rn.ot-t';;;;;'BucLEsFoR'wHILErIFELSEIF ....'..""" "' "'' ""''' ""2

INDICE

Euler....... 6.9.2 Método de Heun....... 6.9.3 Método de las series de Taylor...... 6.9.1 Método de

...............:.....

5

...... 86

......81 ...... 88

ECUACIONES EN DIFERENCIAS CON VALORES INICIALES, VALORES EN ....................... 95 LA FRONTERA Y EN DERIVADAS PARCIALES

7.I SOLUCIÓN NUMÉRICADEECUACIONES DIFERENCIALES .........................95 7.2 ECUACIoNES DIFERENCIALES oRDINARIAS CON VALORES INICIALES ......................... 95 7.3 EcUACToNES DTFERENCTALES oRDINARIAS coN vALoRES EN LA FRoNTERA.............. 99 7.4 EcUACtoNES DIFERENCTALES EN DERTvADAS

pARCIALES...

................. 102

Capítulo

1

INTRODUCCION PRACTICA A MATLAB CÁlCuT,O NUMÉRICO CON MATLAB Podemos usar Matlab como una computadora numérica de gran poEncia. La mayoría de las calculadoras manejan números sólo con un grdo de precisión prefijado, sin embargo Matlab realiza cálculos exactos Grn la precisión que sea necesaria. Además, a diferencia de las *¡ladoras, podemos realizar operaciones no sÓlo con nÚmeros iüvftJuales, sino también con objetos tales como las matrices. La mayoría de los temas del cálculo numérico clásico, son tratados ei¡ este software. Soporta cálculo matricial, estadística, interpolación, ajuste pr mínimos cuadrados, integración numérica, minimización de funciones, pgramación lineal, resolución numérica de ecuaciones algebraicas y frrenciales y un largo etcétera de procesos de análisis numérico que hrrrc viendo altratar los sucesivos temas de este libro. Veamos algunos ejemplos de cálculo numérico con Matlab. (Como ¡p sabemos, para obtener los resultados es necesario teclear Enter una vez ffiritas las expresiones correspondientes al lado del prompt ">>") 1) Podemos calcular sencillamente 4+3 y obtener como resultado 7.

r*+3

8

ECUACIONES DIFERENCIALES CON MATLAB

ans = 7

2) También podemos obtener el valor de 3 elevado a 100, sin haber fijado antes la precisión, basta para ello teclear 3^100. » 3^100

ans =

{

5.1538e+047

I

"j

{

3) También podemos utilizar el comando "format long e" para pasar resultado de la operación anterior a notación científica con 16 dígitos más el exponente.

tI

ú format long e »» 3^100

I

ans =

I

»»

5.1 537 7 5207 3201

1

5e +0

47

4) También podemos trabajar con números complejos. Obtendremos

el resultado de la operación (2+3i) elevado a 10, tecleando la expresión (2+3i)^10,

»

(2 + 3¡)^10

m

ans =

ü

-1.41 5249999999998e+005 -1 .456680000000000e+005i

5) El resultado anterior también se puede obtener en formato corto, utilizando el comando "format short". »»

format short

» (2 + 3¡)^10

ans = -1 .41

f,

52e+005- I .4567e+005¡

l

§ §

l

I

CAPÍTULo 1. INTRoDUCCIÓN PRÁCTICA A MATLAB 9

6) También podemos calcular el valor de la función de Bessel en el punto 11.5 . Para ello teclearemos Besselj(0,11.5). » Besselj(O, 11.5)

ans = 0.2150

7) También podemos realizar integrales numéricas. Para calcular la integral entre 0 y n de la función Seno(Seno(x)) teclearemos la expresión int(sin(sin('x')), 0, pi). » int(sin(sin('x')), 0, pi)

ans = 1235191

1

62052677 122517 9981 3685248*pi

Estos temas serán tratados más a fondo en sucesivos capítulos a lo largo del libro.

cÁr,cur.o SIMBóLICO coN MATLAB

b, ¡

I

Matlab maneja perfectamente el cálculo matemático simbólico, manipula con facilidad y rapidez las fórmulas y expresiones algebraicas y puede realizar la mayoría de operaciones con las mismas. Puede expandir, factorizar y simplificar polinomios y expresiones racionales y trigonométricas, puede encontrar soluciones algebraicas de ecuaciones polinómicas y sistemas de ecuaciones, puede evaluar derivadas e integrales simbólicamente y encontrar funciones solución de ecuaciones diferenciales, puede manipular series de potencias, límites y muchas otras facetas de la matemática algebraica.

r

: I

Para realizar esta tarea, Matlab requiere que todas las variables (o expresiones algebráicas) se escriban entre comillas simples. Cuando Matlab recibe una variable o expresión entre comillas, considera que es simbólica. Veamos algunos ejemplos de cálculo simbólico con Matlab.

ECUACIONES DIFERENCIALES CON MNTLqg

10

1)Podemoselevaralcubolasiguienteexpresiónalgebraica: 2. Esto se consigue tedeándo la siguiente alexpresiÓn: otra expresión gebraica: li-pJro tr(-* r ¡1x*Z¡xürzrg;l El iesultado será

(x+1)(x+2)- (x+2)

>»

expand('((x + 1).(x + 2)'(x + 2)12¡43',

ans =

-x^3-6*x^2'12"x'8

2)Podemosfactorizarelresultadodelcálculodelejemploanterior

+ tecleandá: factor('((x + 1)*(x + 2)-(x 2)"2)^3') »

factor('((¡ + 1)*(x + 2)'(x + 2)42¡n'',

ans = -(x+2)^3

3)PodemosresolverlaintegralindefinidadelafunciÓn

(x^2)Sen'o(x)^2 tecleando: nt('x^2*sin(x)^2"'x') i

»»

int('x^2*sin(x)^2','x')

ans = x^2*(-112*cos(x)*sin(x)+1t2*x|.1l2*x*cos(x)^2+1/4*cos(x)*sin(x)+1/4*x.

l/3*x^3 4) Podemos simplificar el resultado anterior: >» ,,

simplify(int('x^2*sin(x)^2','x'))

ans =

ii:l fl

-112*x^2*cos(x)*sin(x)+l/6*x^3-l/2*x*cos(x)^2+1l4"cos(x)*sin(x)+114*x

t,i

[l

rl rt

5)PodemosexpresarelresultadoanteriorconnotaciÓnmatemática más elegante:

,1 I

» pretty(simplify(int('x^2*sin(x)^2','x')»

fl;E

cnpíruLo 1. rNTRoDUcclóru pnÁcrlcn A MATLAB

11

a-

n:

7.12x cos(x) sin(x) +L/6x 1/4 cos(x) sin(x) + L/4 x

al-

-L/2 xcos(x)

6) Podemos desarrollar en serie de potencias de orden 12 la función )^2, presentando el resultado en forma elegante:

^2*sin(x)^2',12))

68 ,. -Llg x +2/45x I

fxr

7)

Podemos resolver

-L|3LSx

la

10

+

o(x

L2 )

ecuación 3ax-7x^2+x^3=0

(a, es

un

a*x-7*x^2+ x^3 = 0.,'x')

ETII

Ll2* (49-L2*al L/2* (49-L2*al

n n

0l (L/211

lL/2)

1

8) Podemos hallar las cinco soluciones de la ecuación x^5+2x+1=0: rS+2*x+1 ,,,x')

z

5+2*_2+11

Como el resultado obtenido no expresa explícitamente las cinco aplicamos el comando "allvalues":

solve('x^5+2*x+1','x'))

T Éüca

-=

I--?01873568855861 9- . 87 96 97 L97 9298240* i) ü:--7018735688558 61 9+ . 87 9697L97 9298240*il -.48638903s934s4301

-L

12

I I


. .

9450680868231334- . 8545175L44390459*il 94s0680868231334+. 8s45175t44390459*il

Por otra parte, Matlab puede utilizar las librerías del programa Maple V, para trabajar con matemática simbÓlica y poder extender así su campo de ,"óión. De eita forma, se puede usar Matlab para trabajar en temas como la las formas diferenciales, la geometría euclídea, la geometría proyectiva, estadística, etc.

A su vez, tambíén se pueden ampliar temas del cálculo numéri@, utilizano las librerías desde Matlab las librerías de Maple (combinatori4 optimización, teoría de números, etc')

MATLAB Y MAPLE

qtr Matlab se apoya en las librerías del programa Maple . Siempre Matlah' sea necesario acudir á cualquier comado o función de Maple desde se utiliza el comando "maplá" seguido de la sintaxis correspondiente en d posible si entorno Maple. Hay que constatai que esta disponibilidad solo es ;'Tóolbox" de cálculo simbólico "Extended Symbolic Mill se ha instaiado el Toolbox". Para utilizar un comando Maple desde Matlab, Ia sintaxis es

h

siguiente:

maple('si ntaxis-del-comando-en-entorno-Maple') o también: m dpl

e's i ntaxis-d e l-co man d o-en-e nto

rn

o-M

ap

I

e'

Para utilizar una funciÓn Maple con N argumentos desde Matlab sintaxis es la siguiente:

maple('sintaxis-de-la-función-en-entorno-Maple', argumentol, argumento2, ..., argumentoN) Veamos algunos ejemplos al respecto. 1) Podemos calcular el límite de la función (x^3-1)/(x-1) cuando »»

maple('limit((x^3'1 )/(x'1),x=l )')

,--

b

cApÍruLo 1. tNTRoDUccrórrr pRÁcrtca A MATLAB

13

ans = 3

e

e También podría haberse utilizado la siguiente sintaxis:

0

h o,

r maple'limit((x^3-l )/(x-1 ),x=1 )'; ans =

a,

3

2) Podemos calcular el máximo común divisor entre 10.000 y 5.000: ue ¡b,

el

rsi rth

la

r maple('gcd', 10000, 5000) ¿lfls =

5{n0

CNÁTTCOS CON MATLAB Matlab produce gráficos de dos y tres dimensiones, así como mntornos y gráficos de densidad. Se pueden representar los gráficos y listar bs datos, permite el control de colores, sombreados y otras características de los gráficos, también soporta gráficos animados. Los gráficos producidos por Matlab son portables a otros programas. Veamos algunos ejemplos de gráficos con Matlab

l, la

't) Podemos representar

efitre

-xl4 y

ru/4,

la función xSeno(1/x) para x variando tomando 300 puntos equidistantes del intervalo. Ver figura

1.1

r r=linspace(-pi/4,pi/4,300) » y=x.r'sin(1./x)' » plot(x,y) ->1:

;

74


Figura 1.1

2) Podemos dotar al gráfico anterior de las opciones de marco y rejilla, así como de título para el propio gráfico y etiquetas para sus ejes. Ver figura 1.2. » x=l i nspa ce(-pil 4,pil 4,300) » y=x.*s¡n(1./x);

;

» plot(x,y); » grid;

xlabel('Variable independiente X'); ylabel('Variables dependientes Y, Z'); » title('Funciones Seno y Coseno sobre los mismos ejes') »

»»

FuncioFe6

o.E

EEno y Coséño sohre los mismas

ejeE

ü-E

É

ñ4

-q

E n,:

€ü g

-ü.2 -ú,4 r0

E

A.6

.il.{-

-A_2

ú

tr_:

Varisble independaeñIe X

[.4

0.5

O.E

Figura 1.2

3)

Podemos generar

un gráfico de

superficie paru

la

z=Sen(sqrt(yt')+y^2)lsqrt(x^2+y^2), haciendo variar x e y en el inte-. = valores (-7.5,7.5), tomando puntos equidistantes en 5 décimas . Ver flq--.

CAPíTULo 1. INTRoDUCCIÓN PRÁCTICA A MATLAB 15

» x=-7.5:.5:7.5;

» Y=x;

[X,Y]=¡¡ss h g rid(x,y) ; Z=sin(sqrt(X.^2+Y.^2))./sqrt(x. ^2+Y .^21:' »» surf(X,Y,Z) »» »»

v !r

Figura 1.3 Estos gráficos tridimensionales permiten a la perfección hacerse una idea de las figuras en el espacio, y son muy útiles a la hora de determinar visualmente intersecciones entre distintos cuerpos, generación de volúmenes de revolución y desarrollos de todo tipo.

4) Podemos generar un gráfico tridimensional, correspondiente a

la

hélice en coordenadas paramétricas: x=Sen(t), y=Cos(t), z=t. Ver figura 1.4. » t=0:pi/50:1O*pi;

» plot3(si n(t),cos(t),t)

)n

b

Figura 1.4

16


Podemos representar una curva plana dada por sus coordenadre r= Cos(2t).Sen(2t) paratvariandoentre 0yn, tomando puntc equidistantes una centésima en el intervalo considerado. ver figura 1.5.

polares

» t=0:.1:2*pi;

r=sin(2*t).*cos(2*t); » polar(t,r) »»

Figura 1.5

6) Podemos realizar un gráfico de una función considerada como simbólica, utilizando el comando "ezplot". Ver figura 1-6 »»

»r

},='¡,13¡1tn2-11'; ezplot(y,[-5,5])

x"3t{x^2- ) 1

1ú

0

-1rt

'5

I

Figura 1-6

En conceptos.

el correspondiente

capítulo

de gráficos ampliaremos

estos

cnpiruro AS

1. rNrRoDuccróru pnÁclrca A MATLAB 17

}IOTACION GENERAL

OS

Siempre que se utiliza un programa, es necesario conocer las características generales sobre notación, con la finalidad de introducirnos en la práctica con el mismo. La mejor forma de aprender Matlab es utilizar el programa. Cada ejemplo consiste en el input del usuario encabezado con el prompt "»" y la respuesta de Matlab en la línea siguiente. Ver figura 1-7. úmBa¡ds to qet started: intro, omnands fop nore infornátion:

d€mo, help hetp help, uhatsneu, info,

subscribe

A:[123;456;7ACl

Figura 1-7

no

En otras ocasiones, y dependiendo del tipo de entrada (input de usuario) que se le proponga a Matlab, devuelve la respuesta utilizando la expresión "ans= n'. Verfigura 1-8

:- 2*2 lt

Figura 1-8 tos

Es necesario poner atención en la diferencia entre mayúsculas y minúsculas, eltipo de paréntesis o corchetes, la cantidad de espacios y en Ia puntuación (comas, puntos y comas).

18


AYUDAS CON COMAI\DOS Ya hemos visto en el capítulo anterior cómo se obtenía ayuda utilizando los menús desplegables de Matlab.

Pero, además, la ayuda también puede obtenerse a través de comandos (instrucciones o funciones), implementados como objetos de Matlab.

Se puede utilizar la ayuda de acceso inmediato para acceder a diversa infor-mación utilizando el comando help. » help HELP topics:

matlab\qeneraJ matTab\ops

natlab\lanq matlab\elmat matfab\effun matlab\ specfun matlab\natfun matfab\datafun matfab\pofyfun matlab\ funfun matTab\sparfun matTab\graph2d matlab\graph3d natlab\ specqraph matTab\qraphics

natfab\uitools matfab\strfun

matfab\iofun matlab\tinefun mat fab\datatypes matfab\winfun (DDE/ActiveX) matlab\demos tooJbox\ syrnl:o 11c

toolbox\tour toofbox\locaf

GeneraL purpose commands. Operators and speciaf characters, Proqramminq Tanquage constructs. Elementary matrices and matrlx manipulation. Elementary math functions. Specialized nath functions. Matrix functions - numericaf finear algebra. Data anafysis and Fourier transforms. Interpofation and poTynomials. Function functions and ODE sofvers. Sparse matrices. Two dimensi onaJ graphs. Three dimensional graphs. Specialized qraphs. Handle Graphics. GraphicaT user interface toofs. Character strings.

File input/output.

Time and dates. Data types and structures. - ttlindows Operating System Interface Fr,ie-.

trxamples and demonstrations. SynboLic Math Toofbox. MATLAB Tour Preferences.

For more help on directory/topic,

type ,,help topic,,.

como vemos, el comando help muestra una lista de los directori:_-

del programa y de su contenido. También se puede pedir ayuda sob.= cualqu¡er tema del contenido utilizando el comando help tema.

» help ¡nv

CAPÍTULO 1. INTRODUCCIÓN PRÁCTICA A MATLAB 19

:iiliiruda

l.larrix inverse. -',-.''Y) is the inve,se of the square matrix X. - ,.,::ning messaqe -is prin¿€'¿ if X is badJv scafed ot :-==::i y sinqular. -:== a-¿so SLASH, PINV, COND, CONDEST, NNLS, LSCOV.

ide ide

j,¡=-''^ded methods :-=- p sym/inv.m

r lrclp matlab\elfun

era

!-=:=::tary math functions. !:-::::cmetric. Szne.

Hyperbolic sine. Inverse sine. rnverse hyperbolic sine. Cos

ine

.

Hyperbollc cosine. Inverse cosine. I nverse hypet bol i c cosine. tattYL!1L.

HyperboTic tanqent. Inverse tangent. Four quadranL inverse tanqenL. Inverse hyperbolic tanqent. Secant. Hyperbolic secant. Inverse secant. Inverse hyperbolic secant. Cosecant. Hyperbolic casecant. Inverse casecant. Inverse hyperbolic casecant. Cotangent. HyperboTic cotanqent. Inverse catangent. Inverse hyperboTic cotanqent.

Exponential. Naturaf loqarithm. Common (base 70) loqartthm. and dissect - Base 2 loqarlthn ctonos sobre

Base 2 povrer and scafe floating Square raot. Next higher power of 2.

floating

point number.

point

2A


abs

conj imag realunwrap

isreal cpTxpair

-

Alcsofute vafue. Phase angle. Complex conlugate. Complex imaqinary part. Complex reaT part. Unwrap phase angJe. True for reaf array. Sort numbers into complex conjuqate pairs.

Roundinq and remainder. fix f foor ceif round mod

rem

stqn

-

tawards zero. towards minus infinity. towards plus lnfinlty. towards nearest integer. Modulus (siqned remainder after division) Round Round Round Round

Remainder

after division.

Siqnum.

Existe un comando de ayuda sobre una determinada secuencia de caracteres (lookfor secuencia) que permite encontrar todas aquellas funciones o comandos que se refieren a dicha secuencia o la contienen. Este comando es muy útil, bien cuando no se dispone de ayuda directa para la secuencia especificada, o bien para ver la ayuda de todos los comandos que contienen la secuencia. Por ejemplo, si buscamos ayuda para todos los comandos que contienen la secuencia complex, podemos usar el comando lookfor complex para ver en qué comandos Matlab ofrece información al respecto. » lookfor complex ctranspose.m: Z' Complex conjuqate transpose. CONJ Complex conluqate. CPLXPAIR Sort numbers into complex conjuqate pairs. IMAG Complex imaqinary part. REAL Complex real part. CDF2RDF Complex diagonal form to reaf bfock dlaqonaL form. RSF2CStr Real bfock diagonal form to complex diaqonal form. B5ODE Stiff problem, linear with complex eigenvalues (85 of EHL) . CPLXDEMO Maps of functions of a complex variabfe. CPLXGRID Polar coordinate complex qrid. CPLXMAP Plot a function of a complex variabfe. GRAFCPLX Demonstrates complex function pJots in MATLAB. ctranspose.m: *TRANSPOSE Syml:olic matrix complex conjuq_ transpose. SMOKE Complex matrix wlth a ,,smoke rinq,, pseudospectrum.

__

capÍruro

1. TNTRoDUCcTóru

pnÁclcnA

MATLAB 21

COMANDOS DE ESCAPE Y SALIDA AL ENTORNO DOS Existen tres formas de salir desde la ventana de comandos de lilaüab al entorno del sistema operativo MS-DOS para ejecutar tareas tsrnporales.

El comando lorden_dos introducido en la ventana de comandos, rerrnite ejecutar la orden de DOS especificada en ambiente Matlab. Por {prnplo:

r ldir :- ¡-clumen de fa unidad D no tiene etiqueta l- ::fmero de serie def vc¡fumen es 145C-72P2 l-:=ctorio de D : \MATLAB52\bin 13/A3/98 0:76

lde

rdos s los ¡ndo

ta"S BAT :-: _10 DLL :3E-Z BAT ::W-)OL BAT :': -:?TS BAT ;I;: DLL :Tr-DLL .}E{ BAT

in

m,_:_,:-
¡llas nen.

para

al

=,- -

- _-t-

_1;sE DAT DLL

13/03/98 19/01/98 21/08/97 ) 27¿ L3/03/98 i4.992 19/01/98 1.973 19/01/98 25.ABB 1B/12/97 16.896 1B/12/97 2.274 13/03/98 470 13/03/98 66.560 02/05/97

1.872 219.136

10 archivo (s)

2 directorio

(s)

0:76 ..

14:14 bccopts.bat 22:24 cLbs770.dJf 0:28 cnex.bat 14:14 conptool.bat 74:74 df50opts.bat 16:34 fenq.dll 15:34 fnat.dfl 0:28 fmex.bat 0:27 license.dat B:34 w32ssi. dff 11.348.865 bytes 159.383.552 bytes fibres

El comando lorden_dos & se utiliza para ejecutar la orden del DOS especif¡cada en modo background. La orden se ejecuta abriendo una rentana de ambiente DOS sobre la ventana de trabajo de Matlab tal y como se indica en la figura 1-9. Para volver a ambiente Matlab basta con pulsar Eri el ratón en cualquier zona de la ventana de comandos, en cuyo caso se cierra automáticamente la ventana de ambiente DOS. Se puede volver a la uentana del DOS en cualquier momento para ejecutar cualquier orden del sbtema operativo pulsando el icono etiquetado Simbolo de MS-DOS situado en la parte inferior de la pantalla.

22


r

,},

,t; i.r

ir »

; i.. *; B

L

i:

Figura 1-9

El comando orden_dos I se utiliza para ejecutar la orden del Dos con pantalla de Matlab.

anteriores, no sólo pueden ejecutarse tipo de ficheros ejecutables o tareas todo comandos del DOS, sino también batch.

con los tres comandos

El comando dos orden-dos se utiliza también para ejecutar la orden del DOS en modo automático sobre la ventana de comandos de Matlab.

Para salir definitivamente de Matlab, basta con teclear quit en la ventana de comandos y, a continuación, pulsar Enter.

MATLAB Y LA PROGRAMACIÓN Combinando adecuadamente todos los objetos definidos en Matlab, adecuados a las reglas de trabajo definidas en el programa, se puede construir un código de programación muy útil en la investigaciÓn matemática' Los programas consisten habitualmente en una serie de instrucciones en las que Se calculan valores, se les asigna un nombre y se reutilizan en cálculos posteriores.

cnpírulo

1. rNTRoDUcclóru

pnÁcrca

A MATLAB 23

Al igual que en lenguajes de programación como C o Fortran, en se pueden escribir programas con bucles, control de flujo e condicionales. En Matlab se pueden escribir programas tales, es decir, definir una secuencia de pasos estándar a Como en C o en Pascal, se puede realizar un cálculo repetitivo Do, For o While. El lenguaje de Matlab también incluye iones condicionales como lf Then Else. Matlab también soporta funciones lógicas, como And, Or, Not y Xor.

Matlab soporta la programación procedimental (con procesos recursivos, bucles...), la programación funcional y la ión orientada al objeto. Veamos dos ejemplos sencillos de ras. El primero genera lamafriz de Hilbert de orden n ,y el segundo Ios números de Fibonacci.

ndo la matriz de Hilbert de orden n 'ttfr+j-f )'; = l:n

ftri ='l¡¡

EI DOS

a(ii)

cnd

= eval(t);

cutiarse r

tareas

a orden ab.

It en

Ia

Madab,

puede máüca. ¡ en las

álculm

$ll;i=l;

los números de Fibonacci

(i)+f(¡-r)<1000 {i+2)=f(i)+f(i+l);

i=i+l

Capítulo 2

ECUACIONES DIFERENCIALES DE R ORDEN. ECUACIONES EXACTAS, EN VARIABLES SEPARADAS, HOMOGENEAS Y LINEALES MUACIONES DIFERENCTALES DE PRIMER ORDEN número de comandos que implementa Matlab relativos a este tema seguir para ya conocidos métodos de resolución los algebráicos Fo{¡lzrma de ecuación diferencial. También se implementan métodos de aproximados de ecuaciones y sistemas de ecuaciones diferenciales. Etr

rrry elevado, pero sí muy eficiente. De todas formas, es posible

fro E

comando básico para resolución de ecuaciones diferenciales es

Este comando computa soluciones simbólicas de ecuaciones ordinarias y sistemas. Las ecuaciones son especificadas por

simbólicas conteniendo la letra D para denotar diferenciación, o D2, D3,...,etc, para denotar diferenciación de orden 2,3,...,etc. A l¡n de la letra D se sitúa la variable dependiente (que suele ser y), letra no precedida por D es un candidato a variable . Si no se especifica la variable independiente, por defecto es se especifica como variable dependiente, la variable independiente es r es la variable independiente por defecto, y en segundo lugar t.

26


Se pueden especificar condiciones iniciales en ecuaciones adicionales, mediante la forma y(a)=b o Dy(a)=b,...,etc. Si las condiciones iniciales no se especifican, las soluciones de las ecuaciones diferenciales contienen constantes de integración C1 , C2, ..., etc. Los comandos más importantes de Matlab que resuelven ecuaciones diferenciales son los siguientes:

dsolve('ecuación','v') Resuelve la ecuación diferencial siendo v la variable independiente (si no se especifica 'v', la variable independiente es x). Solo devuelve soluciones explícitas dsolve('ecuación','condición_inicial',...,'v') Resuelve la ecuación diferencial sujeta a !a condición inicial especificada

dsolve('ecuación','cond1','cond2',...,'condn','v') Resuelve Ia ecuación diferencial sujeta a las condiciones iniciales especificadas

dsolve('ecuación', 'cond1, cond2,..., condn' , 'v') Resuelve la ecuación diferencial sujeta a las condiciones iniciales especificadas

dsolve('ec1' r'ec?' r..,r'ecn','condl','cond2'r...,'condn','v') Resuelve el sistema diferencial sujeto a las condiciones iniciales especificadas (explícitas)

dsolve('ec1, ec2,..., ecn', 'cond1, cond2,..., condn' , 'v') Resuelve el sistema diferencial sujeto a las condiciones iniciales especificadas

maple('dsolve(ecuación, func(var))') Resuelve la ecuación diferencial, considerando var como variable independiente y func como variable dependiente (devuelve soluciones implícitas)

maple('dsolve({ecuación, cond1,cond2, ....condn}, func(var))')

CAP|TULO 2. ECUACIONES DIFERENCIALES DE PRIMER ORDEN 27

-

Resuelve la ecuación diferencial sujeta a las condiciones iniciales especificadas

leS

]es iles nás

maple('dsolve({ec1, ec2,..,ecn}, {func1(var), func2(var),...,funcn(var)})')

los

Resuelve el sistema de ecuaciones diferenciales especificado (devuelve soluciones implícitas) maple('dsolve(ecuación, func(var),'explicit')') Resuelve la ecuación diferencial, ofreciendo la solución en forma explícita, si es posible

t, l?

A continuación se mostrarán algunos ejemplos. Resolvemos en primer lugar ecuaciones diferenciales de primer oreden y primer grado, sin y con valores iniciales. » pretty(dsolve('Dy = a*y')) exp(a x)

» pretty(dsolve('Df =

Cl

f + sin(t)'))

- l/2 cos(t) - l/2 sin(t) + exp(t) Cl

Las dos ecuaciones anteriores también pueden resolverse de s§uiente forma:

r pretty(maple('dsolve(d iff(y(x),x)= a*y,y(x))')) y(x)

s)

:

exp(a x)

_Cl

» pretty(maple('dsolve(diff(f(t),t)=f+si n(t),f(t))'))

f(t):

- 1/2 cos(t) - I/2 sin(t) + exp(t) _Cl

r pretty(dsolve('Dy

= a*Y', 'y(0) = b'))

exp(a x) b

io var riable

r pretty(dsolve('Df

=

f + sin(t)','f(pi/2) = g'¡¡ exp(t)

- l/2 cos(t) - 1/2 sin(f) + 1/2

**-----

exp(l/2 pi)

Ahora resolvemos una ecuación de segundo grado y primer orden.

la

28

ECUACIONES DIFERENfCIALES CON MATLAB

»»

= 6lss¡re('(Dy)^2 + yn2 = 1', 'y(0) = 0', 's')

1¡

I sin(s)] l-sin(s)l

que también puede resolverse de la siguiente forma: » pretty(maple('dsolve({diff(y(s),s)^2 + y(s)"2 = y(s)

-

sin(s), y(s)

:

l, y(0) = 0}, y(s))'))

- sin(s)

Ahora resolvemos una ecuación de segundo orden y pr¡mer grado' » pretty(dsolve('D2y ='an2ny', 'y(0) = 1, Dy(pi/a) = 0')) cos(a x)

A continuación resolvemos sistemas sin y con valores iniciales. »»

»»

pretty(dsolve('Dx = Y','DY ='x')) x(t) : Cl sin(t) + C2 cos(t), y(t) - Cl cos(t)

- C2 sin(t)

pretty(dsolve('Df = 3*f+4*9', 'Dg = '4*f+3*g')) S6) : - C2 exp(3 x) sin(4 x) + Cl exp(3 x) cos(4 x), f(x)

: Cl exp(3 x) sin(4 x) + C2 exp(j x) cos(4 x)

» pretty(dsolve('Df = 3*f+4*9, Dg = {*f+$*g', 'f(0)=0, g(0)=1')) S(x)

:

exp(i x) cos(4 x),f(x)

-

exp(3 x) sin(4 x)

Este último sistema puede resolverse también de la siguienE forma:

»

pretty(maple('dsolve({diff(f(x),x)= 3*f(x)+4*g(x), diff(g(x),x)=4*f(x)+3*9fiL

f(0)=0, g(Q)=l ), ff(x), g(x)))')) {f(x)

:

exp(3 x) sin(4 x), g(x)

2.2 E,CUACIONES

-

exp(3 x) cos(4 x)}

El\ VARIABLES SBPARADAS

Una ecuación diferencial en variables separadas presenta la fonna f(x)dx = g(y)dy. La resolución de este tipo de ecuaciones es inmediata pon¡endo it(x)Ox = J g1y¡dy +C.

CAPITULO 2. ECUACIONES DIFERENCIALES DE PRIMER ORDEN 29

Si Matlab no resuelve directamente la ecuación diferencial con la ffimción dsolve, entonces se sigue el método algebráico usual, que no p'esenta dificultades especiales para el programa, dada su versatilidad en d cálculo simbólico.

En primer lugar intentamos resolverla directamente. La ecuación plede ponerse de la forma: ado.

Cos j

Ix]

y Ix]

-^j

1 + y[x]

2

r dsolve('Dy=y*cos(x)/(l +y^2)') using ::> ¿to¡r" solution could not befound. @{fuir

'".' Error

Luego la ecuación diferencial no es resoluble con dso/ye.

Aplicaremos

el

método algebráico usual de resolución de las

esuaciones diferenciales en variables separadas (separando variables).

r pretty(solve('i nt(cos(x),x)=¡ nt(( 1 +y ^2) ly,y)'l) asin(Joq¡y¡ +1/2v)

La función solución

¡uiente

será x=asin(log(y)+112 y2), o lo que es lo

músmo:

¡99(x)"

sir(x)=log(y)+112 y2 + C

Ahora hallamos el valor C para x=0

e

y=1.

r C=simple('solve(subs(x=0,y=l,sin(x)=l6g ly)+l l2*y^2+G),C)') ra-

fonna nediata a

--t

l La función solución

será

sin(x)=log(y)+1t2y' - 112

30


De la misma forma se puede resolver cualquier otra

ecuaciÓn

diferencial en variables separadas.

Esta ecuación diferencial también es resoluble directamenE utilizando: »»

pretty(maple('dsolve(diff(y(x),x)=y(x)ncos(x)/(1 loq (y (x)

+y(x)"2),y(x))'))

) + 1/2 y (r)2 - sin (x) :

- C7

2.3 BCUACIONES DIFERENCIALES HOMOCÉXNTS Consideremos una ecuación diferencial general de primer grado y primer orden puesta en la forma M(x,y)dx=N(x,y)dy. Esta ecuaciÓn se die que es homogénea de grado n, si lo son a lavez lasfunciones M y N, es decir, si se cumple simultáneamente: M(tx,tY)=t

n

M(x,Y)

N(tx,tY)=t

n

N(x,Y)

Para este tipo de ecuaciones, el cambio de la variable x por la variable v tal que x=vy, transforma Ia ecuación diferencial inicial (de variables x e y) en otra de variables separadas (de variables u e y). Se resuelve la ecuación en variables separadas y posteriormente se deshace el cambio.

Vamos a resolver la ecuación por el método algebráico ordinario para ecuaciones diferenciales homogéneas. En primer lugar comprobamos si la ecuación es homogénea » maple('m ;=(x,y)-)x^ 2-y n2'l; >» maple('n ¡ =(x,y)->x*y') ; »

facto('m(t*x,t*y)')

t^2*(x-y)*(x+y)

» factor('n(t*x,t*y)')

L

CAPíTULO 2. ECUACIONES DIFERENCIALES DE PRIMER ORDEN 31

on t"2*x*y 'rte

Luego la ecuación es homogénea de grado 2. Para resolverla aplicamos el cambio de variable x=vy.

Antes de realizar las operaciones del cambio de variable, es conveniente cargar la librería difforms, mediante el comando maple('with(difforms)'), que nos permitirá trabajar con formas diferenciales. Una vez cargada esta librería también es conveniente utilizar el comando maple('defform(v=0,x=0,y=Q)'¡, que permite declarar todas

las variables que no van a ser constantes o parámetros en

oy

diferenciación.

lice

» maple('with(difforms)');

CS

»r

maple('defform(v=0,x=0,y=Q)')'

Ahora realizamos ya el cambio de variable x=vy, términos en e¡

rla (de Se

are

d(v) y

ario

y

agrupamos

d(y).

¡fy('su bs(x=v*y,m(x,y)*d (x)+n(x,y)*d(y))')

:rrsrm v^2 *y^ 3 * d(v) +v^3 *y" 2 * d (y) -y"3 * d (v)

» pretty(maple('collect(v^2*y^3*d(v)+v^3*y^2*d(y)-y^3*d(v),{d(v),d(V)})'))

233 (v y -y)d(v)+v

ffi

la

32

y d(y)

Si dividimos la expresión anterior por v3 y3 , y agrupamos términos en d(v) y d(y), tenemos ya una ecuación en variables separadas.

mos

» pretty(maple('collect(((v^2*y^3-y^3)*d(v)+v^3*y^2*d(y))/(v^3*y^3),

{d(v),d(y)})'))

23 (v y -y)d(v)

33v

d(y)

--- + ----

vy Ahora obtenemos simplificados totalmente.

la expresión anterior con todos sus

términos

32


»

pretty(maple('convert(collect(((v^2*y^3-y^3)*d(v)+v^3*y^2*d(y)/(v^3*y^3), {d(v),d(y)}), parfrac,y)')) 2

!:

_ _ _ _1_'_

_1!Y'_

.

1!u_',_

,tY Ahora resolvemos la ecuación en variables separadas

pretty(simple('i nt((v^ 2-1llv ^3,v1+¡

nt(1

/y,y)')) 1

loq(v) +----+loq(y) 2u2 Ahora deshacemos el cambio de variable » pretty(simple('subs(v=ly,log (v)+1 l2lv "2+log(y))')) 2 -Los

(x) *

tlz -!-*'

Ya estamos en condiciones de escribir la solución general de la ecuación diferencial inicial. 2

log(x) + 112

v

----- =

C

2 X

Ahora se representan gráficamente las soluciones de esta ecuación dlferencial. Para ello nos basamos en que el gráfico de la familia de soluciones con parámetro c, es equivalente al gráfico de las curvas de nivel correspondiente a la ecuación sin constante (figura 18-1) [x,y]=mes hgrid(0. I =0.05:1 2=y.^2 J (2*x. ^2)+l og (¡) ; »» contour(2,65) »»

¡y

12,-1 : 0.05 : I

);

cnpírulo

2. EcuActoNEs DtFERENctALEs DE pRtMER oRDEN 33

y^3), 4rl J5

30

lq 2U

15 1ú

5

Figura 2-1

2.4 ECTJ ACIONES DIFEREI\CIALES EXACTAS La ecuación diferencial M(x,y)dx+N(x,y)dy=0 es exacta si se cumple óN/óx = 1Mlóy.Si una ecuación diferencial es exacta, existe una función F tal que su diferencial total dF coincide con la ecuación, o sea:

que

dF=M(x,y)dx+N(x,y)dy

de

la

por lo tanto F(x,y)=c será la familia de soluciones de la

ecuación

diferencial.

A continuación se muestra un ejercicio que sigue el

método

algebráico normal de resolución de este tipo de ecuaciones.

uaclon ilia de vas de

En primer lugar intentamos resolver la ecuación con dso/ye » maple('m ; =(x,y)->'l +y*exp(x*y)+y*cos(x*y)') » maple('n : =(x,y)->1 +x*exp(x*y)+xncos(x*y)') ; » dsolve('m(x,y)+n(x,y)*Dy=Q') ??? Error using::> ¿to¡rn Explicit solution could not

be

;

found.

Se observa que la función dsolve no resuelve Ia ecuación propuesta. Vamos a intentarlo por el método algebráico clásico de resolución de ecuaciones diferenciales exactas.

34


En primer lugar comprobamos si la ecuación diferencial propuesta es realmente exacta » pretty(si mple(d iff('m(x,y)','y'))) exp(y

x) + x y exp(y x) + cos(y x) - x sin(y x) y

» pretty(simple(diff('n(x,y)','x'))) exp(y

x) + x y exp(y x) + cos(y x) - x sin(y x) y

Por ser exacta, su solución se puede hallar de la siguiente forma: » sol ucion I =simplify('i nt(m(x,y),x)+g(y)') solucionl

:

-x + exp (y *x) 1 s in (y*x) + g(y)

Ahora calculamos la expresión de la función g(y), bajo la siguiente condición : diff(int(m(x,y),x)+g(y),y)=n(x,y) » pretty(simpl ify('i nt(m(x,y),x)+g(y)')) -

x + exp(y x) + sin(y x) + g(y)

» pretty(simplify('diff(-x+exp(y*x)+sin(y*x)+g(y),y)')) d

x exp(y x) + x cos(y x)

»»

* -;rfrr,

simplify('solve(x*exp(y*x)+x*cos(y*x)+d¡ff(g(y),y)=n(x,y),diff(g(y),y))')

ans

:

I

Luego 9'(y)=1, con lo que la solución final será, salvo constante: » pretty(simplify('subs(g(y)=int(l,y),-x+exp(y*x)+sin(y.x)+g(y))')) -

x + exp(y x) + sin(y x) + y

Para representar gráficamente la familia de soluciones, graficamos las curvas de nivel de la solución sin constante (figura 18-2) » [x,y]=meshgrid('2*pil3:.2:2* pilSl:' »r z=-X+€XP(y.*x)+s¡ n(y.*x)+y; » contour(2,100)

CAPíTULO 2. ECUACIONES DIFERENCIALES DE PRIMER ORDEN 35

Ha 18 .18 ,14

lr 1D

I

Figura 2-2 De la misma forma se puede resolver cualquier ecuación diferencial reducible a exacta mediante un factor integrante.

fite

2.5 ECUACTONES DIFERENCIALES

LINEALES

Una ecuación diferencial lineal de primer orden es deltipo:

dy/dx + P(x)y = Q(x) con P(x) y a(x) polinomios.

Las ecuaciones diferenciales de este tipo se transforman

en

exactas mediante el factor integrante: J

P(x)ox

e

y su solución general viene dada por la expresión:

hr

J

p(x)ox

(e I

)(le

I e1x¡ox Q(x)dx)

Matlab implementa estas soluciones de las

ecuaciones y y lineales, las diferenciales ofrece siempre cuando la integral del factor para integrante sea evaluable el programa.

»»

pretty(s

i

m

ple(dsolve('x*Dy+t*y=x*s i n (x)'))) sin (x) cos (x) C1 - 6 sin (x) - cos(x) + 3 ------ + 6 -----* + -------x23

36


Otra forma más eficiente de resolución de la ecuación diferencial, que incluye soluciones implícitas, es la siguiente: » pretty(simplify('dsolve(x*diff(y(x),x)+3*y1x)=x*si n(x),y(x))'))

32' x cos(x) - 3 x sin(x) + 6 sin(x) - 6 x cos(x) Y l"/

-

---------;

_C7

al,

Capítulo

3

ECUACIONES DIFERENCIALES DE ORDEN ST]PEROR TRANSFORMADADE LAPLACE Y TIPOS ESPECIALES DE ECUACIOI{ES

3.1 ECUACIONES ORDINARIAS DE ORDEN SUPERIOR Las ecuaciones diferenciales lineales ordinarias de orden

n

tienen la siguiente forma general:

(k) n I ak(x) y (x)=a0(x)y(x)+a1

(n)

(x)y'(x)+¿l(x)y"(x)+....+ an(x)y (x)=f(x)

k=0

Si la función f(x) es idénticamente nula la ecuación se llama y si f(x) no es Ia función cero, la ecuación se llama no

homogénea

homogénea. Si las funciones ai(x) (i=1...n) son constantes, la ecuación se denomina en coeficientes constantes.

Un concepto de gran interés aquí va a ser el de conjunto de funciones linealmente independientes. if1(x), f2(x),....,fn(x)) es un conjunto de funciones linealmente independientes si para algún x del intervalo común de definición, el determinante wronskiano de las funciones es no nulo. El determinante wronskiano de las funciones, en cualquier punto x de su campo de definición común, se define de la siguiente forma:

38


f1(x) f1'(x) f1"(x)

f2(x) f3(x).......,............,.fn(x) (x) f3'(x).....................fn'(x)

f2'

f2"

(x) f3"(x)....................fn"(x) = W(x)

(n-1) (n-1) (n-1)

(n-1) (x)

f1 (x) f2 (x) f3 (x).................fn

Matlab dispone del comando maple('Wronskian') que permite calcular la matriz wronskiana de un conjunto de funciones. Su sintaxis es:

maple('Wronskian(V,variable)')

Calcula la matriz wronskiana correspondiente al vector de funciones V de variable independiente x

Un conjunto S = {f1(x),.....,fn(x)} linealmente independiente de

soluciones no triviales de la ecuación lineal homogénea de orden n: (n) a0(x)y(x)+¿1 (x)y'(x)+¿2(x)y"(x)+....+

an(x)y (x) = 0

se llama conjunto fundamental de soluciones de la ecuación.

Si las funciones ai(x) (i=1....n) son continuas en un intervalo abierto fundamental de soluciones § = {fi(x)} en i.

l, entonces la ecuación homogénea tiene un conjunto

Además, la solución general de la ecuación homogénea vendÉ dada por la función: n

f(x)

=

¡

ci

fi(x)

donde {ci} es un conjunto de constantes arbitrarias

i=0

La ecuación:

2nni

a0+a1 m+a2m + ...... +anm

=I

ai

m=0

i=0

denomina ecuación característica de la ecuación diferencid homogénea en coeficientes constantes. Las soluciones de esta ecuación característica van a determinar las soluciones de la ecuación diferencid lineal homogénea general en coeficientes constantes.

se

t-

CAPíTULO 3. ECUACIONES DIFERENCIALES DE ORDEN SUPERIOR 39

»»

funciones=maple('vector([exp(x),xnexp(x), x^2*exp(x)])')

funciones

lite

;tI

i

tal Ex ,de I

i

* I up (*), x*exp (x), x^2 exp (x)J

>r

W=maple('Wronskian(funciones,x)')

W:

k*p6),

x"2*exp(x)l

x*exp(x),

[exp(x), exp(x)+x*exp(x), 2*x*exp(x)+x"2*exp(x)J kxp (x),

2 *exp (x)

+x*exp(x),

2

*exp(x)

+

4*x*exp (x) +x"2 *"*p (*)l

» pretty(determ(W)) 3

2 exp(x)

Ya tenemos el valor del wronskiano, que evidentemente es siempre distinto de 0, luego el conjunto de funciones es linealmente independiente. Prto

ld"

3.2 ECUACIONES LINEALES DE ORDEN SUPERIOR

HOMOGENEAS EN COEFTCIENTES CONSTA¡ITES.

I

ldrá

La ecuación diferencial lineal homogénea de orden n: (n) aO(x)y(x)+a1 (x)y'(x)+¿2(x)y"(x)+....+ an(x)y (x) = 0

se dice que es de coeficientes constantes cuando las funciones ai(x) (i=1,..,n) son todas constantes (no dependen da la variable x). La ecuación: 2

a0+a1 m+a2m + ...... +anm

=I

nt ai m=0

i=0 hcial ErcN i

t*rl

denomina ecuación característica de la ecuación diferencial homogénea en coeficientes constantes. Las soluciones de esta ecuación característica (m1 , fr2, ...mn) van a determinar las soluciones de la ecuación diferencial lineal homogénea general en coeficientes constantes.

se

40


Si los mi (i=1..n) son todos distintos, una solución general de

la

ecuación homogénea en coeficientes constantes es: m1

y(x) =c1

x

m2x mnx ...... +cn e

e + c2e +

c1, c2,...., cfl son constantes arbitrarias que determinan la familia de soluciones de la ecuación diferencial.

Si algún mi es una raíz de multiplicidad k de la

ecuación

característica, determina los siguientes k términos de la solución:

mix

mix

cie +c(i+1)xe

+c(i+2)x

2 mix

e

+........ +c(i+k)x e

k mix

Si existe alguna raíz compleja mj= a+bi de la ecuación característica, entonces también existe la raíz compleja conjugada a-bi. Estas dos raíces determinan simultáneamente el siguiente par de términos de la solución generalde la ecuación homogénea: cj

e

Cos(bx)+c(i+1)e Sen(bx)

Matlab aplica directamente este método de obtención de las soluciones de las ecuaciones lineales homogéneas en coeficientes constantes, al utilizar el comando dsolve o maple('dsolve').

» pretty(dsolve('3*D2y+2*Dy-5*y=0')) Cl

exp(x) + C2 exp( 5/j x)

» pretty(dsolve('2*D2y+5*Dy+5*y=0','y(0)=0, Dy(Q)=l

t/2 2/15

1/2

3 5

1/2

/!'¡¡

1/2

exp( 5/4 x) sin(l/4

3 5

x)

» pretty(simple(dsolve('9*D4y-6*D3y+46*D2y-6*Dy+37*y=g'¡¡¡

cnpírulo

I l¡)

3. EcuActoNES DtFERENctALES DE oRDEN supERtoR 41

sin(2x) + C2 exp(l/i x) cos(2 x) + C3 sin(x) + C4 cos(x)

solución, es evidente que la ecuación pares de soluciones complejas conjugadas. tiene dos

Si nos fijamos en la

rmilia

^4-6¡kx^ 3+46*x^2 -6*x+37=0')

&

ecr¡acilim

También puede resolverse la ecuación de la forma siguiente:

ple('dso Ive(9.d iff(y(x),x$4) -6*d iff(y(x),x$3)+46.d iff(y(x),x$2),x)+37*y(x)=0,y(x))'))

ecuacIlft gada a{üü" nte Par 3

xt sin(2 x)+ _C2 exp(l/3 x) cos(2 x)+ _C3 sin(x)+ C4 cos(x)

ACIONES 1\O HOMOGÉNEAS CON COEFICM,NTES IONSTAI\TES, VARIACION DE PARAMETROS Dada la ecuación lineal no homogénea en coeficientes constantes: (k)

y

(n) (x)=a0(x)y(x)+at (x)y'(x)+¿l(x)y"(x)+....+ an(x)y (x)=f(x)

ú1(x), y2(x),.....,yn(x)) un conjunto de soluciones

linealmente

de la correspondiente ecuación homogénea: n

)+a1 (x)y'(x)+a2(x)y"(x)+....+ an(x)y (x)= 0

Una solución particular de la ecuación no homogénea viene pot. n

I

ui(x) yi(x) donde las funciones ui(x) se obtienen como sigue:

É1

f(x) Wi (y1 (x), y2(x),.....,yn(x))

W

(y1

(x), y2(x),.....,yn(x))

42


Wi (y1(x), y2(x),....,yn(x)) es el determinante de la matriz que resulta al

sustituir la i-ésima columna del Wronskiano W(y1(x), y2(x),.....,yn(x)), por el transpuesto del vector (0,0,.....,0, 1 ).

La solución general de la ecuación no homogénea viene dada por la solución de la homogénea más la soluciÓn particular de la no homogénea, y tendrá la forma:

mnx x m2x +c2e+......+cne + Yp(x)

m1

y(x)=c1 e

supuestas las raíces mi de la ecuación característica de la homogénea todas distintas. Si no lo fueran, nos remitimos a la forma general de la solución de una ecuación homogénea, que ya hemos estudiado antes.

Vamos a seguir el camino algebráico del método de variaciÓn de parámetros para resolver la primera ecuación. Consideramos la ecuaciÓn característica de la homogénea para obtener un conjunto de solucione linealmente independientes. r»

solve('ma2¡4*¡¡+1 3=0')

ans

:

[-2+i*ü t-2-3*vl

maple('f :=x->x*cos(3*x)^2') ; » maple('y1 :=x->exp(-2*x)*cos(3*x)'); >» maple('y2 : =x-)€xp(-2*x)*si n(3*x)') ; » maple('W¡=¡-)Wrorskian([yl (x),y2(x)],x)') » pretty(si mplify(maple('det(W(x))')))

»»

j

exp(

;

4 x)

Ya tenemos el wronskiano no nulo, lo que indica que las son linealmente independientes. Ahora calculamos las funciones i=1,2 >»

maple('Wl :=x->array([[0,y2(x)],[1,diff((y2Xx),x)]l)');

» pretty(simplify(maple('det(W1 (x))')))

\-

cnpÍruro

ta al or el

-

exp(

3. EcUACtoNES DtFERENCIALES DE oRDEN supERtoR 43

2 x) sin(3 x)

» maple('W2:=x>array(tly1 (x),01,[diff((y1 )(x),x), 1 ]l)'); » pretty(si m pl ify(maple('det(W2(x))'))) exp(

Jada

lno

2

x) cos(3 x)

Ahora calculamos

la

solución particular

de la ecuación no

homogénea.

r maple('ul:=x->factor(simplify(int(f(x).det(Wf

(x))/det(W(x)),x)))');

»

maple('ul(x)')

I

14652300*exp(2*x)*(129285*cos(9*x)*x-6084*cos(9*x)-28730*sin(9*x)*xI 87850*sin(3*x)*x-361 25*sin(3*x)) 1 3*sin(9*x)+28 I 775*cos(3*x)*x-86700*cos(3*x)-

inea

bla

I 30

EFI F: ,I W,Li E-:-

E:I

»

maple('u2(x)')

|

&J rde

» maple('u2:=x->factor(simplify(int(f(x).det(W2(x))/det(W(x)),x)))');

l.' I 4 6 5 2 3 00*exp(2 *x) * ( 5 6 3 5 5 O*cos (3 *x) *x+ I

)60100*sin(3*x)+287j0*cos(9*x)*x+13013*cos(9*x)+

ción

>»

nes

» maple('yp(x)')

+ B 4 5 3 2 5 *s in(3 *x) *x129285*sin(9*x)*x-6084*sin(9*x))

08 3 7 5 *cos (3 *x)

maple('yp :=x-)factor(si mpl ify(y1 (x)*u 1 (x)+y2(x)*u2(x)))')

;

*x)"2 + 2 4/ I I 0 5 *cos (3 *x) * 5 *x*cos (3 *x)"2 + 1 3 4 3 6/ I 2 2 I 0 2 5 *cos (3 *cos *x) *sin(3 *x) *x) *x+ + 5 4/ I I 05 *x- 2 I I 68/ I 2 2 I 02 5 (3 sin(3 3 8 5 2/I 2 2 I 0 2 5 -

/

23 110

Luego ya podemos escribir la solución general de la no homogénea »

*y 1 (xl+ c2* y 2{x) +yp (x} m a p I e ('y : =x .> s ¡ m p I ify( c I }' ) ; 0,",'"om b i ne(y(x),tris)')

:,: *exp

*x) *cos (3 *x)

*exp(2 *x) *sin(j *x)-2

3 /2 2 I 0*x*cos (6*x) + 1/26*x+6718/1221025*cos(6*x)-2/169+12/1105*x*sin(6*x)+ 1926/1221025*sin(6*x)

c1

nes

i(x)

(- 2

+ c2

Ahora representamos gráficamente un conjunto de soluciones, para determinados valores de c1 y c2 (figura 18-3) » fplot(simpl ify('su bs(c1 =-5,c2=-4,y(x))'), [-1, I ])

» hold on » fplot(si mpl ify('s ubs(c1 =-5,c2=4,y(x))'), [-1, I ]) » fplot(simplify('subs(c1 =-5,c2=2,y(x))'),[-1, I ])

44


» fplot(simplify('subs(cl =-5,c2='2,y(x))'),['1,1 ]) » fplot(simplify('subs(c1 =-5,c2=-l,y(x))'),['1,1]) » fplot(simplify('subs(cl =-5,c2=1,y(x))'),[-1,1 ]) » fp lot(s i m pl ify('su bs (c1 =5,c2='l,y(x))'), [-1, 1 ]) » fplot(simplify('subs(c1 =5,c2=-1,y(x))'),['1,1]) >» f pl ot(s i m pl ify('s u bs(c I =5, c2=-2,y(x))'), ['1, 1 ]) » fp lot(si m pl ify('s u bs (c1 =5,c2=2,y(x))'), ['1, 1 ]) » fp lot(s i m p ! ify('s u bs (c1=5,c2=4,y(x))'), ['1, 1 ]) » f plot(s i m pl ify('s u bs (c1 =5,c2=4,y(x))'), [-1, 1 ] ) » fplot(simplify('subs(c1 =0,c2='4,y(x))'),['1, I ]) » fp lot(s i m pl ify('s u bs (c1 =0,c2=4,y(x))'), ['1, 1 ]) »» f p I ot(s i m p I ify('s u bs (c1 =0,c2=2, y(x) )' ), [-1, I ] ) » fplot(simplify('su bs(c1 =0,c2='2,y(x))'), ['1, I ]) » fplot(simpl ify('subs(c1 =0,c2=-1,y(x))'), ['1, 1 ]) » fplot(si mplify('subs(c1 =0,c2=1,y(x))'), ['1, I ])

Figura 3-1

Para la segunda ecuación diferencial aplicamos directamente DSolve y tenemos la solución. » pretty(si mple(dsolve('D2y-2*Dy+y=sxr(x).log(x)')))

1/4exp(x) (2 log(x)x

22 *3x

+4C7+4C2x)

3.4 ECUACIONES NO HOMOCÉXNAS CON COEFICM,NTES

VARIABLES. ECUACIOI{ES DB CAUCI{Y.EULER

CAPíTULO 3. ECUACIONES DIFERENCIALES DE ORDEN SUPERIOR 45

La ecuación lineal no homogénea en coeficientes variables:

n k(k) I ak x y(x)=a0y(x)+a1

n (n) y'(x)+a2 y"(x)+....+ x x an x y (x)=f(x)

k=0

se denomina ecuación de Cauchy-Euler.

Matlab resuelve directamente este tipo de ecuaciones con

el

comando dsolve o maple('dsolve').

» pretty(simple(dsolve('x^3*D3y+1 6*x^2*D2y +f $*a*Dy+1 25*y=0'))¡ C1 + C2 sin(3 log(x)) x + C3 cos(3 log(x)) x

,t 3.5 TRANSFORMADA DE

LAPLACE

Sea una función f(t) una función definida en el intervalo [0,co). La transformada de Laplace de f(t) es la función de variable s:

*oo

_St

L{f}(s)=i e

f(t)dt,

0

Decimos que f(t) es la transformada inversa de Laplace de F(s) si se cumple: -1

,mente

L{fXs) = F(s) y escribimos

L {F(s)i

(t) = f(t)

Matlab aporta los comand os maple('laplace') y maple('invlaplace) que permiten calcular la transformada de Laplace de una expresión respecto a una variable y la transformada inversa. su sintaxis es la siguiente:

maple('laplace(expresión, t, s)')

iTES

Calcula la transformada de Laplace de la expresión con respecto a t. La transformada es de variable s maple('invlaplace(expresión,s,t)')

46


Calcula la transformada inversa de Laplce de la expresión con respecto a s. La transformada inversa es de variable t Veamos algunos ejemplos. » pretty(maple('laplace(t^(3/2)-exp(t)+sinh(a*t), t, s)'))

;

1/2

piTa 3/4 ----5/2 s-7

+ -------

2

2

» pretty(maple('invlaplace(s^2/(s^2+an2¡n(3/2), s, t)')) - t BesselJ(l, a t) a + BesselJ(0, a t)

La transformada de Laplace y su inversa se utilizan para resolver determinadas ecuaciones diferenciales. El método a seguir es calcular las transformada de Laplace de cada término de la ecuación para obtener una ecuac¡ón diferencial en transformadas de Laplace y resolverla en func¡ón de las transformadas de Laplace. Finalmente se halla la solución de Ia ecuación inicial calculando la transformada inversa de Laplace de la solución de la ecuación en transformadas de Laplace. No obstante, Matlab aporta la opción 'laplace' en el comando maple('dsolve'), que tuerza al programa a resolver la ecuación diferencial utilizando transformadas de Laplace. La sintaxis es la siguiente:

maple('dsolve(ecuación, func(var),'laplace')

En primer lugar calculamos las funciones transformadas de Laplace miembros de la ecuación diferencial, y aplicamos las condiciones iniciales.

de los dos

maple('L:=s->laplace(d iff(y(x),x$2)+2*d iff(y(x),x)+4*y(x),x,s)') » pretty(simplify('subs(y(O)='¡,(D(y))(0)=r,L(s))'))

>»

2

;

Japlace(y(x), x, s) s - s - 3 + 2 Taplace(y(x), x, s) s + 4 laplace(y(x), x, s)

CAPÍTULO 3, ECUACIONES DIFERENCIALES DE ORDEN SUPERIOR 47

» maple('Ll : =s->laplace(x-exp(-x),x,s)') » pretty(simpl ify('Ll (s)'))

;

s-7+s

"' ," * ,, A continuación integramos la ecuación diferencial en transformadas de Laplace » pretty(simplify(maple('solve(L(s)=¡1 1"¡, laplace(y(x),x,s))')))

y(0)+3s solver hr las er una unciÓn

¡dela

33 y(A)+D(y) (0) s +D(y)(A)s "' ,"*r,

+ 2 y(A)

2

s +s+1

("2 *2s+4)

Ahora sustituímos en la solución las cond¡ciones iniciales dadas. » maple('TL: =s-)solve(L(s)=11 1s),laplace(y(x),x,s))') » pretty(si mpl ify('su bs(y(0)={,(D(y))(0)=1,TL(s))'))

de la +4s

manÓ xencid

u' ,"*r,

32 +2s

;

+s+7

("2 *2s+4)

Ya hemos obtenido la transformada de Laplace de la solución de la ecuación diferencial inicial. Para calcular la solución de la ecuación inicial calcularemos la transforma inversa de Laplace de la función obtenida en el pso anterior.

r[3Plm snos tÉ

n maple('TL0 :=s->simpl ify('subs(y(0)=1,(D(y))(0)=l,TL(s))')') r solucion=simple(maple('i nvlaplace(TL0(s),s,x)')) ; n petty(solucion) 7

It

7/4 x - 7/8 - -------- + 5/B 3 exp (x)

sin

1/2 (3 x) exp (x)

1/2 3

35

;

1/2

cos

24

(3

x)

exp (x)

Ya tenemos la solución de la ecuación diferencial inicial. También se podía haber resuelto directamente mediante:

pretty(s

i

m ple(m

aple('dsolve({d iff(y(x),x$2)+2.d iff(y(x),x)+4*y(¡)=

),y(0)=1, D(y)(0)=1 ),y(x), laplace)')))

48


1/2

y(x) :7/4

3.6 POLINONIIOS

1/2

1/2

1 x) 3 sin(3 x - 1/8 - -------- + 5/8 ----3 exp (x) exp (x)

+

_:: _

"_"_:_

24

!:_

_ _ _

_:)_ _

exp (x)

ORTOGOI\ALES

Dos funciones distintas f(x) y g(x) se d¡ce que son ortogonales en un intervalo [a,b] si su producto interno es 0, es decir si ln

J f1x)g(x):0 ' ,n

ejemplo de familia ortogonal pueden ser las funciones

fn(x)=sen(nx)

y

g"(x)=cos(nx), n=1,2,3,... en el intervalo [-rc,r].

Matlab aporta una amplia lista de polinomios ortogonales, que van a ser muy útiles a la hora de resolver ciertas ecuaciones diferenciales no lineales de orden superior. Las funciones simbólicas que permiten el trabajo con estos polinomios son las siguientes:

T(n,x)

Polinomios de Chebychev de primera especie.

U(n,x)

Polinomios de Ghebychev de segunda especie.

P(n,x)

Polinomios de Legendre

H(n,x)

Polinomios de Hermite.

L(n,x)

Polinomios de Laguerre.

L(n,a,x)

Polinomios generalizados de Laguerre.

P(n,a,b,x)

Polinomios de Jacobi.

G(n,m,x)

Polinomios de Gegenbauer.

Veamos ahora su relación con las ecuaciones diferencialesPrecisamente esta relación es la que permite hallar soluciones de determinadas ecuaciones no lineales de orden superior. Para usar estre funciones es necesario ejecutar antes maple('with orthopoly').

3.6.1 Polinomios de Chebychev de primera y segunda especic Los polinomios de Chebychev de primera especie son las soluciones de la ecuación diferencial:

l

L

CAPÍIULO 3. ECUACIONES DIFERENCIALES DE ORDEN SUPERIOR 49

(1-

^')

y'-x'y-t

r'y :

0 n: 0,1,2,...

Su ortogonalidad viene dada por:

T.(x)T,(x)

ú-í'

dx:0 m+n

Los polinomios de Chebychev relación de ortogonalidad :

de segunda especie cumplen

la

1

?

J

unlx¡um(x)(1-x2)1/2

dx=g

m+n

-1

3.6.2 Polinomios de Legendre a

Son soluciones de la ecuación diferencial de Legendre:

D

o (1 -

^')

y'-2x'yt

n(n+ 1) y

:

0 1

Su ortogonalidad viene dada por la relación

I pn(r)pr(x)dx=0

m*n

3.6.3 Polinomios asociados de Legendre Las soluciones de la ecuación diferencial:

(l-*') y'-2x'y+{n(n+f t

9tv:O l-x-

son las llamadas funciones asociadas de Legendre

3.6.4 Polinomios de Hermite les. de stas

Las soluciones de la ecuación diferencial de Hermite:

y'-2x'y*2ny:

O¡Error! Argumento

de modilicador no especificado.

se conocen como polinomios de Hermite. BCre

ones

., _* Su ortogonalidad viene dada por, I Hn(x)Hm(x)e dx=0 m+n -1

¡Error! Argumento de modilicador no especificado.

50


3.6.5 Polinomios generalizados de Laguerre Las soluciones de la ecuación diferencial general de Laguerre:

xy"+(a+1-x)Y'+oY=0 se conocen como polinomios generalizados de Laguerre. 1


I

J

Ln(x)Lm(x)

x" e-'dx=o

m+n

-1

3.6.6 Polinomios de Laguerre Las soluciones de la ecuación diferencial de Laguerre:

x'y+(1-x)y'iny:0 Se conocen como polinomios de Laguerre. Se trata del caso particular a=0 de los polinomios generalizados de Laguerre

3.6.7 Polinomios de Jacobi 1


I

pn(r)prn(x)(1-x)" (1+x)b dx=O m+n

-1

3.6.8 Polinomios de Gegenbauer Su estructura es la siguiente: I

f(a+

G(n,a,x): (-2 )"

1

;)r(n+2a)

n!r(2 a)r(a+ r+

]Xr

-

*')^)

t9l" dx

(1 -

¡2)"*"-i

» pretty(simple(maple('T(7,x)'))) 7

64x I

I

I

t_

-772x

+56x

3

-7x

CAPiTULO 3. ECUACIONES DIFERENCIALES DE ORDEN SUPERIOR 51

» pretty(s¡mple(maple('P(6,x)')))

23t 6 ---x-*--x t6

315

105 +---x

16

2

-5/16

16

» pretty(simple(maple('H(5,x)')))

53 32x -16Ax +720x » pretty(si m ple(maple('L(5,x)')))

1-5x+5x 3.7 FIINCIONIES DE

{

234 -5/3x

+5/24x -7/720x

AIRY Y BESSEL

Reciben el nombre de funciones de Airy las soluciones linealmente independientes de la ecuación diferencial lineal de segundo orden: Y" - xY =

0

(ecuación deAirY)

Reciben

el

nombre

de funciones de Bessel las

soluciones

linealmente independientes de la ecuación diferencial:

y" + y'lx

* (k'

- n2 lx2 ¡y = g (ecuación de Bessel)

Reciben

el nombre de funciones de Bessel modificadas las

soluciones linealmente independientes de la ecuación diferencial:

y" + y'lx- (k' *

n2

/x) y = 0 (ecuaciÓn de Bessel modificada)

Matlab implementa al respecto las siguientes funciones simbólicas:

Ai(z) y Bi(z)

Dan soluciones independientes de Ia ecuación diferencial de Airy.

BesselJ(n,z) y BesselY(n,z) Dan soluciones independientes de !a ecuación diferencial de Bessel.

Bessell(n,z) y BesselK(n,z) Dan soluciones independientes de la ecuación diferencial de Bessel modificada.

52

ECUACIONES DIFERENGIALES CON MATLAB

Se trata de la ecuación diferencial de Bessel para n=112. Obtendremos dos soluciones linealmente independientes de la siguiente forma: » pretty(simple(maple('BesselJ

(1

/2,x)')))

2

sin (x)

1/2 pix

1/2

» pretty(simple(maple('BesselY(1 /2,x)'))) 1/2

2

cos (x)

1/2

1/2

MGS

Capítulo 4

ECUACIONES DIFERENCIALES POR

METODOSAPROKMADOS 4.1 BCUACIONES DE ORDEN Y GRADO SUPERIOR A UNO,

LINEALES YNO LINEALES, NMTONOS APRO)ilMADOS Cuando los métodos conocidos de resolución de de ecuaciones diferenciales y sistemas no ofrecen solución, se suele recurrir a métodos de aproximación. Dentro de los métodos aproximados los hay simbólicos y numéricos. Los simbólicos permiten obtener soluciones algebráicas aproximadas, y el más representativo es el de la series de Taylor. Los numéricos permiten obtener conjuntos de puntos de las soluciones, que posteriormente pueden ajustarse por cualquier método algebráico (interpolación, regresión,....) a la curva que mejor les aproxime. Dicha curva será una solución aproximada de la ecuación diferencial. Entre de los métodos numéricos más comunes se halla el método de Runge Kuta. El uso más corriente de los métodos aproximativos es la resolución de ecuaciones y sistemas de ecuaciones diferenciales no lineales de orden y grado superiores a uno, siempre que la solución exacta no sea posible de obtener mediante otros métodos.

4.2 EL VTÉTOOO DE LAS SERIES DB TAYLOR

Este método ofrece soluciones polinómicas aproximadas de ecuaciones diferenciales generales, y se basa en el desarrolo de funciones

54


en series de potencias de Taylor. Matlab ofrece la opción 'series'

del

comando maple('dsdve'), que permite obtener este tipo de soluciones. Su sintaxis es la siguiente:

maple('dsolve(ecuación, func(var),'series')

Además existe el comando maple('powsolve') que resuelve ecuaciones diferenciales lineales por series de potencias, y cuya sintaxis es la siguiente: maple('powseries[powsolve](ecuación, condl,....,condn)') Con el comando maple('conveñ(polynom)') se puede convertir una solución complicada a polinomio ordenado en potencias de la variable.

pretty(s i mp le(m apl e('dsolve(4*x^2*d iff(y(x),x$2)+4*x*d iff(y(x),x)+ (x^2-1 )*Y(x)=0,Y(x),series)')))

y(x) :

2 46 t/2 (1 - 1/24 x + 7/1920 x + O(x )) Cl x 6

1/2

1/2

»» pretty(simplify(maple('convert(_C1*x^(112).Í1-1124"x^2+111920*x^4+ o(x^6))+-G2*(1lx^(112)*!og(x)*(o(x^6))+'ll¡r(112).(1''ll8*x^2+11384*x^4+ O(x^6))),polynom)'))) 3

5

1/7920 (7920 _C1 x - 80 _C1 x + _C1 x + 1920

+

1og(x) o(x ) +

6 + 7920 _C2 o(x )) /

x O(x ) 2

6

192A _C2

5

_C1

1920

_c2 - 240 _c2 x +5

C2x

1/2 x

pretty(maple('dsolve({y(x).d iff(y(x),x$2)+diff(y(x),x)^2+ I =0, Y(0)=1,D(y)(0)=1 ),y(x),series)')) »»

CAPÍTULO 4. ECUACIONES DIFERENCIALES POR MÉTODOS APROXIMADOS 55

* u

y(x) :7+x-x

2

+x -3/2x

4

+5/2x

6

+ A(x

)

e

b

t¿l

»

pretty(s¡mple(maple('dsolve({diff(x(t),t$2)+diff(y(t),t)a*x+12=0, diff(y(t),t$2)-l0.diff(x(t),t)-y(t)+/=0,x(0)=1 ,Y(0)=1 ,D(x)(0)=l ,D(y)(0)=1), {x(t),y(t)},series)'))) 4 5 6 i 387 2 {y(t) :7 + t + 2 t + (- 43/2 + 20,/3 x) t - 3/2 t + (--- - 3 x) t + A(t ), x(t) :1 + t + (- 13/2 + 2 x) t

tr

- 2/3 t +(43/8 - 5/3 x) t

456

+ 3/1A t

+ O(t)]

»

pretty(simple(maple('dsolve({diff(x(t),t$2)+2*diff(x(t),t)+2*diff(y(t),t)+ 3*diff(z(t),t)+x(t)=1,diff(y(t),t)+diff(z(t),t)-x(t)=0,diff(x(t),t)+z(t)=g), {x(t),y(t),2(t)}, series)'))) {x(t) = x(0) + D(x) (0) t + (- D(x) (a) - 1/2 x(A) + 1/2) t 3

+ (1/2 D(x) (0) + 7/3 x(0) - 1/3) t + (- 1/6 D(x) (0) - 1/8 x(A) + 1/8) t

I

)

._:

56 + (1/24 D(x) (0) + 7/30 x(0) - 1/30) t + o(t ), z(t):

2

z(0) +x(0) t +7/2 D(x)(0) t + (- 7/3D(x)(0) - 1/6x(0) + 1/6) t

4+ ,4+

+ (7/8 D(x) (0) + 1/12 x(A) - 1/12) t

4

56 + (- 1/30 D(x) (0) - 1/40 x(0) + 1/40) t + o(t ), y(t):y(0)

2

+x(0) t+ 1/2 D(x)(0) t + (- 1/3 D(x)(0) - 7/6 x(0) + 1/6) t

+ (1/8 D(x) (0) + 1/12 x(0) - 7/72) t

4

+ (- 1/3A D(x) (0) - 1/40 x(0) + 7/40) t + a(t

6

)

56


vrÉro»o

4.3 EL

DE RUNGE-KUTTA

El método de Runge-Kutta permite obtener un conjunto de puntos de la solución aproximada de una ecuación diferencial. Maple aporta la opción numeric del comando maple('solve') que posibilita el cálculo de soluciones numéricas aproximadas de ecuaciones diferenciales. su sintaxis es:

maple('dsolve(ecuación, func(var),'numeric')

maple('f := dsolve({3*diff(y(x),x$2)^2 = diff(y(x),x$3)*diff(y(x),x), W''',o(yX0)=r,(o@@2XvX0)=1 ), v(x), numeric)') »»

f ::

proc (x)' dsolve/numeric/result 2'

(x, 3 B 7 9004,

I j | ) end

Ahora calculamos varios puntos de la función solución, para hacer una representación gráfica de la misma (figura 18-4).

»

[maple('f(-O.3)'),maple('f(-0.2)'),maple('f(-O.1)'),maple('f(0)'), maple('f(0. I )'), maple('f(0.2)'),maple('f(0.3)')l an§:

{x:

-.j,y(x)

:

.23s08893s9260396}{y(x)

.4045548849869109,

x:

.605s728090006967j{y(x) .3000400000000000j

- l){y(*)

:

»>

:

-

x: -.2} {y(x) 0} {x - .1, y(*) : .z} {y(x) : .A675444679682489, x

.3167840433732281,

.5000000000000000,

.7254033307s97474,

x:

x

-

-

y=1.2350889359260396,.31 67 8404337 32281,.4045548849869i 09,.S, .6055728090006967,.7254033307 597 47 4,.867 544467 9682a5il; » plot((-0.3:.1 :0.3),y)

CAPíTULO 4. ECUACIONES DIFERENCIALES POR MÉTODOS APROXIMADOS 57

u.9 0,a

tc

0.7

rla

0.6

de Su

0.4

0,5

D:

{:

O,?L

I lx),

nl

¡.1

t:1

03

Figura 4-1 Se observa que la figura es de tipo parabólico, con lo que podemoa ajustar la nube de puntos a una parábola. Dicha parábola será la solución simbólica aproximada de la ecuación diferencial. » pretty(vpa(poly2sym(polyfit((-O.3:.

1

:

0.3),y,2))))

2

.5747427827483573

x

+ 1.041293962469090 x + .4991457846921903

Ya tenemos la función polinómica y(x) solución de la ecuación D)),

D,-5,

Capítulo 5

SISTEIUAS DE ECUACIONES DIFERENCIALES Y ECUACION{ES EI\ DIFERENCIAS FIMTAS 5.1 SISTEMAS DE BCUACIONBS LINEALES HOMOCÉWNAS

COI{ COEFICM,NTES CONSTANTES Matlab resuelve directamente este tipo de sistemas sin más que utilizar el comando dsolve o maple('dsolve') con la sintaxis ya conocida Un sistema de ecuaciones diferenciales escrito como X'(t) = A X(t), presenta una solución general de la forma: )"i t

n

x(t)= T ciVie i=1

siendo i7,k) los autovalores correspondientes a los autovectores {Vk} de matriz A del sistema (k=1 ,2,...n) y todos los i,k distintos.

Si

existiese algun )"k =ak+bki, número complejo, genera

siguiente componente de la solución general:

t ¡.k t ck1 Wk1 e +Ck2Wk2e )"k

donde:

wk1 = % (vk + Vk)cos(bkt) + i¡2 (Vk - vk)sen(bk t)

la

60


Wk2 = itz

(ik-

Vk)Cos(bkt) - 1t2 (Vk + Vk)Sen(bk t)

Vk es el autovector correspondiente al áutovalor ik y Vk es su conjugado' si existe un fui de multiplicidad m>1, genera una parte de la solución general, de la forma: k }.i t 2 Xil ).it ),it ci

e Vi +c(i+1)te V(i+1)+c(i+2)t

e

V(i+Z¡+.'..+c(i+k)t

e

VK

» pretty(dsolve('Dx='5*x+3*y,Dy=-2*x'1 0*y','t')) y(t)

:

Ct exP( 7 0 +

*(t)

:

- 3/2 Cl exP(

C2 exP( 8 t),

7 r)

- C2 exP( 8 t)

También se puede usar la siguiente sintaxis:

»

pretty(maple('dsolve({diff(x(t),t)=-5*x(t)+3*y(t),diff(y(t),t)=-2*x(t}{x(t),Y(t)})'))

10.y(t)),

{y(t)

:

x(t)

:

-Cl - i/2

exP(

-Cl

7

t) +

-C2

exp( 7 t) -

exP( 8 t)'

-C2

exp(

I

t))

5.2 SISTEMAS DE ECUACIONES LINIEALES 1\O

HOMO'

CÉNSAS CON COEF'ICM,NTES CONSTANTES Ahora vamos a considerar sistemas de ecuaciones diferenciales no + homogéneas con coeficientes constantes de la forma X' = A X F(t)'

sea X = o(t)c la solución general del sistema homogéneo X',=AX-

Una solución particular del sistema no homogéneo es: -1

Xp =o(t)J o(t) F(t) dt

La solución general del sistema no homogéneo será X=O(t)c+Xp, con lo que tendremos la soluciÓn general:

CAPÍTULO 5. SISTEMAS DE ECUACIONES DIFERENCIALES 61

-1

X=

o(t)c

+ o(t)j o(t) F(0 dt

Este método es la generalización a sistemas de ecuaciones del método de variación de parámetros para las ecuaciones simples. Matlab ofrece directamente la solución de estos sistemas con el comando dsolve o maple('dsolve'), siempre y cuando pueda obtener la forma algebráica de las soluciones de las integrales que aparecen en la solución.

» pretty(simple(dsolve('Dx-Dy=s¡p(-t), Dy+§*¡+l*y=5i

¡(t+[)',

'x(0)=xs,r1g¡=yo','t'))) y(t)

:

7/50 sin(3) cos(t)

+ 7/50 cos(3) sin(t) + 5/6 exp( t)

+ l/50 sin(3) sin(t) - 1/50 cos(3) cos(t) - 5/7 +

5/7 yo - 5/7 xo

+ (- 7/50 sin(3) + 2/7 yo + 1/50 cos(3) - 5/42 + 5/7 xo) exp(7f)

5.3 ECUACIONI-ES EN DIFERENCIAS

FIMTAS

Matlab habilita la función maple('rsolve'), que permite resolver ecuaciones en diferencias finitas y recurrentes en general. Su sintaxis es:

maple('rsolve({ecuación, condiciones_iniciales}, función)') o (

» pretty(maple('rsolve({y(m+l )=m*y(m)+(m+1 )!,y(1 )=2},y)')) 1/2GAMaA(m) (n +m+2) p,

» pretty(maple('rsolve({y(2*n)={*y(n)+5,y(1 )=a},y)'))

62


loq (n) ______ +

2 2 log(2) a n + n (- 20/3 (1/4)

7

+ 5/3)

» pretty(maple('rsolve({y(n+2)-3*y(n+1)+!*y(n)=4^n,y(0)=1 ,y({ )=1},y)'))

+1/64

4/3-1/22

» pretty(maple('rsolve({x(n)-n*x(n)*x(n+l )=x(n+1 ),x(0)=l },x)')) 2

n+2 » maple('rsolve({x(n+2)-2*x(n+l )+5*¡(¡)=cos(3*n),x(0)=1,¡(l )=1 },x)') ans

:

*cos(l)-25*cos(n1 /2*(1 -2*ü/'tt+ 1 /2+(l +2*i)/Yt+(-39+cos(n-2)n3+22*cos(n-5)^3I 9*cos(n-4)"3 +32*cos(n-3)^3-25*cos(n-2)^3 2)^3*cos(j)+100*cos(n-2)"3*cos(1)n2+100*co,t(n-2)^3*cos(3)^2 +10*cos(n-3)/'3+cos(1)^2+10*cos(n-3)^3*cos(3)n2-72*cos(n3)n3 *cos(3)-7 2*cos (n-3)^3 *cos (l )-49*cos (n-4)^3 xcos( I ) + I 20*cos (n-4)^3 *cos (3)^2 + I 20*cos (n-1)"3*cos( I )^2-49*cos(n4)^3*cos(3) -60*cos(n-5)^3acos(3)-60+cos(n-5)^3*cos(l)+l69acos(n-2)^3+cos(l)*cos(3)-120+cos(n-2)^3 *cos(3)*cos(l)^2120*cos(n-2)n3*cos(3)^2*cos(1)-290*cos(n-3)^3*cos(3)"2*cos(l) +200*cos(n-3)^'3*cos(3)"2*cos(l)"2+208*cos(n3)^3*cos(1)acos(3)-290*cos(n-3)^3*cos(3)* cos(l)"2-100*cos(n-4)n3+cos(3)*cos(1)^2+265+cos(n-4)"3*cos(l)*cos(3)100*cos(n-4)"3 *cos(3)n2*cos(1)+50*cos(n-5)^3*cos(l)*cos(3)+100*cos(n-2)^3+cos(3)^2*cos(l)^2) /(2 5*cos(3)"2 + 36*cos (l )+cos(3)-3 0*cos (3)-3 0*cos(3)"24cos( I ) + 2 5 *cos ( l)"2-30*cos( I )30*cos(3)*cos(1)"2+25+25*cos(3)^2*cos(l)^2)+l/4+(-30*i*(l+2*r^n-35+(l-2*i)^n*cos(3)^2 +27a(l-2*i)^n*cos(1)-25*(12*i)^n*cos(l)^2-3 5 *(l +2*i)^nacos(3)^2-37+i*(l -2*i)^n*cos(3)-25*(l +2*i)^n*cos(l)^'2+ 2 I *(l-2*i)^n*cos(3)^2*cos(1 )+23*(1 2*i)^n*cos(3)*cos(112-20*(1 -2+i)^nacos(3)^Z*cos(1)^2+25*i+( I -2*i)^n*cos(3)^2+ I 5 *i+(l -2*i)^n*cos(l)^2-20*(l +2*r^n *cos(3)^2*cos(l )"2- I 5*¡*(l +2*i)^n*cos(1)"2+374i*(l +2+i)t'n*cos(3)+39*i*(l + 2*i)^n*cos(l)-t9*¡+( I 2*i)nn*cos(l )+ 2 I * (1 + 2*i)^n*cos(3)^ 2*cos( I )+4 ¡ +¡ 1 ¡ 2 *i)^n*cos (3)+ 2 7 + (l +2 *i)^n *cos( I ) +48*i*( l 2+i))t*cos(l )*cos(3)+ I 0*i*(l -2*i)^n*cos(3)^2*cos(l )^2-25*i*(l +2*i)^n*cos(3)"2+ 19*i*(l +2+i)^n*cos(3)*cos(l )^248*i*(l+2*i)nnacos(1)*co.t(3) +33*i*(l+2*i)"n*cos(3)^2+cos(l)-24*(1+2*i)^n*cos(l)*cos(3)+23+(l+2*i)^n*cos(3)*cos(l)^240*(1 -2*i)^rt-40*(1 +2*i)^n+4 I *(1 -2*i)^nacos(3)-l 0*i*(1 +2*i)^n*cos(3)^2*cos(1)"2-3 3*i+( I -2*i)"n*cos(3)^2*cos(t )+30*i*(1 2*i)^n-24*(l-2*r^n+cos(l)*cos(3)-19*i*(1-2*í)"n*cos(3) *cos(l)"2)/(5*cos(l)"2-6*cos(l)+5)/(5-6+cos(3)+5*cos(3)^2)-3/2+C 2*cos(n-2)*cos(l )-2*cos(n-2)'r I O*cct.s(n-2)*cos(1 )^2-5*cos(n-3)*cos(l )+cos(n-3))/(5*cos(l )^2-6*cos(1 )+ 5)-3/8*(í*(1 2*i)^n+cos(l )"2-2*(l -2*i)^n*cos(l )^2-4*i*(l -2*i)^n*cos(l)+ 3*i*(l -2*i)^n-4*(1 -2*i)^n+2*(l -2*i)^n*cos(l )2*(1 +2*i)^n*cos(1 )^2-¡*(l +2*i)^n*cos(1 )^2+4*i*(l +2*i)^n+cos(l )-4*(l +2*i)^n3 *i* ( I + 2+ ¡)/'n+ 2*( I + 2*i)^n+cos(1 ))/(5 *cos (l )^ 2-6*co,s( I )+ 5)

Ahora intentamos simplificar el resultado no trivial.

»

pretty(simple(maple('evalf(rsolve({x(n +2)-2*x(n+1)+5*¡1¡¡=cos(3*n),

x(0)=1,¡11 ¡=1 ),x))'))) .4373424522418379 (1. - 2. i)

+ .4373424522418384 (1. + 2. i)

-

CAPÍTULO 5. SISTEMAS DE ECUACIONES DIFERENCIALES 63

- .2069717398607332 cos(n)

9418421 760747166 cos (n)

sin

(n)

2

4405817940094277

3

- .7714545418874307 sin(n)

cos(n) sin(n) n

.0626567574570487 .

i (1.

8765270805583898 cos (n)

+

2.

i)

t .0626567574514487 i (1. - 2. i)

180378847j807544 sin(n)

Capítulo 6

CALCULO I{LMERICO EN MATLAB. APLICACIONES A LAS ECUACTONES DIFERENCIALES 6.1

MATLAB Y LA PROGRAMACION

MATLAB puede utilizarse como un lenguaje de programación de alto que nivel incluye estructuras de datos, funciones, instrucciones de control de flujo, manejo de entradas/salidas e incluso programación orientada a objetos.

Los programas de MATLAB suelen escribirse en ficheros denominados M-ficheros. Un M-fichero no es más que código MATLAB (scripts) que simplemente ejecuta una serie de comandos o funciones que aceptan argumentos y producen una salida. Los M-ficheros se crean utilizando el editor de texto, tal y como ya se vio en el Capítulo 2. 6.2 BDITOR DE TEXTO Para crear un nuevo M-fichero se utiliza el Editor/Debugger, que se activa haciendo clic en el botón D de la barra de herramientas de MATLAB o mediante File+ New-+ M-file en el escritorio de MATLAB (Figura 6-1) o en la ventana de comandos (Figura 6-2). El Editor/Debugger se abre con un fichero en blanco en el cual crearemos el M-fichero, es decir, un fichero con código de programación MATLAB (Figura 6-3). Cuando se trate de

66


abrir un M-fichero ya existente se utiliza File-+ Open en el escritorio de MATLAB (Figura 6-1) o, alternativamente, se puede utilizar el comando Open en la ventana de comandos (Figura 6-2). También se puede abrir el Editor/Debugger haciendo clic con el botón derecho del ratón en el interior de la ventana Current Directory y eligiendo New -+ M-file en el menú emergente resultante (Figura 6-4). La opción Open de este menú abre un M-fichero existente. Se pueden abrir varios M-ficheros simultáneamente, en cuyo caso aparecerán en ventanas diferentes.

I Dr\...12\worklmatri¿1,trc 2 Dr\,..ltinnnEelBtrubond.m

I

D:\,..\Finan.elÉsud¡sr,m

4

Dit.,,llinihEeldmrtizÉ,m

Figura 6-l

Figura 6-2

w fr-#rrumm fÍ * couuicarions ror ú-$cnnt.ar stscet rdiW ii;$nuo aE$ieil.iutr 'j lenffie +'],out*.e"

otld' ToúIbo* :

r*á#W;*itiJiiit

Figura 6-3

: tu

Figura 6-4

CAPÍULO 6 CALCULO NUMERICO EN MATLAB. APLICACIONES A LAS ECUACIONES... 67

La Figura 6-5 muestra las funciones de los iconos del Editor/Debugger. ,{}ru ruete

trf-áehe*o

Ahrir EÉ-$iah*ro

v§.

*ub sxiitsat*

s,,rqqfi+zg:+d-tei*s

ee r¡n S;[-§s&ere

A X!*i& EÉ&et*¡l

É¡¡*ráer&á-€i*hsrc Empriroir l{-ñe&.ero

+tusr pr:uto d* esrt* p*re de¡*lggs flrec'fulerfttJ paato de mrt* IEI;"'li*n Fe:Et*.{**§et§ I F1*m*u liaee e*Ml d+ h{-Éiehero I I &so d* eiecueioa e Esá fuaeióe

I

I 'I

.Peu¡aáe¡*reaF;ireiar 1Ei**raÉhs**nrs*errarÉmrt*

I I I t

=i*;'ta;

!* ir

i

#*\.e@*

**6:1*;'

;

Figura 6-5

6.3 SCRIPTS Los scripts son el tipo de M-fichero más sencillo posible. Un script no tiene argumentos de entrada ni de salida. sencillamente está formado por instrucciones MATLAB que se ejecutan secuencialmente y que podrían submitirse igualmente en serie en la ventana de comandos. Los scripts operan con datos existentes en el espacio de trabajo o con nuevos datos creados por el propio script. cualquier variable que se cree mediante un script permanecerá en el espacio de trabajo y podrá utilizarse en cálculos posteriores después de finalizar el script.

A continuación se presenta un ejemplo de script que genera varias curvas en polares representando pétalos de flores. Una vez escrita la sintaxis del script en el editor (Figura 6-6), se guarda en la librería de trabajo (work) y simultáneamente se ejecuta, haciendo clic en er botón JE o utilizando la opción save and run del menú Debug (o presionando FS). Para pasar de un gráfico al siguiente basta con pulsar ENTER.

Figura 6-6

68


Figra6-7

Figura

6-9

Figura 6-8

Figura 6-10

6.4 FUNCIONES Y M.FICHEROS. FUI{CTION,

BVAL Y

FEVAL Ya sabemos que MATLAB dispone de gran variedad de funciones para usar en el trabajo cotidiano con el programa. Pero, además, el programa ofrece la posibilidad de definir funciones a medida. El camino más usual para definir una función a medida es escribir su definición en un fichero texto, denominado M-fichero, que será permanente y que, por lo tanto, permitirá el uso posterior de la función siempre que se requiera.

MATLAB es habitualmente utilizado en modo comando (o interactivo), en cuyo caso se submite un comando que se escribe en una

CAPíTULO 6 CALCULO NUMERICO EN MATLAB. APLICACIONES A LAS ECUACIONES... 69

única línea sobre la ventana de comandos y se procesa de inmediato. pero MATLAB también permite la ejecución de conjuntos de comandos en modo batch, en cuyo caso se submiten secuencialmente un conjunto de comandos escritos previamente en un fichero. Este fichero (M-fichero) ha de ser almacenado en disco con la extensión ".m" en el camino de subdirectorios de MATLAB, utilizando cualquier editor ASCII o Ia subopción M-file de Ia opción New del menú File de ra barra superior de menús, la cual nos lleva a un editor de texto que permitirá escribir las líneas de comandos y guardar el fichero con un determinado nombre. La opción open M-File del menú File de la barra superior de menús permite editar cualquier M-fichero preexistente.

Para ejecutar un M-fichero basta con teclear su nombre (sin extensión) en modo interactivo sobre la ventana de comandos y pulsar Enter. MATLAB interpreta secuencialmente todos los comandos o sentencias incluidos en las diferentes líneas del M-fichero y los ejecuta. Normalmente no aparecen en pantalla los riterales de los comandós que MATLAB va interpretando, salvo que se active er comando echo on, y sólo se van viendo los resultados de las ejecuciones sucesivas de los comandos interpretados. Normalmente, el trabajo en modo batch es útil cuando se procesan conjuntos muy largos de comandos de escritura tediosa y con propensión a cometer errores, pero la mayor utilidad se presenta en Ia automatización de procesos. Además, en las líneas de un M-fichero se pueden introducir textos explicativos y comentarios, empezando cada línea al efecto por el símbolo %. con el comando Help se accede al texto explicativo de un M-fichero. El comando function permite la definición de funciones a medida en MATLAB, constituyendo una de la aplicaciones más útiles de los Mficheros. La sintaxis de este comando es la siguiente:

nction pa rámetros_sa ida = nom bre_fu nción (pa rámetros_e ntrada ) cuerpo de la función

fu

I

una vez que la función ha sido definida, se guarda en un M-fichero para su uso posterior. Es Útil también introducir algún texto explicativo en la sintaxis de la función (entre o/o), al cual se accederá con el comando de ayuda Help. cuando los parámetros de salida son más de uno, se sitúan entre corchetes y separados por comas. Si los parámetros de entrada son más de uno, se separan por comas. El cuerpo de la función es la sintaxis que la

7O


define, y debe incluir comandos o instrucciones que asignen valores a los parámetros de salida. Cada comando o instrucción del cuerpo suele ir en una línea que finaliza con una coma o con un punto y coma en caso de que se definan variables (para evitar repeticiones en las salidas al ejecutar la función). La función se guarda en el M-fichero de nombre nombre_función.m.

A continuación vamos a definir la función funl(x)=¡n3-2x+cosx, creando el correspondiente M-fichero de nombre fun1.m. Para introducir dicha sintaxis en MATLAB se selecciona la subopciÓn M-file de la opciÓn New del menú File de la barra superior de menús (o se hace clic en el botón D Oe la barra de herramientas de MATLAB), la cual nos lleva al editor de texto MATLAB Editor/Debugger que permitirá escribir las líneas de comandos relativas a la sintaxis de la funciÓn, tal y como se indica en la Figura 6-11.

?.T:.rlt i.:iúe áe s1B fu.rié* p=x^3-:*x+cqE (xi ¡

Bieplt

;;::;t: Figura 6-11

Para guardar de forma definitiva dicha sintaxis en MATLAB se selecciona la opción Save del menú File de la barra superior de menús del MATLAB Editor/Debugger, la cual nos lleva a la caja de diálogo Guardar de la Figura 6-12, mediante la cual podemos guardar nuestra funciÓn con el nombre deseado y en el subdirectorio que se indique como ruta en el campo Nombre de archivo. También puede hacerse clic en el botón 'lü o utilizar la opción Save and run del menú Debug. Se recomienda guardar las funciones como ficheros de nombre igual al nombre de la función y en el subdirectorio de trabajo por defecto de MATLAB C :\MATLAB6pI\work. Precisamente es la situación que presenta el programa por defecto al intentar guardar cualquier función a medida y al utilizar el icono JE o la opción Save and run del menú Debug para ejecutar y guardar simultáneamente el M-fichero actual del editor.

cnpílulo

6

cAlcuLo NUMERtco

EN MATLAB. ApLtcActoNEs A LAS

EcuActoNEs... 71

Figurra6-12

Una vez definida y guardada la función anterior en un M-fichero, se puede utilizar desde la ventana de comandos. Por ejemplo, para hallar el valor de la función en 3nl2 escribimos sobre la ventana de comandos usual de MATLAB:

) funl (3*pí/2) ans e

95.2274

atr

le el

ol

Para pedir ayuda sobre la función anterior (suponiendo que fue previamente introducida como comentario al elaborar el M-fichero que Ia define) se utiliza el comando help, como sigue:

o

>) heJ-p funl Definición de una función simple

al

La evaluación de una función en sus argumentos (o parámetros de entrada) también puede realizarse a través del comando feval, cuya sintaxis es la siguiente:

le ta al ar

Evalúa lafunción F (M-fichero F.m) en los argumentos especificados argl , arg2,...,argn

72


Como ejemplo construimos un M-fichero de nombre ecuacion2.m que contiene la función ecuacion2, cuyos argumentos son los tres coeficientes de la ecuaciÓn de segundo grado ax2+bx+c - 0 y cuyas salidas son sus dos soluciones (Figura 6-13)'

Si ahora queremos resolver la ecuación x2+2x+3 = 0 utilizando feval' escribimos en la ventana de comandos lo siguiente:

>> [x1 rx2]=¡sr"J-('ecuaciorr2',1,2,31

x7: -1.0000 + 7.4742i

x2: -1.0000 - 1.4142i También se Puede resolver la ecuación de segundo grado como sigue:

>> [x1 ,x2] =ss¡r"cion2 (l ,2,31

x7: -1.0000 + 1.4742i .-a

-

-1.0000 - 7.4742i

CAPíIULO 6 CALCULO NUMERICO EN MATLAB. APLICACIONES A LAs ECUACIONES... 73

Si pedimos ayuda sobre la función ecuación2 tenemos lo siguiente:

>> help ecuacion2

Esta

función

resuefve

fa

ecuación

de

segundo qrado

ax^2+bx+c:0

cuyos coeficientes son a, b y c (parámetros de entrada) y cuyas sofuciones son x7 y x2 (parámetros de salida) La evaluación de una función en sus argumentos (o parámetros de entrada) cuando son cadenas se realiza a través del comando eval, cuya sintaxis es la siguiente: eval(expresión)

Ejecuta la expresión cuando es una cadena

Como ejemplo evaluamos una expresión carácter que define matiz mágica de orden 4.

la

>> n=4; )) eval ( ['M' num2str (n) ' = magic (n) ' ] )

M4: 162313 511108 97672 474757 6.5

VARIABLES LOCALES Y GLOBALES

Normalmente, cada función de MATLAB definida como un M-fichero contiene sus variables como variables locales, es decir, como variables que tienen su efecto en el inter¡or del M-fichero, separadamente de otros M-ficheros y del espacio de trabajo base. No obstante, es pos¡ble definir var¡ables en el interior de M-ficheros de modo que tengan efecto simultáneamente en el interior de otros M-ficheros y en el propio espacio de trabajo base. Para ello es preciso definir las variables como globales con el comando GLOBAL, cuya sintaxis es la siguiente: GLOBAL xy 2...

Define las vqriables x, y, 2... como globales

Cualquier variable definida como global en el interior de una función

es accesible separadamente para el resto de las funciones y para

el

espacio de trabajo base línea de comandos. Si la variable global no existe, la primera vez que se define en la sentencia GLOBAL, se inicializará como

74


la matriz vacía. Si ya existe una variable con el mismo nombre que la que se está definiendo como global, MATLAB emite un mensaje de peligro y cambia el valor de la variable a enlazar como global. Es conveniente declarar una variable como global en cada función que necesite acceso a ella, y también en la línea de comandos, para tener acceso a ella desde el espacio de trabajo base. En el interior de las funciones el comando GLOBAL se sitúa al principio (antes de cualquier ocurrencia de la variable). Como ejemplo, supongamos que se quiere estudiar el efecto de los coeficientes de interacción cr y B en el modelo de Lotka-Volterra:

t1 = 3'1-a}r3t,

i

z=

-j

z+

P-u1J2

Para ello creamos la función lotka en el M-fichero lotka.m de la Figura 6-14.

Figura 6-14

Si posteriormente sobre la línea de comandos escribimos

lo

siguiente:

)) global ALPIA BETA ALPHA = 0.01 BETA = 0.02

ya se podrán utilizar estos valores globales para o y B en el interior del M-fichero lotka.m (sin tener que volver a especificarlos). por ejempro, podrá realizarse una representación (Figura 6-15) con la sintaxis siguiente:

>> [t,y] = ode23('lotka',0,10, [1; 1]); plot(t,y¡

cAPrruLo 6 CALCULO NUMERTCO EN MATLAB.

APLTCACTONES

A LAS ECUACTONES... 7s

Figura 6-15

6.6 TIPOS DE DATOS En MATLAB hay 14 tipos de datos diferentes que se representan en la Figura 6-'16. ARFIAY

funetion handle

int8. uint§, intl §, uintl Sn

srngJ-e

doubl-e I

int33, uint32

Epar§É

Figura 6-16

A continuación se detallan los distintos tipos de datos: Tioo de dato

Eiemplo

single

Descripción

3*10^38

Precisión numérica simple. Requiere menor almacenamiento que la doble precisión, pero es menor. Este tipo de datos no debe usarse en operaciones

double

3*10^300 5+6i

Doble precisión numérica. Es el tipo de variable más común en MATLAB

matemáticas.

76


speye (5)

sparse

int8, uint8, int16, uint16, int32. uint32

uintS (¡nagic (3) )

Matriz disoersa con doble precisión. Enteros con y sin signo y con 8, 16, y 32

bits de longitud. Posibilitan usar cantidades enteras con manejo eficienfe de memoria. Este tipo de datos no debe us

ars e en ooeraciones matemátic as.

char

'HeIIo'

Caracteres (cada carácter tiene una

ce11

{17 'heI1o' ef'e(2)}

lonsitud de 16 bits). Celda kontiene datos de similar tamaño)

structure user class java class

a.daY = L2 a.color = 'Red'

a.mat = maqic(3) inline ('sin (x) ') java. awt. Erame

Esfructura (contiene celdas de similar tamaño) Clase MATLAB (creada con funciones)

Clase

Jwa

(defrnida en API o propia con

.Iava)

function

@humps

Maneja funciones MATLAB. Puede ser

pasada en una lista de argumentos

handle

evaluada con

Y

fevaf.

6.7 CONTROL DE FLUJO: BUCLES FOR, ELSETF

WHILE E IF

El uso de funciones recursivas, cond¡cionales y definidas a trozos

es muy habitual en matemáticas. Para la definición de este tipo

de natural, la definiciÓn eS Como manejo de bucles. el necesario eS funciones de las funciones se hará a través de M-ficheros.

6.7.1El bucle FOR MATLAB dispone de su propia versión de la sentencia DO (definida en la sintaxis de la mayoría de los lenguajes de programaciÓn). Esta sentencia permite ejecutar de forma repetitiva un comando o grupo de comandos varias veces. Por ejemplo:

': *r

i=1:3, x(i)=O,

end

CAPÍIULO 6 CALCULO NUMERICO EN MATLAB. APLICACIONES A LAS ECUACIONES... 77

La forma general de un bucle FOR es la siguiente:

for

variable

:

expresión

comandos

end

El bucle siempre empieza con la cláusula for y termina con

la

cláusula end, e incluye en su interior todo un conjunto de comandos que se separan por comas. Si algún comando define una variable, se finaliza con punto y coma para evitar repeticiones en la salida. Normalmente, los bucles se utilizan en la sintaxis de M-ficheros. Veamos un ejemplo (Figura 6-17):

Figtra 6-17 En este bucle hemos definido la matriz de Hilbert de orden (m,n). Si su contenido como un M-fichero de nombre matriz.m, podremos construir cualquier matriz de Hilbert posteriormente ejecutando

guardamos

el

M-fichero

y

especificando los valores para las variables

(dimensiones de la matriz) como se indica a continuación:

)) M=matríz(4,51 M7-

0000

0.5000 0.3333 0.2500

0.

s000

0.3333 0.2500 0.2000

0.3333 0 .2s00 0.2000 0.1567

0-2500 0.2000 0.1667 0.1429

0-2000 0.1667 0.7429 0.7250

my

n

78


6.7.2 El bucle

WHILE

MATLAB dispone de su propia versión de la sentencia WHILE definida en la sintaxis de la mayoría de los lenguajes de programación. Esta sentencia permite ejecutar de forma repetitiva un comando o grupo de comandos un número determinado de veces mientras se cumple una condición lógica especificada. La sintaxis general de este bucle es la siguiente: while condición comandos

end

El bucle siempre empieza con la cláusula while seguida de una condición, y termina con la cláusula end, e incluye en su interior todo un conjunto de comandos que se separan por comas y que se ejecutan mientras se cumple la condición. Si algún comando define una variable, se finaliza con punto y coma para evitar repeticiones en la salida. Como ejemplo escribimos un M-fichero (Figura 6-18) que se guarda con el nombre while1.m, cuya ejecución permite calcular el mayor número cuyo factorial no excede a 10100.

Figura 6-18

Ahora ejecutamos el M-fichero.

)

whi1e1

70

CAPÍTULO 6 CALCULO NUMERICO EN MATLAB. APLICACIONES A LAS ECUACIONES... 79

6.7.3 El bucle IF ELSEIF ELSE END MATLAB, al igual que la mayoría de los lenguajes de programación

estructurada, también incorpora la estructura lF-ELSEIF-ELSE-END. Mediante esta estructura, se pueden ejecutar secuencias de comandos si se cumplen determinadas condiciones. La sintaxis del bucle es la siguiente: if condición comandos

end

En este caso se ejecutan los comandos si la condición es cierta. Pero la sintaxis de este bucle puede ser más general. if condición comandosl else comandos2

end

En este caso se ejecutan los comandosl si la condición es cierta, y se ejecutan los comandos2 si la condición es falsa.

Las sentencias lF, al igual que las sentencias FOR, pueden ser anidadas. Cuando se anidan varias sentencias lF se utiliza la sentencia ELSEIF, cuya sintaxis general es la siguiente: if condiciónl comandosl elseif condición2 comandos2

elseif condición3 comandos3 'else

end

En este caso se ejecutan los comandosl si condición1 es cierta, se ejecutan los comandos2 si condición1 es falsa y condición2 es cierta, se ejecutan los comandos3 si condición'l y condición2 son falsas y condición3 es cierta, y así sucesivamente.

80


La sintaxis anidada anterior es equivalente, pero más rápida de ejecución, a la sintaxis sin anidar siguiente: if condiciónl comandosl else

if condición2 comandos2 else

if condición3 comandos3 else

end end end

Veamos como ejemplo (Figura 6-19) elsel.m siguiente:

el

M-fichero

de

nombre

Figura 6-19

Al ejecutarlo obtendremos el tipo de número (negativo, par o impar) para un valor de n especificado:

)) else1 (8) , elsel (5) , else1 (-10) n es par

CAPÍTUIO 6 CALCULO NUMERICO EN MATLAB. APLICACIONES A LAS ECUACIONES... 81

n es impar A: n es neqativo 6.7.4 SWITCH Y CASE La instrucción switch ejecuta ciertas sentencias basadas en el valor de una variable o expresión. su sintaxis básica es la siguiente: switch expresión (escalar o cadena) case

valorl

sentencias case

% Ejecttasi expresión esvalorl

valor2

sentencias '% Ejectta si expresión es valor2

otherwise sentenciass %o Ejeatta si expression no cumple %o

ningún caso

end continuaciÓn se Presenta un ejemplo de una función que devuelve -1, 0, 1 u otro número según sea la entrada (Figura 6-20).

A

Figura 6-20

82


Después ejecutamos el ejemplo anterior.

>> casel (25)

otro vafor

>> casel (-1) menos uno

6.7.5 COI{TINUE La instrucción continue pasa el control a la iteración s¡guiente en un bucle for o while en el cual aparece ignorando las restantes instrucciones en el cuerpo del bucle. A continuación se muestra el M-fichero continue.m (Figura 6-21) que cuenta las líneas de cÓdigo en el fichero magig.m ignorando las líneas blancas y los comentarios.

Figura6-21

Seguidamente ejecutamos el M-fichero.

)) continuel 25 fineas

6.7.6 BREAK La instrucción break finaliza la ejecución de un bucle for o while en el cual aparece continuando la ejecución en la siguiente instrucción fuera del bucle. A continuación se muestra el M-fichero breakl.m (Figura 6-22) que lee las líneas de código en el fichero fft.m saliendo del bucle cuando se encuentre Ia primera línea vacía.


Figra6-22 Después ejecutamos el M-fichero.

>> breakl

Discrete Fourier transforn. is the discrete Fourier transform (DFT) of vector X. For matrices, the FFT operation is appTied to each column. For N-D arrays, the FtrT operation operates on the fitst non-singleton

';FFT

2; Z * Z 6

* 3 z

Z 3

FFT(X)

dlmension.

is the N-po1nt FFT, padded with zeros 1f X ¡as than N points and truncated if it has more.

FFT(X,N)

-¿ess

FFT(X,[],DIM) or trFT(X,N,DIM) applies the FFT operation across the dlmension DIM.

z

& tror Tenqth N input vector x, the DET is a length N vector X, Z with efements ZN X(k) : sum x(n)*expr-i*2*pi*(k-1)*(n-7)/N), 1 <:k <:N. Z n:7 & * The inverse DFT (conpüted by IFFT) is qiven by ZN x(n) : (1/N) sun X(k)*exp( j*2*pi*(¡<-1)*¡n-1)/N), 7 (: n <: N. & k:1 Z Z

*

See

afso IFFT, FFT2. IFtrT2, FETSHIFT.

84

ECUACIONES DIFERENCIALES.CON MATLAB

6.7.7 TRY ... CATCH no ocurra un Las instrucciones entre try y catch se ejecutan mientras del error' La error. La instrucción lasterr sá utiliza para ver loa causa sintaxis general de la instrucción es la siguiente: tty, instrucción, instrucción, catch,

instrucción, ..

.,

instrucción. end

6.7.8 RETURN y La instrucción return finaliza la secuencia actual de comandosse continuación A devuelve el control a la funciÓn invocada o al teclado' pl""'"ntu un ejemplo (Figura 6-23) que calcula el determinante de una 1' matizque no es vacía. Si-la matriz es vacía Se obt¡ene el valor

Figura6-23

vacía' Posteriormente ejecutamos la funciÓn para una matriz no

¡¡ [=[-1,-1,! ;L,0,L;L,L,Ll A:

CAPÍTLO 6 CALCULO NUMERICO EN MATLAB. APLICACIONES A LAS ECUACIONES... 85

-7-71 707 117 >> detl (A) ans 2

Ahora aplicamos la función a una malriz vacía.

>>

F[]

D-

tl >> detl (B) ans : .L

6.8 SUBFUNCIONES Las funciones definidas mediante M-ficheros pueden contener código para más de una función. La función principal en el M-fichero se denomina función primaria, que es precisamente la función que invoca el M-fichero. Pero adicionalmente puede haber subfunciones colgando de la función primaria y que sólo son visibles para dicha función primaria o para otra subfunción dentro del mismo M-fichero. Cada subfunción comienza con su propia línea de definición de función. En la Figura 6-24 se presenta un ejemplo.

Figtra6-24

86

ECUACIONES DIFERENCIALESCON MATLAB

Las subfunciones mean y median calculan la media y la mediana de la lista de entrada. La función primaria newstats determina la longitud de la lista y llama a las subfunciones pasándoles la lista de longitud n. A la hora de ejecutar la función principal, basta darle su entrada correspondiente (una lista de valores para los que se calculará su media y su mediana) y las llamadas a las subfunciones Se realizan automáticamente tal y como se ve a continuación.

>> [media, mediana] = neytstats media :

( [10

,20,3,4,5,67)

I mediana

:

5.5000

6.9 ECUACIONES DIFERENCIALES ORDINARIAS

MEDIANTE MÉTODOS DE CÁLCULO I{UMÉRICO Obtener soluciones exactas de ecuaciones diferenciales ordinarias

no eS una tarea sencilla. Existen diferentes métodos para obtener resultados numéricos aproximados para las soluciones de ecuaciones diferenciales ordinarias, entre los que destacan el método de Euler, el método de Heun, el método de las series de Taylor, los métodos de Runge-Kutta (implementados en el mÓdulo básico de MATLAB), el método Adams-Bashforth-Moulton, el método de Milne y el método de Hamming.

de

6.9.1Método de Euler Sea [a,b] el intervalo en el que queremos resolver la ecuaciÓn diferencial y'=f(t,y) con y(a)=y0. Dividimos el intervalo [a,b] en M subintervalos del mismo tamaño usando la partición dada por los puntos tk ¡¡=(b-ayM. El método de Euler obtiene las =a+kh k=0,1, ..., M "on aproximaciones a la solución de la ecuaciÓn diferencial mediante la iteración yk+1 = yk+hf(tk, yk) k=0,1, ..., M-1. El método de Euler puede implementarse mediante el M-fichero de

laFigura6-24.

i I

I

l L


Figwa6-24

6.9.2 Método de Heun Sea [a,b] el intervalo en el que queremos resolver la ecuación diferencial y'=f(t,y) con y(a)=yO. Dividimos el intervalo [a,b] en M subintervalos del mismo tamaño usando la partición dada por los puntos tk =a+kh k=0,1, ..., M con h=(b-a)/M. El método de Heun obtiene las aproximaciones a la solución de la ecuación diferencial mediante la iteraciÓn yk+1 = yk+h(f(tk, yk)+ f(tk+1, y[+ f(tk, yk)))/2 k=0,'1, ..., M-1.

El método de Heun puede implementarse mediante el M-fichero de la Figura 6-25.

88

ECUACIONES DIFERENCIALESCON MATLAB

Figlora6-25

6.9.3 Método de las series de Taylor

sea [a,b] el intervalo en el que queremos resolver la ecuaciÓn diferencial y'=f(i,y) con y(a)=y0. Dividimos el intervalo [a,b] en M

subintervaloi Oel mismo tamaño usando la partición dada por los puntos tk =a+kh k=0,1, ..., M con ¡=(b-a)/M. El método de las series de Taylor (vamos a considerarlo de orden 4) obtiene las aproximaciones a la solución de la ecuaciÓn diferencial evaluando y', Y", y"' e y"" mediante sus respectivos desarrollos en series de Taylor de orden 4. El método de las series de Taylor puede implementarse mediante el M-fichero de la Figura6-26.

t-

CAPíTULO 6 CALCULO NUMERICO EN MATLAB. APLICACIONES A LAS ECUACIONES.., 89

Figxa6-26 Como ejemplo resolvemos la ecuación diferencial y'(t) = (t-y)12 en y(0)=1 mediante los métodos de Euler, Heun y series de Taylor. con [0,3] Comenzaremos definiendo la función f(t,y) mediante el M-fichero de la Figura 6-27 .

Figxa6-27 La solución de la ecuación por el método de Euler en 100 pasos se calculará como sigue:

)) E=euler('dif1',O,3,1,100)

90


E_ . 0. 0. 0

0 0 06000000000000 09000000000000 03 00000000000

1.00000000000000 0.98500000000000 0.97067s00000000 0.95707487500000

0.12400000000000 0.15000000000000 0.18000000000000

0. 944009 6s1B7500

,.

7.56377799005970 1 .58307732020827 7. 60252525040509 7. 62213737764907 7.64790537107428 1.66182673140876

ttooooo0000000 2.88000000000000 2.91000000000000 2.94000000000000 2.97004000000000 3.00000000000000

0

.

9316495

0

7096B

I

0.97992476449042

Esta solución puede graficarse (Figura 6-28) como sigue:

>> plot

(E ( :

,2)

)

Figura 6-28

La solución de la ecuación por el método de Heun en 100 pasos se calculará como sigue:

CAPÍruLO 6 CALCULO NUMERICO EN MATLAB. APLICACIONES A LAS ECUACIONES... 91

)) H=heun('dif1' r0r3r1,100) n0 0 060a0000000000 0.09000000000000

. 0. 0

03 00000000000

0.12000000000000

2.88000000000000 2.91000000000000 2.94000000000000 2.97000000000000 3.00000000000000

1.00000000000000 0.9

B53

3 750000000

0.97133997296875 0.95799734007443 0 - 94530002967496

1 . 590822 093 794 54

7.67023972987327 7.62987497089478 7.64954529740884 1 . 66942 I 5 60I B2 99

La solución mediante el método de las series de Taylor exige definir previamente la función df = [y' y" y"' y""] mediante el M-fichero de la Figura 6-29.

Figtura6-29

La ecuación diferencial se resolverá por Taylor mediante la sintaxis:

>> T=taylor ( 'df ' ,0 ,3, 1 ,100)

0

.

03 00000000000

0 0

0.06000000000000 0

.

0900000000000

0

1.00000000000000

0.98533581882813

0.971i366006828i 0 . 957 9924455544i 0-

88000000000000 2.91000000000000 2-

94529360082516

1.59078327648360 1.61020109213866

92

ECUACIONES DIFERENCIALE§ CON MATLAB

2.94000000000000 2.97000000000000 3.00000000000000

1. 62977645599332 7. 64950702245046 1. 66939048087422

Comenzamos definiendo la función f(t,y) mediante el M-fichero de la Figura 6-30.

Figura 6-30

A continuación resolvemos la ecuación diferencial por el método de Euler en 20 pasos mediante la sintaxis siguiente:

)>

E=eu]-er

('díf2 ' ,0,0

.8 ,L ,201

L-

0

1.00000000000000

0.04000000000000 0.08000000000000 0.12000000000000 0.16000000000000 0.20000000000000 0.2400000000a000 0 .28000000000000 0. s2000000000000 0.36000000000000 0.40000000000000 0.44400000000000 0.48000000000000 0.52000000000a00 0.56000000000000 0.60000000000000

1.04000000000000 1.08332800000000 7.73052798222336 7.78222772296696 1.23915821852503 7.30217874274655 7.372309s2120649 7.45077485808625 1.53906076564045 1.63899308725380 7 .75284502085643 7. 8834B764 7542 0

I

2 - 03460467627982 2 .27100532382941

2.419097105s0949

CAPÍIULO 6 CALCULO NUMERICO EN MATLAB. APLICACIONES A LAS ECUACIONES.,. 93

0. 64000000000000 0.68000000000000 0. 72000000000000

0.76000000000000 0.80000000000000

2.66757777657970 2.96859267586445 3.339590s006230s s.8064408s566367 4.4091045090

79 9

9

La solución puede graficarse (Figura 6-31) como sigue:

plot

(E

(: ,21I

Figura 6-31

Capítulo 7

ECUACIONES EN DIFERENCIAS CON VALORES IMCIALES, VALORES EN LA FRONTERA Y EN DERIVADAS PARCIALES 7.1 SOLUCIÓN NUMÉRICA DE ECUACIONES

DIFERENCTALES MATLAB dispone de comandos en su módulo básico que permiten

resolver numéricamente ecuaciones diferenciales ordinarias (ODEs), ecuaciones diferenciales algebraicas (DAEs) y resolución de problemas con valores en la frontera oDEs. También es posible la resolución de sistemas de ecuaciones diferenciales con valores en la frontera y ecuaciones diferenciales en derivadas parciales parabólicas y elípticas. 7

.2 ECU ACIONES DTFERENCIALBS ORDINARIAS

CON VALORES INICIALES Una ecuación diferencial ordinaria contiene una o más derivadas de

la variable dependiente y respecto a la variable independiente t. una

ecuación diferencial ordinaria de primer orden con valores iniciales para la variable independiente puede representarse como:

96


.3': f(É,3)

y(f6l:

3,n

Se puede generalizar el problema anterior para el caso de que y sea un vector y = (y1 , Y2, ..., yn). Los comandos del módulo básico de MATLAB relativos a ecuaciones diferenciales ordinarias y algebraicas con valores iniciales se presentan en la tabla siguiente: Comando ode45 ode23

odel13 odel5s

Clase de problema que resuelve. método Ecuaciones diferenciales comunes por el Ecuaciones diferenciales comunes por el Ecuaciones diferenciales comunes por el

Ecuaciones diferenciales ordinarias

y

numérico v sintaxis método de Runse-Kutta método de Runse-Kutta método de Adams

algebraicas mediante NDFs

BDFs) ode23s

Ecuaciones diferenciales ordinarios por el método de Rosenbrock

ode23t

Ecuaciones diferenciales ordinarias

ode23tb

trapezoidal E cuaciones diferenciales ordinarias mediante TR-B DF 2

y

algebraicas

por la

regla

La sintaxis común para los 7 comandos anteriores es la siguiente:

Y] : solver (odefun, tspan, y0) IT,Y] : solver (odefun,tspan, y0, opciones) IT,Y] : sofver (odefun,tspan, y0ropciones,p1,p2. . .) ITrYrTE¡YE, IE] : solver (odefun,tspanry0ropciones) IT,

En la sintaxis anterior solver puede ser cualquiera de los comandos ode45, ode23, odel 13, ode15s, ode23s, ode23t u ode23tb.

El argumento odefun evalúa el lado derecho de la ecuación diferencial o sistema escrita en la forma y'=f(t,y) o M(t,y)y'=f(t,y), donde M(t,y) se denomina matriz de masa. El comando ode23s sólo puede resolver ecuac¡ones con matriz de masa constante. Los comandos ode15s y ode23t pueden resolver ecuaciones con matriz de masas singular y ecuaciones diferenciales algebraicas. El argumento tspan es un vector que especifica el intervalo [10, tf] de integración (para obtener soluciones en valores específicos de t, todos crecientes o decrecientes, se usa tspan = [10, t1, .,., tfl).El argumento y0 especifica un vector de condiciones iniciales. Los argumentos p1, p2,... son parámetros opcionales que se

CAPÍTULO 7. ECUACIONES DIFERENCIALES CON VALORES INICIALES... 97

pasan a odefun. El argumento opciones especifica opciones adicionales de integración mediante el comando odeset que pueden encontrarse en el manual del programa. Los vectores T e Y presentan los valores numéricos de las variables independiente y dependiente encontrados para las soluciones.

Como primer ejemplo hallamos las soluciones en el intervalo 10,121 del sistema de ecuaciones diferenciales ordinarias siguiente:

rr-

., ,^\

1 .?'o : -3,r J'r

3'*

:

fi

-vr(O)=0 J2(0) = I ó - t* -0"51srse ss(ül = 1 -,2-,

/

Para ello definimos una función de nombre sistemal en un Mfichero, con la finalidad de almacenar en ella las ecuaciones del sistema. La forma de definir esta función es partir de un vector columna con tres filas vacío, para asignarle posteriormente tres componentes que constituyen la sintaxis de las tres ecuaciones (FiguraT-1).

Figura

7-l

A continuación resolvemos el sistema tecleando en la ventana de comandos lo siguiente: >> [T,Y] = ode45(Gsistemal,lO L2),[0 1 1])

98


.0001 0.0001 0

0.0002 0.0002 0.0005

17.6736 71-7424 77 - 8772

12.0000

Y-

0313 0..0627

1.0000 1 .0000 1 .0000 1.0000 1 .0000 1 .0000 1.0000 1.0000 1.0000 7.0000 1.0000 1.0000 1 .0000 0.9999 0.9998 0.9997 0.9995 0.9980

1.0000 7.0000 1.0000 1.0000 7.0000 7.0000 1.0000 1.0000 1.0000 1.0000 1.0000 7-0000 1.0000 1.0000 0.9999 0.9998 0.9997

0.8594 0.7257 0.5228 0 .2695 -0 - 0778 -0 .29i6 -0 - 4098 -0.5769 -0.6735 -0.6987

-0.5705 -0.6876 -0.8524 -0.96i1 -0.9990 -0.9540 -0.9702 -0.85s9 -0.7874 -0-7128

0.7894 0.8552 0.9287 0.9815 0.9992

0

.0007 0 .0001

0

0

.0002

0-

0002

0.0005 0 .0007 0.0070 0 .0012 0 .002s 0.0037 0.0050 0.0062 0 .0725 0.0188 0-

0257

0-

0-

9990

0.976i 0-

9548

0.9279 0.8974 0.8650

Para interpretar mejor los resultados, la solución numérica anterior puede graficarse (Figura 7-2) mediante la siguiente sintaxis:

CAPiTULO 7. ECUACIONES DIFERENCIALES CON VALORES INICIALES... 99

FiguraT-2

7.3 ECUACIOI\ES DIFERENCIALES ORDINARIAS CON

VALORES EN LA FRONTERA MATLAB también perm¡te resolver ecuaciones diferenciales ordinarias con valores en la frontera. Los valores en la frontera no son más que condiciones que especifican una relación entre los valores de Ia solución en los puntos inicial y final de la variable independiente. El problema más sencillo de este tipo es el sistema de ecuaciones

.y'

:

f{s,y}

donde x es la variable independiente, y es la variable dependiente e

y' es la derivada de y respecto a x (y'=dy/dx). Además la solución en

el

intervalo [a, b] ha de satisfacer: sLY{n},}i&}J

:

ü

Una forma más general de este tipo de ecuaciones diferenciales puede expresarse como sigue:

1OO ECUACIONES DIFERENCIALES CON MATL,AB

.y"

:

f{:,:,3,,p}

stv(s],"v(á3,s)

:

ü

donde el vector p está formado por parámetros que han de ser determinados simultáneamente con la solución a través de las condiciones de valores en la frontera.

El comando que resuelve estos problemas es bvp4c, cuya sintaxis es la siguiente:

soI sol sol

bwp4c (odefun,

bcfun, solinit)

bvp4c (odefun, bcfun, bvp4c (odefun, bcfun,

solinlt, options) solinit, options ,p7,p2.

En la sintaxis anterior odefun es una función que evalúa las ecuaciones diferenciales f(x,y) y que puede tener la forma siguiente: dydx dydx dydx dydx

odefun odefun odefun odefun

(x, (x, (x, (x,

y) y,p1 ,p2, . . .) y. parametros

)

y.parametros,p7,p2,

..

.)

En la sintaxis de bvp4c el argumento bcfun es una función que computa el residuo en las condiciones en la frontera. Su forma es la siguiente:

res res res res

bcfun bcfun bcfun bcfun

(ya, (ya, (ya, (ya,

yb)

yb,pL,p2, . . .) yb, parametros ) yb, parametro s, p1-, p2,

El argumento solinit es una estructura con los campos x (nodos ordenados de la malla inicial de modo que las condiciones en la frontera son impuestas mediante a = solinit.x(1) y b=solinit.x(end) e y (suposiciones iniciales para la solución de modo que a = solinit.x(1) y b = solinit.x(end) es una suposición para la solución en el nodo solinit.x(i). El comando vpin¡t fija solinit con la sintaxis solinit = vpinit(x,y).

Como ejemplo resolvemos

la ecuación

orden:

y"+

L?l

=

ü

diferencial de segundo

CAPÍIULO 7. ECUACIONES DIFERENCIALES CON VALORES INICIALES,.. 101

cuyas soluciones han de satisfacer las condiciones en la frontera:

"y{t}} = ü *'{4} = -É El problema anterior es equivalente al siguiente:

3'r' = #E ss,

:

-lsrl

Consideramos como suposición inicial una malla de 5 puntos igualmente espaciados en [0,4] y valores constantes y1(x)=0 e y2(x)=Q. ¡¿ sintaxis siguiente recoge estas suposiciones:

>> so]-init = bwpinit(linspace(0,4,5), t1 0]); Codificaremos en MATLAB la función y las condiciones frontera mediante los M-ficheros de las Figuras 7-3 y 7-4.

en

la

FiguraT-4

FiguraT-3

La solución de la ecuación la ofrece la sintaxis siguiente:

)) sol = bwp4c(Gt¡roode,Gtwobc,solinit)

,'

que puede graficarse (Figura 7-5) usando el comando bvpval como sigue:

))y = brpval (so1 ,linspace >>plot(x,y(1 ,:)) ;

(0

,4) I ;

702


r§

Í

"*

Figura 7-5 7

.4 E,CU ACIONES DIFERENCIALES EN DERIVADAS

PARCIALES El módulo básico de MATLAB dispone de funciones que permiten resolver ecuaciones y sistemas de ecuaciones diferenciales en derivadas parciales con valores iniciales en la frontera. La función básica para el cálculo de las soluciones es pedepe, y la función básica para evaluar dichas soluciones es pdeval. La sintaxis de la función pedepe es la siguiente:

sol : sol : sol : pdepe

pdepe (m, pdefun, icfun, bcfun, xmesh, tspan) pdepe (m,pdefun, icfun,bcfun,xmesh, tspan, options) (m,

pdefun, icfun, bcfun,

>«nesh,

tspan, options, p1, p2 . . . )

El parámetro m vale 0, 1 o 2 según la naturaleza de la simetría (bloque, cilíndrica o esférica, respectivamente). El argumento pdefun es la función que define las componentes de la ecuación diferencial, icfun es la función que define las condiciones iniciales, bcfun es la función que define las condiciones frontera, xmesh y tspan son los vectores [x0, x1,...,xn] y [t0, t1,...,tf] que especifican los puntos en los cuales se requiere la solución (n,f>3), options especifica las opciones de cálculo de las soluciones (RelTol, AbsTol, NormControl, lnitialStep, y MaxStep para especificar tolerancia relativa, tolerancia absoluta, tolerancia en norma, paso inicial y paso máximo, respectivamente) y p1,p2,... son parámetros a pasar a las funciones pdefun, icfun y bcfun. La forma de la ecuación diferencial en derivadas parciales es:

CAPíTULO 7. ECUACIONES DIFERENCIALES CON VALORES INICIALES... 1.03

*] ñuár

\ r, ","dr/ "f", donde

alxlb y

áu\\ f áu\ x:r*a -á ¡l.TmÍr{.r, r, H, -., JJ + 8f f,, f, il, 1- j , \ d:\ cl§lr óxr t0
t = t0 y para todo x las

componentes solución satisfacen las condiciones iniciales:

u{x,

t6} :

tsn{s}

y para todo t y cada x=a o x=b, Ias componentes solución satisfacen ecuaciones frontera de la forma: F{n, f, r*} + g{r,

f}d'*,r, L

u,

}J dx!

=

O

Además, se tiene que xmesh(1)=a, xmesh(end)=b, tspan(1)= t0 y tspan(end)= tf. Por otra parte pdfun halla los términos c, f y s de la ecuación diferencial en derivadas parciales, de modo que:

Ic, frs]

: pdefun(x,tru,dudx)

De modo similiar icfun evalúa las condiciones iniciales mediante

u : icfun

(x)

Por último, bcfun evalúa los valores b y c de las condiciones en la frontera mediante:

Ip1,ql,pr,er] : bcfun (x1,u1,xr,ur,t) Como primer ejemplo resolveremos parciales (xe[0,1] y t>0) siguiente:

9Eu

,L.-=

á /Eu\ _ ás -t \drl I

dt

satisfaciendo la condición inicial ax(r,

:

fi! : sjñfin

y las condiciones en la frontera

:

la

ecuación en derivadas

104

ECUACIONES DIFERENCIALFJ CON MATLAB

r*[0, *] =

Ü

du-

-+-11,f):0 tt€-+ ds'

Comenzamos definiendo las funciones en M-ficheros de las figura 7-6 a7-8.

FiguraT-7

Figura 7-6

Figura 7-8

Una vez definidas las funciones de apoyo, se define la función que resuelve el problema mediante el M-fichero de la FiguraT-9.

CAPíTULO 7. ECUACIONES DIFERENCIALES CON VALORES INICIALES... 105

FigraT-9 Para ejecutar la solución (Figuras 7-10 comandos de MATLAB se utiliza la sintaxis:

y T-11), en la ventana

)) pdexl

Figura 7-10

Figura 7-11

de

106


Como segundo ejemplo resolvemos

el sistema de ecuaciones

diferenciales en derivadas parciales (xe[0,1] y t>0) siguiente:

--2 du,

+o[ :

¡.pa

da

i]

__ü

d¡'

-F6r-u2\

dt¿n e)-tt *_'.f; "# +Fqur-a2) = ü.t7o ox.

F(v

) = expi5.?§3] - exp{-11.4$3}

satisfaciendo las condiciones iniciales:

ur(r,0)=

1

ug(r,0)=

0

y las condiciones en la frontera: dra.,

dx =+(0,

f) = 0

f} = 0 ur(1, f) = I rcg{f},

r): o P,r, (,.tr Para utilizar convenientemente la función pdepe, el sistema puede escribirse como: r'1

lil _ ul,,J_-

Lrl

,rl"rl

afo.m+¡a,,rraxi] .[-rru,-u2r'l aTlo.rrotaurrar{ [^r't*, -,,rr]

La condición izquierda en la frontera puede escribirse como:

Irl - lr] * [o.m+{a*,rr"il LI¿rJ

=

LrJ fo.lzotaul,rartl

lo] Lrl

y la condición derecha en la frontera puede escribirse como:

CAPITULO 7. ECUACIONES DIFERENCIALES CON VALORES INICIALES... 107

h, -

L

,'l lo]

fo"m+1aur.a.fl= fol

, J- LrJ'. lo.rzotar,rra,{

t-l

Comenzamos definiendo las funciones en M-ficheros de las Figuras 8-39 a 8-41.

FiguraT-12

Figura 7-13

FigwaT-14 Una vez definidas las funciones de apoyo, se define la función que resuelve el problema mediante el M-fichero de la Figura 7-15.

108


Figura 7-15

Para ejecutar la solución (Figuras 7-16 comandos de MATLAB se utiliza la sintaxis:

y 7-17), en la ventana

>> pdex4

FigluraT-16

FigxaT-17

de

CAPÍTULO 7. ECUACIONES DIFERENCIALES CON VALORES INICIALES... 109

(

Comenzamos definiendo una función de nombre vdp100 en un Mfichero, con la finalidad de almacenar en ella las ecuaciones del sistema. La forma de definir esta función es partir de un vector columna con dos filas vacío, para asignarle posteriormente tres componentes que constituyen la sintaxis de las tres ecuaciones (FiguraT-18).

Figura 7-18

A continuación resolvemos el sistema y graficamos la solucióñ !1 = y1(t) relativa a la primera columna (Figura 7-19) de la salida tecleando lo siguiente en la ventana de comandos:

>> [T,Y] = ode15s(0v$1000, [0 >>

plot(T,Y(:,1),'-')

3000]

,12

oll;

110


Figura 7-19

FiguraT-20

De forma similar graficamos la solución y2 = yZ(t) (Figura 7-20) mediante la sintaxis:

)} plot(T,Yl:,2),'-')

La ecuación dada es equivalente al sistema de diferenciales de primer orden siguiente:

-?1 =

-vZ

3g' = -{A -g

q ffiE g"n}3t

con las siguientes condiciones en la frontera:

3r{ü}-1 = 0 rttÜl : ü }g(a) = 0

ecuaciones

CAPíTULO 7, ECUACIONES DIFERENCIALES CON VALORES INICIALES... 111

Para representar el sistema usamos la función del M-ficheros de la Figura 7-21, para representar las condiciones en Ia frontera usamos la función del M-fichero de la Figura z-22, y mediante la función del M-fichero de la FiguraT-23 realizamos suposiciones para la solución inicial.

FiguraT-21

FigraT-22

FigtraT-23 La solución inicial para l"=1s y 10 puntos igualmente espaciados en n] se calcula mediante la sintaxis MATLAB siguiente: [0,

>> Ia¡nbda = 15;

so].init = bwpinit(Iinspace (0,pi,LOl,Gmat4init,larnbda),. La soluciÓn numérica del sistema se calcula mediante la sintaxis

siguiente:

>> soI = bvp4c(G¡nat4ode,Gmat4bc,solinit),. Para graficar su primera componente en 100 puntos igualmente espaciados en el intervalo [0, n] se utiliza la siguiente siniaxis:

)) xint = linspace(0,pi); Sxint = bwpval_(sol ,xint),. plot(xint,sxint(1, : ) )

T12


axis ( [0 pi -1 1.1] xJ-abe1 ( 'x'

)

)

ylabel('solucion y') El resultado se presenta en la FiguraT-24.

FigaraT-24

La ecuación general anterior es

equivalente

al

sistema de

ecuaciones lineales de primer orden siguiente: n,l y 1

:

3'E

= P{1-3r-be--Yr

Y6 E)

cuyas ecuaciones pueden definirse para p=1 mediante el M-fichero de la Figura 7-25.

capirulo

7. EcuActoNES DtFERENcTALES coN vALoRES tNtc¡ALES... 113

FigxaT-25 Si consideramos valores iniciales y1=2 e y2=0 en el intervalo [0,20], se puede resolver el sistema mediante la sintaxis MATLAB siguiente:

>>

[t,y] = ode45(Gvdpl, [0 20] ,f2;

L-

0

0.0000 0.0001 0.0001 0.0001 0.0002 0.0004 0.0005 0.0006 0.0.012

19.9559 19.9780 20.0000 de 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000

7.8729 1.9358

0

-0.0001 -0.0001 -0.0002 -0.0002 -0.000s 1.0366 0.7357

071

TLA


7.9787 2.0046 2.0096 2.0133 2.0158 2.0172

0.4746 0.2562 0.1969 0.7473 0.0892 0.0404

Podemos graficar las soluciones (Figura 7-26) mediante la sintaxis.

plot(t,y( i ,Ll ,'-' ,t,y (: ,2) ,'--') ) xlabel('tiempo t') >>

>> yJ-abe1 ('solucion y' ) >} legend('y_1 ' ,'y_2')

lt

,tt

iI

i.

FiguraT-26

Para resolver el sistema general con el parámetro p se define el sistema mediante el M-fichero de la FiguraT-27.

fuic,§,iüB:

:drdt ='rdpa (qyrmul

dydc = ty{2)

;

uu+

{I-y{I)

"2} *y(2)

:*

-y(l¡ l,

r,t'

::=}I FigwaT-27

r* l

CAPíTULO 7. ECUACIONES DIFERENCIALES CON VALORES INICIALES... 115

Ahora podemos graficar la primera solución (Figura 7-28) el intervalo [0,1500] para y1=2 e y2=O

correspondiente a p=1000 en mediante la sintaxis siguiente:

>> [t,y] = ode15s(0vdp2, [0 1500], [2; 0], [],1000); )) xIabel('tiempo t') >> ylabeI ('solucion y_l')

Si queremos graficar la primera solución para otro valor del parámetro, por ejemplo p=100, en el intervalo [0,1500] para y1=2 e y2=0 (Figura 7-29), utilizaremos la sintaxis:

>> >>

[t,y] = odelSs(Gvdp2, [0 1500], [2; 0], il,100);

plot(t,y( | ,Ll ,'-')

FigwaT-28

;

FiguraT-29

tsBN 9781491271537

llffi[illllil]il[ililll

lilflfiflfltl

Ecuaciones Diferenciales Con MATLAB Cesar Perez Lopez

Recommend Documents