F. Palacios 2012

Computación científica y técnica con Maple

Algunos ejemplos ilustrativos F. Palacios Quiñonero Dept. Matemática Aplicada III. Universidad Politécnica de Cataluña Ver. 1.1

Abril, 2012

Este documento se ha editado usando las hojas de trabajo (worksheet) de Maple 15

Contenido 1. Transformada de Laplace 2. Resonancia 3. Resolución de EDOs con dsolve( ) 4. Sistema de primer orden con excitación sinusoidal 5. Resolución numérica 6. EDO de segundo orden con amortiguamiento 7. Dígitos, programación

1 Transformada de Laplace > with inttrans addtable, fourier, fouriercos, fouriersin, hankel, hilbert, invfourier, invhilbert, invlaplace, invmellin, laplace, mellin, savetable Transformada de Laplace > f d cos w t f := cos w t > F d laplace f, t, s s F := 2 2 s Cw

(1.1)

(1.2) (1.3)

Transformada inversa > invlaplace F, s, t cos w t

(1.4)

Función definida por tramos > fd

2

t !1

K1

t !3

1

otherwise

f :=

2

t !1

K1

t !3

1

otherwise

(1.5)

Representación gráfica > plot f, t =K2 ..4, color = blue, thickness = 3, gridlines 1

F. Palacios 2012

2 1

K2

K1

0

1

2 t

3

4

K1 > F d laplace f, t, s Ks K3 s 2 K3 e C2 e F := s

(1.6)

> g d invlaplace F, s, t g := 2 K3 Heaviside K1 Ct C2 Heaviside K3 Ct > f1 d convert g, piecewise 2 t !1 f1 :=

undefined

t=1

K1

t !3

undefined

t=3

1

3 !t

(1.7)

(1.8)

> plot f1, t = 1 ..4

1 0.5 0

2

3

4

t K0.5 K1 Resolución de EDOs lineales usando Laplace 2 .. > edo d x Cω0 x = sin w t edo :=

2 d2 x t Cω0 x t = sin w t 2 dt

(1.9)

> edoT d laplace edo, t, s 2

edoT := s2 laplace x t , t, s KD x 0 Ks x 0 Cω0 laplace x t , t, s =

w s Cw2 2

(1.10)

> alias X = laplace x t , t, s 2

F. Palacios 2012

X

(1.11)

> edoT w s Cw2

2

s2 X KD x 0 Ks x 0 Cω0 X =

(1.12)

2

Condiciones iniciales en notación matemática usual . > cis d x 0 = v0, x 0 = x0 cis := x 0 = x0, D x 0 = v0

(1.13)

> edoT2 d eval edoT, cis w s Cw2

2

edoT2 := s2 X Kv0 Ks x0 Cω0 X =

(1.14)

2

Resolvemos respecto de X > solX d solve edoT2, X solX :=

w Cv0 s2 Cv0 w2 Cs3 x0 Cs x0 w2 2

2

s Cw

2

s

(1.15)

2 Cω0

Anti-transformamos > ft d invlaplace solX, s, t ft := x0 cos ω0 t C

sin w t 2 ω0 Kw2

C

sin ω0 t

2

w Kω0 v0 Cv0 w2 2

ω0 w

(1.16)

2 Kω0

2 Resonancia > fr d limit ft, w = ω0 2

1 2 x0 cos ω0 t ω0 C2 sin ω0 t v0 ω0 Csin ω0 t Kcos ω0 t t ω0 fr := 2 2 ω0 > fr1 d collect fr, t, cos, sin 1 cos ω0 t t 1 fr1 := K Cx0 cos ω0 t C 2 2 ω0

2 ω0 v0 C1 sin ω0 t 2

ω0

Gráficos interactivos > plots :-interactiveparams 'plot', (2.2), t = 0 .. 25, labels = t, "" , = 0 .. 10, x 0 = 0 .. 10

(2.1)

(2.2)

