Integração Numérica

Existem diversos métodos comumente utilizados para encontrar soluções aproximadas para sistemas de EDOs. Neste documento introduziremos os mais usados assim com exemplos do seu uso no Sage.

Método de Euler

No método de Euler aproximamos a EDO de primeira ordem por uma série de pontos no tempo a intervalos regulares de tamanho $h$ definido pelo usuário. Logo, $t_n=t_0 + n h$ . Assumindo um EDO da forma $\frac{dy}{dt}=f(t,y(t))$ , a solução aproximada pelo método de Euler é $y_{n+1} = y_n + hf(t_n,y_n)$ .

Abaixo temos um exemplo interativo para um sistema relativamente simples: $\frac{dx}{dt} = sin(x)$

In [7]:

%display typeset

In [2]:

def tab_list(y, headers = None):
    '''
    Converte uma lista em uma tabela html
    '''
    s = '<table border = 1>'
    if headers:
        for q in headers:
            s = s + '<th>' + str(q) + '</th>'
    for x in y:
        s = s + '<tr>'
        for q in x:
            s = s + '<td>' + str(q) + '</td>'
        s = s + '</tr>'
    s = s + '</table>'
    return s

var('x y')
@interact
def euler_method(y_exact_in = input_box('-cos(x)+1.0', type = str, label = 'Solução Exata = '), y_prime_in = input_box('sin(x)', type = str, label = "y' = "), start = input_box(0.0, label = 'Val. inicial de x: '), stop = input_box(6.0, label = 'Valor de parada de x: '), startval = input_box(0.0, label = 'Valor inicial de y: '), nsteps = slider([2^m for m in range(0,10)], default = 10, label = 'Número de passos: '), show_steps = slider([2^m for m in range(0,10)], default = 8, label = 'Numero de passos mostrados na tabela: ')):
    y_exact = lambda x: eval(y_exact_in)
    y_prime = lambda x,y: eval(y_prime_in)
    stepsize = float((stop-start)/nsteps)
    steps_shown = min(nsteps,show_steps)
    sol = [startval]
    xvals = [start]
    for step in range(nsteps):
        sol.append(sol[-1] + stepsize*y_prime(xvals[-1],sol[-1]))
        xvals.append(xvals[-1] + stepsize)
    sol_max = max(sol + [find_local_maximum(y_exact,start,stop)[0]])
    sol_min = min(sol + [find_local_minimum(y_exact,start,stop)[0]])
    show(plot(y_exact(x),start,stop,rgbcolor=(1,0,0))+line([[xvals[index],sol[index]] for index in range(len(sol))]),xmin=start,xmax = stop, ymax = sol_max, ymin = sol_min)
    if nsteps < steps_shown:
        table_range = range(len(sol))
    else:
        table_range = list(range(0,floor(steps_shown/2))) + list(range(len(sol)-floor(steps_shown/2),len(sol)))
    html(tab_list([[i,xvals[i],sol[i]] for i in table_range], headers = ['step','x','y']))

Campos Vetoriais e o Método de Euler

In [8]:

x,y = var('x,y')
from sage.ext.fast_eval import fast_float
@interact
def _(f = input_box(default=y), g=input_box(default=-x*y+x^3-x),
      xmin=input_box(default=-1), xmax=input_box(default=1),
      ymin=input_box(default=-1), ymax=input_box(default=1),
      x_ini=input_box(default=0.5), y_ini=input_box(default=0.5),
      tam_passo=(0.01,(0.001, 0.2, .001)), passos=(600,(0, 1400)) ):
    ff = fast_float(f, 'x', 'y')
    gg = fast_float(g, 'x', 'y')
    passos = int(passos)
    
    points = [ (x_ini, y_ini) ]
    for i in range(passos):
        xx, yy = points[-1]
        try:
            points.append( (xx + tam_passo * ff(xx,yy), yy + tam_passo * gg(xx,yy)) )
        except (ValueError, ArithmeticError, TypeError):
            break
    equilibria = solve([f,g], [x,y])
    xnull = plot(solve(f,y)[0].rhs(),(x,xmin,xmax), color='red', legend_label="Nuliclina de X")
    ynull = plot(solve(g,y)[0].rhs(),(x,xmin,xmax), color='green', legend_label="Nuliclina de Y")
    starting_point = point(points[0], pointsize=50)
    solution = line(points)
    vector_field = plot_vector_field( (f,g), (x,xmin,xmax), (y,ymin,ymax) )

    result = vector_field + starting_point + solution +xnull+ ynull

    pretty_print(html(r"$\displaystyle\frac{dx}{dt} = %s$ <br>$\displaystyle\frac{dy}{dt} = %s$" % (latex(f),latex(g))))
    result.show(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax)

