12 - EQUACOES-DIFERENCIAIS.sagews

%auto
typeset_mode(True, display=False)

EQUAÇÕES DIFERENCIAIS

Equações Diferenciais de Primeira Ordem

Equações Diferenciais de Primeira Ordem Simbólicas

O comando desolve calcula a "solução geral" a uma equação de 1ª ou 2ª ordem utilizando o Maxima. Sintaxe:

desolve(EDO, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False)

EDO - expressão da equação diferencial

dvar - variável dependente

ics - (optional) condições iniciais de contorno

Para uma equação de primeira ordem, especifique as condições iniciais de [x0 e y(x0)].
Para uma equação de segunda ordem,especifique as condições iniciais de x0 , [y(x0) e y'(x0)]
Para uma solução de contorno de segunda ordem , especifique as condições de x e de contorno de y na forma [x0,y(x0),x1,y(x1)]

ivar - (optional) variável independente

show_method - se verdadeiro, Sage retorna o par [solução, método], onde o method é o método que tem sido usado para obter uma solução (Maxima usa a seguinte ordem de equações de primeira ordem: linear, separável, exata (incluindo exata com fator de integração), homogênea , Bernoulli, generalizada homogênea).

contrib_ode - (optional) Se verdadeiro, desolve permite resolver Clairaut, Lagrange, Riccati e algumas outras equações. Isso pode levar um longo tempo e é assim desativado por padrão. Ex: Resolver a seguinte equação diferencial

$y' - 2xy = 0$

reset()

y = function('y')(x)
y

y(x)

edo1 = y.derivative(x,1) - 2*x*y
edo1

$-2 \, x y\left(x\right) + D[0]\left(y\right)\left(x\right)$

desolve(edo1,y)

$C e^{\left(x^{2}\right)}$

desolve( y.derivative(x,1) - 2*x*y , y)

$C e^{\left(x^{2}\right)}$

Utilizando show_method = true. Sage retorna o par [solução, método], onde o method é o método que tem sido usado para obter uma solução (Maxima usa a seguinte ordem de equações de primeira ordem: linear, separável, exata (incluindo exata com fator de integração), homogênea , Bernoulli, generalizada homogênea).

desolve(edo1,y,show_method = true)

[$C e^{\left(x^{2}\right)}$, linear]

desolve( y.derivative(x,1) - 2*x*y  , y, show_method = True)

[$C e^{\left(x^{2}\right)}$, linear]

Equação diferencial com valores iniciais

Ex: Resolver a seguinte equação diferencial no intervalo [0,3], sendo y(0) = 1

$y' - x^2 + y = 0$

reset()

y = function('y')(x)

edo2 = y.derivative(x,1) - 2*x*y
edo2

$-2 \, x y\left(x\right) + D[0]\left(y\right)\left(x\right)$

$-2 \, x y\left(x\right) + D[0]\left(y\right)\left(x\right)$

sol2 = desolve(edo2, y, ics = [0, 1], show_method = true)
sol2

[$e^{\left(x^{2}\right)}$, linear]

Mais um exemplo

$y'(x) = 2*x*y$

y = function('y')(x)

sol3 = desolve ( y.derivative(x,1) == 2*x*y , y , ics =[0 , 1], show_method=true)
sol3

(sol3[0](1)).n()

Error in lines 1-1 Traceback (most recent call last): File "/usr/local/lib/python2.7/dist-packages/smc_sagews/sage_server.py", line 957, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> NameError: name 'sol3' is not defined

Gráfico da solução da equação diferencial de primeira ordem

type(sol3)
sol3[0]

<type 'list'>

$e^{\left(x^{2}\right)}$

plot(sol3[0],(x,0,1), figsize = (4,3), gridlines="minor")

Campo de Direções

Ex: Resolver a seguinte equação diferencial no intervalo [0,1], sendo y(0) = 1

$y'(x)=2xy$

Obter o campo de direções vetoriais

var('y')

plot_campo = plot_slope_field(2*x*y, (x,0,1), (y,0,5), figsize = (4,3))
plot_campo

$y$

Resolver a equação diferencial $y'(x)=f(x,y)$

y = function('y')(x)

edoC1 = y.derivative(x,1) - 2*x*y  == 0

solC1 = desolve(edoC1,y)
print(solC1)

_C*e^(x^2)

Plotar o gráfico com os valores de C0 e C1

# substituir o valor da constante por 1 e 2

solP1 = solC1.subs(_C = 1)
solP2 = solC1.subs(_C = 2)

plot_campo + plot(solP1, (x,0,1)) + plot(solP2, (x,0,1), color = 'red', figsize = (4,3))

Campo com solução exata

# Resolver a equação diferencial

y = function('y')(x)

edoC2 = y.derivative(x,1) - 2*x*y  == 0

solC2 = desolve(edoC2, y , ics = [0,1])
print(solC2)

e^(x^2)

# Testando a solução
solC2(0.5)

$1.28402541668774$

# Plotar a solução com o campo gerado anteriormente

plot_campo + plot(solC2, (x,0,1), color = 'red', figsize = (4,3))

Equação diferencial de primeira ordem numérica

Método de Euler

O método de Euler é empregado para encontrar numericamente a solução da EDO na forma $y ' = f ( x , y )$ , com valor inicial $y(xo) = k$. Os resultados são apresentados em uma tabela, no qual a primeira coluna temos os valores de x incrementados por h e na segunda coluna temos o novo valor de y. Sintaxe

eulers_method(f, x0, y0, h, x1, algorithm='table')

f: função equação diferencial

x0: condição inicial de x

y0: condição inicial de y(x0)

h: incremento do vetor x.

x1: extreo final do vetor x

Ex: Resolver a seguinte equação diferencial no intervalo [0,3], sendo y(0) = 1

$y' = 2*x*y$

Solução exata para y(0.5) = 1.28402541668774

reset()

var ( 'x , y ')

# intervalo
a = 0
b = 0.5

# condições iniciais y(x0)
y0 = 1

#passo de h
h = 0.05

# função
fDif  = 2*x*y

($x$, $y$)

eulers_method( fDif , a , y0 , h, b )

x y h*f(x,y) 0 1 0.000000000000000 0.0500000000000000 1.00000000000000 0.00500000000000000 0.100000000000000 1.00500000000000 0.0100500000000000 0.150000000000000 1.01505000000000 0.0152257500000000 0.200000000000000 1.03027575000000 0.0206055150000000 0.250000000000000 1.05088126500000 0.0262720316250000 0.300000000000000 1.07715329662500 0.0323145988987500 0.350000000000000 1.10946789552375 0.0388313763433312 0.400000000000000 1.14829927186708 0.0459319708746832 0.450000000000000 1.19423124274176 0.0537404059233794 0.500000000000000 1.24797164866514 0.0623985824332572

# menor h
h = 0.02

eulers_method(fDif, a , y0 , h , b)

x y h*f(x,y) 0 1 0.000000000000000 0.0200000000000000 1.00000000000000 0.000800000000000000 0.0400000000000000 1.00080000000000 0.00160128000000000 0.0600000000000000 1.00240128000000 0.00240576307200000 0.0800000000000000 1.00480704307200 0.00321538253783040 0.100000000000000 1.00802242560983 0.00403208970243932 0.120000000000000 1.01205451531227 0.00485786167349890 0.140000000000000 1.01691237698577 0.00569470931112030 0.160000000000000 1.02260708629689 0.00654468535230009 0.180000000000000 1.02915177164919 0.00740989275587416 0.200000000000000 1.03656166440506 0.00829249331524051 0.220000000000000 1.04485415772030 0.00919471658793867 0.240000000000000 1.05404887430824 0.0101188691933591 0.260000000000000 1.06416774350160 0.0110673445324167 0.280000000000000 1.07523508803402 0.0120426329859810 0.300000000000000 1.08727772102000 0.0130473326522400 0.320000000000000 1.10032505367224 0.0140841606870047 0.340000000000000 1.11440921435924 0.0151559653152857 0.360000000000000 1.12956517967453 0.0162657385873132 0.380000000000000 1.14583091826184 0.0174166299575800 0.400000000000000 1.16324754821942 0.0186119607715108 0.420000000000000 1.18185950899093 0.0198552397510477 0.440000000000000 1.20171474874198 0.0211501795778589 0.460000000000000 1.22286492831984 0.0225007146810851 0.480000000000000 1.24536564300093 0.0239110203456178 0.500000000000000 1.26927666334654 0.0253855332669309

Método de Euler Interativo

reset()
var('x y')

@interact
def Euler_interativo(df_dx = input_box(default= '2*x*y'),y0 =input_box(default= '1'),   h = slider(0.01, 0.1) ):

    a = 0.0
    b = 0.5
    points = eulers_method(df_dx,a,y0,h,b, algorithm = 'none')

    point2d(points,gridlines="minor", figsize = (4,3)).show()

($x$, $y$)

Interact: please open in CoCalc

Método de Rage-Kutta

O resultado de 'dsolve_rk4' é uma lista contendo os resultados de x, y, y' e assim por diante.Sintaxe

desolve_rk4(de, dvar, ics=None, ivar=None, end_points=None, step=0.1, output='list')

de : Equação diferencial a ser resolvida na forma y' = f (x, y)

dvar : variável dependente (variável simbólica definida)

ics: Condições iniciais na forma [x0 , y(x0)]

ivar: variável independente (declarada como função da variável independente)

end_point Intervalo. Se vazio, é usando end_points=ics[0]+10

step : (Opcional) tamanho do passo de integração, padrão é 0.1

output: (Optional) saída do resultado. Pode ser: 'list', 'plot' ou 'slope_field'. O padrão é 'list'

reset()

Solução Numérica $y' - 2*x*y == 0$. Colocando na forma correta para o rk4 $~fvar = 2*x*y$.

Solução exata para y(0.5) = 1.28402541668774

var('x,y')

edo_rk4 = 2*x*y
h = 0.1

sol_rk4 =desolve_rk4( edo_rk4 , y , ics = [0,1] , ivar = x , end_points=[0 , 0.5], step = h)
sol_rk4

($x$, $y$)

[[$0$, $1$], [$0.1$, $1.01005016667$], [$0.2$, $1.04081076977$], [$0.3$, $1.09417426548$], [$0.4$, $1.1735108136$], [$0.5$, $1.28402525566$]]

Gráfico da solução

desolve_rk4(edo_rk4 , y , ics = [0,1] , ivar = x , end_points = [0,1], step = h, gridlines="minor", output='plot', figsize = (4,3))

Gráfico com campo de direções

desolve_rk4(edo_rk4 , y , ics=[0,1] , ivar=x , end_points=[0,1], step = h, output='slope_field', figsize = (4,3))

Apesar do manual confirmar que a forma da equação diferencial de primeira ordem deve ser na forma " fvar = f(x,y)", podemos colocar de forma implícita que na grande maioria dos casos será resolvida.

y = function('y')(x)

edoRK = y.derivative(x,1) - 2*x*y  == 0

desolve_rk4(edoRK, y , ics=[0,1] , ivar=x , end_points=[0,0.5], step = h)

Error in lines 3-3 Traceback (most recent call last): File "/usr/local/lib/python2.7/dist-packages/smc_sagews/sage_server.py", line 957, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> NameError: name 'h' is not defined

Equações Diferenciais de Segunda Ordem

Equações Diferenciais de Primeira Ordem Simbólicas

reset()

Solução por dsolve()

Resolver a equação diferencial: $y''(x)- 4y'(x) - x = 0$

y = function('y')(x)
y

$y\left(x\right)$

edo_2ord = y.derivative(x,2) - 4*y.derivative(x,1) - x
edo_2ord

$-x - 4 \, D[0]\left(y\right)\left(x\right) + D[0, 0]\left(y\right)\left(x\right)$

desolve(edo_2ord , y)

$-\frac{1}{8} \, x^{2} + K_{1} e^{\left(4 \, x\right)} + K_{2} - \frac{1}{16} \, x - \frac{1}{64}$

Solução por Laplace

y = function('y')(x)

edo_2ord_lap = y.derivative(x,2) - 4*y.derivative(x,1) - x

desolve_laplace(edo_2ord_lap,y )

$-\frac{1}{8} \, x^{2} + \frac{1}{64} \, {\left(16 \, D[0]\left(y\right)\left(0\right) + 1\right)} e^{\left(4 \, x\right)} - \frac{1}{16} \, x + y\left(0\right) - \frac{1}{4} \, D[0]\left(y\right)\left(0\right) - \frac{1}{64}$

Equação diferencial de segunda ordem simbólica com valores iniciais

Solução por desolve()

Resolver a equação diferencial: $y''(x)- 4y'(x) - x$ com $x0 = 0,~ y(x0) = 1 ~e~ y'(x0) = 0$

y = function('y')(x)

edo_2ord =  y.derivative(x,2) - 4*y.derivative(x,1) - x

edo_2ord

$-x - 4 \, D[0]\left(y\right)\left(x\right) + D[0, 0]\left(y\right)\left(x\right)$

# ics =[ x0 , y(x0), y'(x0)]

sol_2ord = desolve ( edo_2ord , y , ics = [0, 1 , 0],show_method = true)
sol_2ord

[$-\frac{1}{8} \, x^{2} - \frac{1}{16} \, x + \frac{1}{64} \, e^{\left(4 \, x\right)} + \frac{63}{64}$, variationofparameters]

Solução por Laplace

y = function('y')(x)

edo_2ord_lap = y.derivative(x,2) - 4*y.derivative(x,1) - x

desolve_laplace(edo_2ord_lap,y,ics =[0,1,0] )

$-\frac{1}{8} \, x^{2} - \frac{1}{16} \, x + \frac{1}{64} \, e^{\left(4 \, x\right)} + \frac{63}{64}$

Gráfico da solução da equação diferencial de segunda ordem

plot(sol_2ord[0],(x,0,1), figsize = (4,3), gridlines="minor");

Campo de direções

reset()

y = function('y')(x)

edo_2ord = y.derivative(x,2) - 4*y.derivative(x,1) - x

sol_2ord = desolve(edo_2ord,y, ics = [0,1,0])
sol_2ord

$-\frac{1}{8} \, x^{2} - \frac{1}{16} \, x + \frac{1}{64} \, e^{\left(4 \, x\right)} + \frac{63}{64}$

plot_campo = plot_slope_field(sol_2ord, (x,0,1.5), (y,0,10), figsize = (4,3))
plot_campo

plot_campo + plot(sol_2ord, (x,0,1.5), figsize = (4,3))

# testando a solução
sol_2ord(1.5).n()

$6.91294989832399$

Equação diferencial de segunda ordem numérica

Para resolver uma equação de segunda ordem numérica é necessário utilizar o recurso de redução da equação a uma de primeira ordem e dessa forma é criado um sistema de equações diferenciais.

Sistema de equações diferenciais

Sistema de Equações diferenciais Simbólico

reset()

ODEs de primeira ordem

$y'(x) + y - 1 ==0$

$w'(x) - x + 1 ==0$

reset()

w = function('w')(x)
y = function('y')(x)

ode1 = w.derivative(x,1) + y - 1 == 0
ode2 = y.derivative(x,1) - x + 1 == 0

desolve_system([ode1, ode2], [y,w], ivar = x)

[$y\left(x\right) = \frac{1}{2} \, x^{2} - x + y\left(0\right)$, $w\left(x\right) = -\frac{1}{6} \, x^{3} + \frac{1}{2} \, x^{2} - \frac{1}{3} \, x {\left(y\left(0\right) - 1\right)} - \frac{2}{3} \, x y\left(0\right) + \frac{2}{3} \, x + w\left(0\right)$]

ODEs de segunda ordem

reset()

w = function('w')(x)
y = function('y')(x)

ode1 = w.derivative(x,2) + y - 1 == 0
ode2 = y.derivative(x,2) - x + 1 == 0

desolve_system([ode1, ode2], [y,w], ivar = x)

[$y\left(x\right) = \frac{1}{6} \, x^{3} - \frac{1}{2} \, x^{2} + x D[0]\left(y\right)\left(0\right) + y\left(0\right)$, $w\left(x\right) = -\frac{1}{120} \, x^{5} + \frac{1}{24} \, x^{4} - \frac{1}{6} \, x^{3} D[0]\left(y\right)\left(0\right) - \frac{1}{4} \, x^{2} {\left(y\left(0\right) - 1\right)} - \frac{1}{4} \, x^{2} y\left(0\right) + \frac{1}{4} \, x^{2} + x D[0]\left(w\right)\left(0\right) + w\left(0\right)$]

Sistema de Equações diferenciais Simbólico com valores iniciais

ODEs de primeira ordem

w = function('w')(x)
y = function('y')(x)

ode1 = diff(w,x) + y - 1 == 0
ode2 = diff(y,x) - x + 1 == 0

desolve_system([ode1, ode2], [y,w], ivar = x, ics=[0,1,2])

[$y\left(x\right) = \frac{1}{2} \, x^{2} - x + 1$, $w\left(x\right) = -\frac{1}{6} \, x^{3} + \frac{1}{2} \, x^{2} + 2$]

ODEs de segunda ordem

reset()

w = function('w')(x)
y = function('y')(x)

ode1 = w.derivative(x,2) + y - 1 == 0
ode2 = y.derivative(x,1)  + 1 == 0

desolve_system([ode1, ode2], [y,w], ivar = x, ics=[0,1,3])

[$y\left(x\right) = -x + 1$, $w\left(x\right) = \frac{1}{6} \, x^{3} + x D[0]\left(w\right)\left(0\right) + 3$]

Sistema de equações diferenciais numérico

ODEs de primeira ordem

Runge Kutta 45

desolve_system_rk4(des, vars, ics=None, ivar=None, end_points=None, step=0.1)

des : [ode1, ode2]

vars : variaveis dependentes

ivar : (optional) should be specified, if there are more variables or if the equation is autonomous and the independent variable is missing

ics : condições iniciais na forma [x0, y01, y02, y03,....]

end_points - the end points of the interval

    if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a)
    if end_points is None, we use end_points=ics[0]+10
    if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b)