, ω 0 = 0 .. 30, v 0

3

F. Palacios 2012

0.8 0.6 0.4 0.2 0 5

K0.2

10

15

20

25

t

K0.4 K0.6 K0.8

3 Resolución de EDOs con dsolve( ) > edo 2 d2 x t Cω0 x t = sin w t 2 dt

(3.1)

> f

> edo2 d

2

t !1

K1

t !3

1

otherwise

(3.2)

2 d2 x t Cω x t =f 0 dt2 2

edo2 := > sol d dsolve

d x t dt2

2 Cω0 x

t =

2

t !1

K1

t !3

1

otherwise

(3.3)

. edo2, x 0 = v0, x 0 = x0 sin ω0 t v0 ω0

C

sin ω0 t v0

sol := x t =

ω0 sin ω0 t v0 ω0

C

cos ω0 t

C

cos ω0 t

cos ω0 t

2 ω0 2

x0 ω0 K2 2

ω0 2

x0 ω0 K2 2

ω0

2

x0 ω0 K2

C

1 2

ω0

C

K

1 2

ω0

C

C

2

2

ω0

3 cos ω0 K1 Ct

3 cos ω0 K1 Ct ω0

2

2

ω0 K

2 cos ω0 K3 C 2

ω0

4

F. Palacios 2012

4 Sistema de primer orden con excitación sinusuoidal . > edo d y Ca y = sin w t d y t Ca y t = sin w t dt

edo :=

(4.1)

> sol d dsolve edo, y 0 = y0 sol := y t = eKa t y0 C

w a Cw2 2

C

Kcos w t w Csin w t a a2 Cw2

(4.2)

> fy d rhs sol fy := eKa t y0 C

w a Cw2 2

C

Kcos w t w Csin w t a a2 Cw2

(4.3)

> vp d a = 5, w = 360, y0 = 1 vp := a = 5, w = 360, y0 = 1

(4.4)

> fp d eval fy, vp fp :=

25997 K5 t 72 1 e K cos 360 t C sin 360 t 25925 25925 25925

(4.5)

> fp1 d evalf fp fp1 := 1.002777242 eK5. t K0.002777242044 cos 360. t C0.00003857280617 sin 360. t > plot fp1, t = 0 ..3

(4.6)

1 0.8 0.6 0.4 0.2 0 0

1

2

3

t Observamos que no hay respuesta oscilante (ver fórmula)

5 Resolución Numérica . > edo d y C

t y = sin t 1 Ct2 edo :=

d ty t y t C = sin t dt 1 Ct2

(5.1)

> s d dsolve edo, y 0 = 0 t

sin _z1 s := y t =

1 C_z12 d_z1

0

(5.2) 2

1 Ct dsolve( ) no obtiene solución exacta; pasamos a resolución numérica > sn d dsolve edo, y 0 = 0 , numeric (5 53)

F. Palacios 2012

sn := proc x_rkf45 ... end proc

(5.3)

t = 1.5, y t = 0.728338957113725

(5.4)

> sn 1.5 > with plots animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fieldplot, fieldplot3d, gradplot, gradplot3d, implicitplot, implicitplot3d, inequal, interactive, interactiveparams, intersectplot, listcontplot, listcontplot3d, listdensityplot, listplot, listplot3d, loglogplot, logplot, matrixplot, multiple, odeplot, pareto, plotcompare, pointplot, pointplot3d, polarplot, polygonplot, polygonplot3d, polyhedra_supported, polyhedraplot, rootlocus, semilogplot, setcolors, setoptions, setoptions3d, spacecurve, sparsematrixplot, surfdata, textplot, textplot3d, tubeplot > odeplot sn, t = 0 ..5

(5.5)

1.0 y 0.6 0

1

2

3

4

5

t Una segunda EDO . 1 > edo2 d y C y = sin t 1 Ct2 d y t y t C = sin t dt 1 Ct2 > sn2 d dsolve edo2, y 0 = 0 , numeric sn2 := proc x_rkf45 ... end proc edo2 :=