In [ ]:

Campo Vetorial com Runge-Kutta Fehlberg

In [9]:

var('x y')
@interact
def _(fin = input_box(default=y+exp(x/10)-1/3*((x-1/2)^2+y^3)*x-x*y^3), gin=input_box(default=x^3-x+1/100*exp(y*x^2+x*y^2)-0.7*x),
      xmin=input_box(default=-1), xmax=input_box(default=1.8),
      ymin=input_box(default=-1.3), ymax=input_box(default=1.5),
      x_start=(-1,(-2,2)), y_start=(0,(-2,2)), error=(0.5,(0,1)),
      t_length=(23,(0, 100)) , num_of_points = (1500,(5,2000)),
      algorithm = selector([
         ("rkf45" , "runga-kutta-felhberg (4,5)"),
         ("rk2" , "embedded runga-kutta (2,3)"),
         ("rk4" , "4th order classical runga-kutta"),
         ("rk8pd" , 'runga-kutta prince-dormand (8,9)'),
         ("rk2imp" , "implicit 2nd order runga-kutta at gaussian points"),
         ("rk4imp" , "implicit 4th order runga-kutta at gaussian points"),
         ("bsimp" , "implicit burlisch-stoer (requires jacobian)"),
         ("gear1" , "M=1 implicit gear"),
         ("gear2" , "M=2 implicit gear")
      ])):
    f(x,y)=fin
    g(x,y)=gin

    ff = f._fast_float_(*f.args())
    gg = g._fast_float_(*g.args())

    #solve
    path = []
    err = error
    xerr = 0
    for yerr in [-err, 0, +err]:
      T=ode_solver()
      T.algorithm=algorithm
      T.function = lambda t, yp: [ff(yp[0],yp[1]), gg(yp[0],yp[1])]
      T.jacobian = jacobian([f,g],[x,y])
      T.ode_solve(y_0=[x_start + xerr, y_start + yerr],t_span=[0,t_length],num_points=num_of_points)
      path.append(line([p[1] for p in T.solution]))

    #plot
    vector_field = plot_vector_field( (f,g), (x,xmin,xmax), (y,ymin,ymax) )
    starting_point = point([x_start, y_start], pointsize=50)
    xnull = plot(solve(f,y)[0].rhs(),(x,xmin,xmax), color='red', legend_label="Nuliclina de X")
    ynull = plot(solve(g,y)[0].rhs(),(x,xmin, xmax), color='green', legend_label="Nuliclina de Y")
    show(vector_field + starting_point + sum(path)+xnull, aspect_ratio=1, figsize=[8,9])

In [5]:

var('x y z')

sol = solve(-x*y^3 - 1/12*(4*y^3 + (2*x - 1)^2)*x + y + e^(1/10*x), y)
sol[0].simplify_full()