step : (optional, default: 0.1) the length of the step

Resolver a equação diferencial: $y'' = 4y'' + x$

Para resolver essa equação diferencial devemos primeiramente fazer uma substituição de variável para reduzir a equação diferencial para uma de primeira ordem e dessa forma criar um sistema de EDOs. Mudança de variavel:

w1 = y
w2 = y'
w3 = y''

Assim, derivando as variáveis acima, teremos um sistema de equações na forma:

w1' = w2
w2' = f(x,w1,w2) = 4*w2 + x

reset()

var('x, w1 , w2 , y')

sol_ode2_rk = desolve_system_rk4([w2, 4*w2 + x ],[w1 , w2], ics = [0,1,0],ivar = x, step = 0.1, end_points = 2)

($x$, $w_{1}$, $w_{2}$, $y$)

x_vet = [i for i, j, k in sol_ode2_rk] # valores do eixo x

y1_sol = [j for i, j, k in sol_ode2_rk] # valores de y solução

y2_sol = [k for i, j, k in sol_ode2_rk] # valores de y' solução

#plotando a sulução de y
line2d(zip(x_vet,y1_sol), gridlines="minor", figsize = (4,3))

Resolver o oscilador de Van der Pol escrito na forma:

$dx/dt$ = y

$dy/dt = u*y(1-x^2) - x$

Sendo $u$ uma constante

reset()

var('t, x, y')

u = 1.0
h = 0.1

Van_der_Pol = [y, -x + u * y * (1 - x^2)]

#ics=[t0, x0, y(x0)]

solOVP = desolve_system_rk4(Van_der_Pol, [x, y], ivar = t , ics = [0, 2, 0] , end_points = 20, step = h)

($t$, $x$, $y$)

x_vet = [i for i, j, k in solOVP] # valores do eixo x

y1_sol = [j for i, j, k in solOVP] # valores de y solução

y2_sol = [k for i, j, k in solOVP] # valores de y' solução

# plotando a solução de y
line2d(zip(x_vet,y1_sol), figsize = (4,3))

# plotando as duas soluções
line2d(zip(y1_sol,y2_sol), figsize = (4,3))

Solução por desolve_odeint

desolve_odeint(des, ics, times, dvars, ivar=None)

des : right hand sides of the system

ics : initial conditions

times : a sequence of time points in which the solution must be found

dvars : dependent variables. ATTENTION: the order must be the same as in des, that means: d(dvars[i])/dt=des[i]

ivar : independent variable, optional.

var('t, x, y')

u = 1.0
h = 0.15
times = srange(0,1,0.1)
h = 0.01

f = y
g = -x + u * y * (1 - x^2)

Van_der_Pol = [f,g ]

sol_odeint = desolve_odeint(Van_der_Pol,[2,0],times,[x,y] )
sol_odeint

($t$, $x$, $y$)

[[ 2. 0. ] [ 1.99093338 -0.17265481] [ 1.96695252 -0.30072111] [ 1.93183209 -0.39739534] [ 1.88817772 -0.47284973] [ 1.83771915 -0.53452342] [ 1.78155292 -0.58775061] [ 1.72032135 -0.63637467] [ 1.65433584 -0.68324222] [ 1.58365714 -0.73057114]]

points((sol_odeint))

** Solução por Método de Euler**

eulers_method_2x2(f, g, t0, x0, y0, h, t1, algorithm='table')

f, g : equações diferenciais

t0, x0, y0 : condições iniciais

h: fragmento

t1 : passo

t, x, y = PolynomialRing(QQ,3,"t,x,y").gens()

f = y
g = -x + u * y * (1 - x^2)

u = 1.0
h = 0.1

passo = 2

# entre as funções e h use t0, x0, y0
eulers_method_2x2(f,g ,0, 2, 0, h, passo)

t x h*f(t,x,y) y h*g(t,x,y) 0 2 0.000000000000000 0 -0.200000000000000 0.100000000000000 2.00000000000000 -0.0200000000000000 -0.200000000000000 -0.140000000000000 0.200000000000000 1.98000000000000 -0.0340000000000000 -0.340000000000000 -0.0987064000000000 0.300000000000000 1.94600000000000 -0.0438706400000000 -0.438706400000000 -0.0723362114537600 0.400000000000000 1.90212936000000 -0.0511042611453760 -0.511042611453760 -0.0564170690908084 0.500000000000000 1.85102509885462 -0.0567459680544568 -0.567459680544568 -0.0474201128039363 0.600000000000000 1.79427913080017 -0.0614879793348505 -0.614879793348505 -0.0429591798438796 0.700000000000000 1.73279115146532 -0.0657838973192384 -0.657838973192384 -0.0415425733257929 0.800000000000000 1.66700725414608 -0.0699381546518177 -0.699381546518177 -0.0422868199436481 0.900000000000000 1.59706909949426 -0.0741668366461825 -0.741668366461825 -0.0447016096559828 1.00000000000000 1.52290226284808 -0.0786369976117808 -0.786369976117808 -0.0485498375252795 1.10000000000000 1.44426526523630 -0.0834919813643088 -0.834919813643088 -0.0537624039206839 1.20000000000000 1.36077328387199 -0.0888682217563772 -0.888682217563772 -0.0603879146563158 1.30000000000000 1.27190506211561 -0.0949070132220087 -0.949070132220087 -0.0685624118267037 1.40000000000000 1.17699804889360 -0.101763254404679 -1.01763254404679 -0.0784879392213772 1.50000000000000 1.07523479448892 -0.109612048326817 -1.09612048326817 -0.0904097653298254 1.60000000000000 0.965622746162107 -0.118653024859799 -1.18653024859799 -0.104579981304186 1.70000000000000 0.846969721302307 -0.129111022990218 -1.29111022990218 -0.121189207486986 1.80000000000000 0.717858698312089 -0.141229943738917 -1.41229943738917 -0.140237042092511 1.90000000000000 0.576628754573173 -0.155253647948168 -1.55253647948168 -0.161294573586846 2.00000000000000 0.421375106625005 -0.171383105306852 -1.71383105306852 -0.183090349285228

ODEs de segunda ordem

%%% FIM EQUAÇÕES DIFERENCIAIS %%%