(5.6)

(5.7)

> sn2 3 t = 3., y t = 1.49438163398571

(5.8)

> odeplot sn2, t = 0 ..5

1.4 1 y 0.6 0

0

1

2

3

4

5

t >

6 EDO de segundo orden con amortiguamiento .. . 2 > edo d x C2 ζ ω0x Cω0 x = sin w t 6

F. Palacios 2012

d2 x t C2 ζ ω0 dt2 . > sol d dsolve edo, x 0 = 0, x 0 = 0 sol := x t

d x t dt

edo :=

Kζ C ζ2 K 1 2

2 4

2

ω t 0

2

2

2

ζ K1 w

3

2

w2 ζ Kw2 C2 ω0 ζ K3 ζ ω0 Cω0 K2 ζ K ζ C ζ2 K 1

1 K 2

ω t 0

e 2

2

2

(6.1)

(6.2)

e

1 = 2

2

Cω0 x t = sin w t

2

2

2

2

ζ K1 ω0 C2 ζ ω0

ζ K1

ω0

w 2

2

ζ K1 Kω0 Cw2

ζ K1 ω0 2 ζ ω0 C2 ζ ω0 2

ω0 Kw2 sin w t K2 cos w t w ζ ω0

C

4

2

2

ω0 C 4 ζ K2 w2 ω0 Cw4

Caso subamortiguado . > sol1 d dsolve edo, x 0 = 0, x 0 = 0 Kω t ζ 0

e

sol1 := x t = K Kω t ζ 0

2e

C

2

Kζ ω0 Cω0 w 2 ζ ω0 Kω0 Cw2

2

2

4

2

2

2

2

2

Kζ ω0 Cω0 t

sin

2

2

2

2

2

4

2

4

ω0 4 w2 ζ ω0 K6 ζ ω0 w2 Cζ w4 Kw4 C2 ω0 w2 Cζ ω0 Kω0 2

2

2

Kζ ω0 Cω0 t ζ ω0 w

cos 2

2

2

2

assuming 0 ! ζ, ζ ! 1, ω0 O 0

4

2

4 ζ ω0 w2 Cω0 K2 ω0 w2 Cw4

(6.3)

2

ω0 Kw2 sin w t K2 cos w t w ζ ω0

C

4

2

2

ω0 C 4 ζ K2 w2 ω0 Cw4

> sol2 d collect sol1, exp, sin, cos sol2 := x t = K

2

2

2

2

2

4

2

2

2

2

2

2

2

2

Kζ ω0 Cω0 w 2 ζ ω0 Kω0 Cw2

Kζ ω0 Cω0 t

sin

2

2

4

2

4

ω0 4 w2 ζ ω0 K6 ζ ω0 w2 Cζ w4 Kw4 C2 ω0 w2 Cζ ω0 Kω0 2

C

2

2

Kζ ω0 Cω0 t ζ ω0 w

2 cos 2

2

4

Kω t ζ 0

e

2

2 2 4 4 ζ ω0 w Cω0 K2 ω0 w Cw

(6.4)

2

C

ω0 Kw2 sin w t 4

2

2

ω0 C 4 ζ K2 w2 ω0 Cw4

2 cos w t w ζ ω0

K

4

2

2

ω0 C 4 ζ K2 w2 ω0 Cw4

Podemos extraer subexpresiones seleccionando y arrastrando 2

> K

2

2

2

2

2

2

2

2

2

2

2

ω0 Kζ ω0 w 2 ζ ω0 Kω0 Cw2

ω0 Kζ ω0 t

sin

2

4

2

2

4

2

4

ω0 Kw4 Cζ w4 K6 ζ ω0 w2 C4 w2 ζ ω0 C2 ω0 w2 Kω0 Cζ ω0 K

sin 2 ω0

2

2

2

Kζ ω0 Cω0 t 2 4

4w ζ

2

2

2

2

2

2

Kζ ω0 Cω0 w 2 ζ ω0 Kω0 Cw2

2 2 2 2 ω0 K6 ζ ω0 w2 Cζ w4 Kw4 C2

2 4 2 4 ω0 w2 Cζ ω0 Kω0

(6.5)

7

F. Palacios 2012

> CodeGeneration 'Matlab' '(6.5)' cg = -sin(sqrt(-Zeta ^ 2 * omega(0) ^ 2 + omega(0) ^ 2) * t) * sqrt(-Zeta ^ 2 * omega(0) ^ 2 + omega(0) ^ 2) / omega(0) ^ 2 * w * (0.2e1 * Zeta ^ 2 * omega(0) ^ 2 - omega(0) ^ 2 + w ^ 2) / (0.4e1 * w ^ 2 * Zeta ^ 4 * omega(0) ^ 2 - 0.6e1 * Zeta ^ 2 * omega(0) ^ 2 * w ^ 2 + Zeta ^ 2 * w ^ 4 - w ^ 4 + 0.2e1 * omega (0) ^ 2 * w ^ 2 + Zeta ^ 2 * omega(0) ^ 4 - omega(0) ^ 4); Podemos obtener la expresión Matlab usando el Menú de Contexto

7 Dígitos, programación 7.1 Resolución ecuaciones > restart Definimos las expresiones y realizamos un gráfico > f d exp Kx f := eKx

(7.1.1)

> g d x2 > plot

g := x2 f, g , x =K2 ..2, color = red, blue

(7.1.2)

7 6 5 4 3 2 1 K2

K1

0

1 x

2

Solución exacta con solve( ) > s d solve f = g 1 2

s := 2 LambertW K

, 2 LambertW

1 2

(7.1.3)

> sf d evalf s sf := K1.588047265 C1.540223501 I, 0.7034674224

(7.1.4)

s1 := 0.7034674224

(7.1.5)

0.7034674224

(7.1.6)

s2 := 0.7034674225

(7.1.7)

> s1 d sf 2 > sf2 Solución numérica con fsolve( ) > s2 d fsolve f = g >

8

F. Palacios 2012

7.2 Método de Newton-Raphson Método de Newton-Raphson para ecuaciones f x = 0 x0 aproximación inicial xj C 1 = xj K

f xj f' xj

Programa del método y ejecución con 3 iteraciones > f d x/eKx Kx2; x0 d 0.7; df d D f ; n d 3; for j from 0 to n K1 do f xj xj C 1 d xj K df xj end do; f := x/eKx Kx2 x0 := 0.7 df := x/KeKx K2 x n := 3 x1 := 0.7034721896 x2 := 0.7034674225 x3 := 0.7034674225

(7.2.1)

eKa Ka2

(7.2.2)

KeKa K2 a

(7.2.3)

0.

(7.2.4)

> f a > df a > f x3

9

F. Palacios 2012

El método de Newton-Raphson es de orden 2, ej C 1 z k ej 2, donde ej = α Kxj, y α es la solución de la ecuación. Esto intuitivamente significa que podemos doblar el número de decimales exactos a cada iteración. > Digits d 50; f d x/exp Kx Kx2; x0 d 0.7; df d D f ; n d 6; for j from 0 to n K1 do f xj xj C 1 d xj K df xj end do; Digits := 50 f := x/eKx Kx2 x0 := 0.7 df := x/KeKx K2 x n := 6 x1 := 0.70347218961269235924710536764659532536847293171448 x2 := 0.70346742250738434708832541148616093662142459848785 x3 := 0.70346742249839165204985060257268964960601811365475 x4 := 0.70346742249839165204981860185990213034292843143945 x5 := 0.70346742249839165204981860185990213034292843103422 x6 := 0.70346742249839165204981860185990213034292843103423

(7.2.5)

K1. 10-50

(7.2.6)

> f x6

10