\renewcommand{\Bold}[1]{\mathbf{#1}}y = -\frac{4^{\frac{2}{3}} {\left(x {\left(i \, \sqrt{3} + 1\right)} \left(-\frac{4 \, x^{3} - 4 \, x^{2} + x - \sqrt{\frac{16 \, x^{7} - 32 \, x^{6} + 24 \, x^{5} - 8 \, x^{4} + x^{3} + 144 \, x e^{\left(\frac{1}{5} \, x\right)} - 24 \, {\left(4 \, x^{4} - 4 \, x^{3} + x^{2}\right)} e^{\left(\frac{1}{10} \, x\right)} - 16}{x}} - 12 \, e^{\left(\frac{1}{10} \, x\right)}}{x}\right)^{\frac{2}{3}} - i \cdot 4^{\frac{2}{3}} \sqrt{3} + 4^{\frac{2}{3}}\right)}}{16 \, x \left(-\frac{4 \, x^{3} - 4 \, x^{2} + x - \sqrt{\frac{16 \, x^{7} - 32 \, x^{6} + 24 \, x^{5} - 8 \, x^{4} + x^{3} + 144 \, x e^{\left(\frac{1}{5} \, x\right)} - 24 \, {\left(4 \, x^{4} - 4 \, x^{3} + x^{2}\right)} e^{\left(\frac{1}{10} \, x\right)} - 16}{x}} - 12 \, e^{\left(\frac{1}{10} \, x\right)}}{x}\right)^{\frac{1}{3}}}

In [6]:

complex_plot(sol[0].subs(x=z), (x,-1, 1), (y,-1,1))

<ipython-input-6-9dff5e30990a>:1: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...)
See http://trac.sagemath.org/5930 for details.
  complex_plot(sol[Integer(0)].subs(x=z), (x,-Integer(1), Integer(1)), (y,-Integer(1),Integer(1)))

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-6-9dff5e30990a> in <module>
----> 1 complex_plot(sol[Integer(0)].subs(x=z), (x,-Integer(1), Integer(1)), (y,-Integer(1),Integer(1)))

~/Downloads/SageMath/local/lib/python3.9/site-packages/sage/misc/lazy_import.pyx in sage.misc.lazy_import.LazyImport.__call__ (build/cythonized/sage/misc/lazy_import.c:4032)()
    358             True
    359         """
--> 360         return self.get_object()(*args, **kwds)
    361 
    362     def __repr__(self):
~/Downloads/SageMath/local/lib/python3.9/site-packages/sage/misc/decorators.py in wrapper(*args, **kwds)
    489                 options['__original_opts'] = kwds
    490             options.update(kwds)
--> 491             return func(*args, **options)
    492 
    493         #Add the options specified by @options to the signature of the wrapped
~/Downloads/SageMath/local/lib/python3.9/site-packages/sage/plot/complex_plot.pyx in sage.plot.complex_plot.complex_plot (build/cythonized/sage/plot/complex_plot.c:6071)()
    402     g = Graphics()
    403     g._set_extra_kwds(Graphics._extract_kwds_for_show(options, ignore=['xmin', 'xmax']))
--> 404     g.add_primitive(ComplexPlot(complex_to_rgb(z_values),
    405                                 x_range[:2], y_range[:2], options))
    406     return g
~/Downloads/SageMath/local/lib/python3.9/site-packages/sage/plot/complex_plot.pyx in sage.plot.complex_plot.complex_to_rgb (build/cythonized/sage/plot/complex_plot.c:3994)()
    113                 z = <ComplexDoubleElement>zz
    114             else:
--> 115                 z = CDF(zz)
    116             x, y = z._complex.dat
    117             mag = hypot(x, y)
~/Downloads/SageMath/local/lib/python3.9/site-packages/sage/rings/complex_double.pyx in sage.rings.complex_double.ComplexDoubleField_class.__call__ (build/cythonized/sage/rings/complex_double.c:5419)()
    334         if im is not None:
    335             x = x, im
--> 336         return Parent.__call__(self, x)
    337 
    338     def _element_constructor_(self, x):
~/Downloads/SageMath/local/lib/python3.9/site-packages/sage/structure/parent.pyx in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9335)()
    896         if mor is not None:
    897             if no_extra_args:
--> 898                 return mor._call_(x)
    899             else:
    900                 return mor._call_with_args(x, args, kwds)
~/Downloads/SageMath/local/lib/python3.9/site-packages/sage/structure/coerce_maps.pyx in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4622)()
    159                 print(type(C), C)
    160                 print(type(C._element_constructor), C._element_constructor)
--> 161             raise
    162 
    163     cpdef Element _call_with_args(self, x, args=(), kwds={}):
~/Downloads/SageMath/local/lib/python3.9/site-packages/sage/structure/coerce_maps.pyx in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4514)()
    154         cdef Parent C = self._codomain
    155         try:
--> 156             return C._element_constructor(x)
    157         except Exception:
    158             if print_warnings:
~/Downloads/SageMath/local/lib/python3.9/site-packages/sage/rings/complex_double.pyx in sage.rings.complex_double.ComplexDoubleField_class._element_constructor_ (build/cythonized/sage/rings/complex_double.c:5977)()
    366                 return t
    367         elif hasattr(x, '_complex_double_'):
--> 368             return x._complex_double_(self)
    369         else:
    370             return ComplexDoubleElement(x, 0)
~/Downloads/SageMath/local/lib/python3.9/site-packages/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression._complex_double_ (build/cythonized/sage/symbolic/expression.cpp:11835)()
   1663             0.5000000000000001 + 0.8660254037844386*I
   1664         """
-> 1665         return self._eval_self(R)
   1666 
   1667     def __float__(self):
~/Downloads/SageMath/local/lib/python3.9/site-packages/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression._eval_self (build/cythonized/sage/symbolic/expression.cpp:10179)()
   1430             return R(ans)
   1431         else:
-> 1432             raise TypeError("Cannot evaluate symbolic expression to a numeric value.")
   1433 
   1434     cpdef _convert(self, kwds):
TypeError: Cannot evaluate symbolic expression to a numeric value.

In [ ]: