Método do shooting não linear¶

Prof. Doherty Andrade

Script baseado no algoritmo de Burden-Faires.

O método do shooting para problemas de contorno de segunda ordem não linear $$ PVF (1) \begin{cases} y^{\prime\prime}= f(x,y,y^\prime), a\leq x \leq b\cr y(a)=\alpha,\quad \mbox{ e }\cr y(b)=\beta. \end{cases} $$

é semelhante ao método do shooting linear, exceto que a solução para problemas não lineares não pode ser expressa como uma combinação linear de soluções de dois problemas de valor inicial. Em vez disso, aproximamos a solução para o problema de contorno pela utilização de soluções de uma sequência de problemas de valor inicial envolvendo um parâmetro $t$.

Esses problemas têm a forma $$ PVI (2) \begin{cases} y^{\prime\prime}= f(x,y,y^\prime), a\leq x \leq b\cr y(a)=\alpha,\quad \mbox{ e }\cr y^\prime(a)=t. \end{cases} $$

Fazemos isso pela escolha dos parâmetros $t = t_k$ de modo a assegurar que $$ \lim_{k \to \infty}y(b,t_k)= y(b)= \beta,$$ em que $y(x,t_k)$ denota a solução para o problema de valor inicial (2) com $t=t_k$ e $y(x)$ denota a solução problema de contorno (1).

O resultado que assegura a existência de solução para o PFV(2) é dado pelo seguinte teorema:

TEOREMA: Suponha que a função $f$ no PVF (2) seja contínua no conjunto $$ D = \{ (x,y,y^\prime); a \leq x \leq b, -\infty < y < \infty , -\infty < y^\prime < \infty \}.$$ e que as derivadas parciais $f_x$ e $f_y$ também sejam contínuas em $D$. Se

(i) $f_x(x,y,y^\prime)> 0 $ em $D$, e

(2i) existir contante $M\geq 0$ tal que $$ \vert f(x,y,y^\prime)\vert \leq M, \forall (x,y,y^\prime )\in D$$

então o PFV (2) tem uma única solução.

Para determinar o parâmetro $t_k$, suponha que a função $f$ satisfaz às condições do teorema acima. Se $y(x,t)$ denota a solução do PVI (2), em seguida determinamos $t$ como $$y(b,t)-\beta = 0.$$ Esta é uma equação não-linear e devemos utilizar um dos métodos conhecidos para resolvê-la.

Se usarmos o método da secante, precisaremos de duas aproximações iniciais $t_0$ e $t_1$ e gerar os termos restantes da sequência por $$t_k= t_{k-1} - \frac{\left(y(b,t_{k-1})-\beta \right)(t_{k-1}-t_{k-2})}{y(b,t_{k-1})-y(b,t_{k-2})}, k=2,3,\ldots.$$

Se usarmos o método de Newton-Raphson para gerar a sequência $t_k$, precisaremos de uma aproximação inciail $t_0$ e gerar os demais termos da sequência por

\begin{equation} \mbox{(MNR) } t_k = t_{k-1} - \frac{y(b,t_{k-1})-\beta }{\frac{dy}{dt}(b,t_{k-1})}, k=1,2,3,\ldots.\end{equation}

Note que este enfoque exige o conhecimento da derivada $\frac{dy}{dt}(b,t_{k-1})$ o que é uma dificuldadeuma vez que uma representação explícita de $y(b,t)$ não é conhecida, conhecendo-se apenas valores $y(b,t_0),$ $y(b,t_1)$,...,$y(b,t_{k-1})$.

Para contornar esta dificuldade vamos derivar com relação $t$ a equação $$y^{\prime\prime}= f(x,y,y^\prime)$$ do PVI (2). Usando que $x$ e $t$ são independentes, $\frac{\partial x}{\partial t} =0$, usando as condições iniciais obtemos que $$\frac{\partial y}{\partial t}(a,t)=0$$ e $$\frac{\partial y^\prime }{\partial t}(a,t)=1,$$ podemos simplificar e obter $$\frac{\partial y^{\prime \prime}}{\partial t} (x,t)= \frac{\partial f}{\partial y} (x,y, y^\prime)\frac{\partial y}{\partial t} (x,t)+ \frac{\partial f}{\partial y^\prime} (x,y, y^\prime) \frac{\partial y^\prime}{\partial t}(x,t) $$ com as condições iniciais $$\frac{\partial y}{\partial t}(a,t) =0$$ e $$\frac{\partial y^\prime}{\partial t}(a,t) =1.$$

Para simplificar a notação, vamos usar $z(x,t)$ para denotar $\frac{\partial y}{\partial t}(x,t)$. Admitindo que a ordem de derivação com relação a $x$ e a $t$ pode ser trocada, a equação acima fica

$$ PVI (3) \begin{cases} z^{\prime\prime}(x,t) = \frac{\partial f}{\partial y} (x,y, y^\prime)z (x,t)+ \frac{\partial f}{\partial y^\prime} (x,y, y^\prime) z^\prime(x,t) , a \leq x \leq b,\cr z(a,t)= 0 , \cr z^\prime(a,t)=1. \end{cases} $$

Assim, o método de Newton-Raphson exige que dois problemas de valor inicial sejam resolvidos a cada iteração, PVI(2) e PVI(3). Então da equação (MNR) fica \begin{equation} t_k = t_{k-1} - \frac{y(b,t_{k-1}) - \beta}{z(b,t_{k-1})}. \end{equation}

Os problemas de valor inicial não podem, em geral, serem resolvidos explicitamente, portanto, há a necessidade de algum método numérico que pode ser Runge-Kutta de ordem 4 para aproximar ambas as soluções que são necessárias no método de Newton-Raphson.

Observe que cada PVI deve ser escrito como um sistema de EDOS de primeira ordem e então utilizar o Runge-Kutta.

Exemplo¶

Considere o PVF dado a seguir:

\begin{cases} y^{\prime \prime} = \frac{1}{8}\left(32+2x^3-yy^\prime \right), 1 \leq x \leq 3\cr y(1) = 17 ,\cr y(3) = \frac{43}{3} \end{cases}

cuja solução exata é $y(x) = x^2+ \frac{16}{x}$.

Vamos aplicar o método do shooting não linear programada a seguir.

In [14]:
import numpy as np  # biblioteca básica de matemática
import sympy as sp
from sympy import Symbol
import scipy.linalg   # SciPy biblioteca de algebra linear
import math
In [15]:
from numpy import zeros, abs

def shoot_nonlinear(a,b,alpha, beta, n, tol, M):

    w1 = zeros(n+1)  
    w2 = zeros(n+1)
    h = (b-a)/n
    k = 1
    TK = (beta - alpha)/(b - a)

    print("i""    x" "      " "W1""      " "W2")
    while k <= M:

        w1[0] = alpha
        w2[0] = TK
        u1    = 0
        u2    = 1

        for i in range(1,n+1):
            x = a + (i-1)*h     #step 5

            t = x + 0.5*(h)

            k11 = h*w2[i-1]     #step 6

            k12 = h*f(x,w1[i-1],w2[i-1])
            k21 = h*(w2[i-1] + (1/2)*k12)
            k22 = h*f(t, w1[i-1] + (1/2)*k11, w2[i-1] + (1/2)*k12)
            k31 = h*(w2[i-1] + (1/2)*k22)
            k32 = h*f(t, w1[i-1] + (1/2)*k21, w2[i-1] + (1/2)*k22)
            t   = x + h
            k41 = h*(w2[i-1]+k32)
            k42 = h*f(t, w1[i-1] + k31, w2[i-1] + k32)
            w1[i] = w1[i-1] + (k11 + 2*k21 + 2*k31 + k41)/6
            w2[i] = w2[i-1] + (k12 + 2*k22 + 2*k32 + k42)/6   
            kp11 = h*u2
            kp12 = h*(fy(x,w1[i-1],w2[i-1])*u1 + fyp(x,w1[i-1], w2[i-1])*u2)
            t    = x + 0.5*(h)
            kp21 = h*(u2 + (1/2)*kp12)
            kp22 = h*((fy(t, w1[i-1],w2[i-1])*(u1 + (1/2)*kp11)) + fyp(x+h/2, w1[i-1],w2[i-1])*(u2 +(1/2)*kp12))
            kp31 = h*(u2 + (1/2)*kp22)
            kp32 = h*((fy(t, w1[i-1],w2[i-1])*(u1 + (1/2)*kp21)) + fyp(x+h/2, w1[i-1],w2[i-1])*(u2 +(1/2)*kp22))
            t    = x + h
            kp41 = h*(u2 + kp32)
            kp42 = h*(fy(t, w1[i-1], w2[i-1])*(u1+kp31) + fyp(x + h, w1[i-1], w2[i-1])*(u2 + kp32))
            u1 = u1 + (1/6)*(kp11 + 2*kp21 + 2*kp31 + kp41)
            u2 = u2 + (1/6)*(kp12 + 2*kp22 + 2*kp32 + kp42)


        r = abs(w1[n] - beta)
        #print(r)
        if r < tol:
            for i in range(0,n+1):
                x = a + i*h
                print("%.2f %.2f %.4f %.4f" %(i,x,w1[i],w2[i]))
            return

        TK = TK -(w1[n]-beta)/u1

        k = k+1


    print("Maximum number of iterations exceeded")   
    return

Exemplo 1¶

In [17]:
def f(x,y,yp):
    fx = (1/8)*(32 + 2*x**3 -y*yp)
    return fx

def fy(xp,z,zp):
    fyy = -(1/8)*(zp)
    return fyy

def fyp(xpp,zpp,zppp):
    fypp = -(1/8)*(zpp)
    return fypp

a = 1         # start point
b = 3         # end point
alpha = 17    # boundary condition
beta = 43/3   # boundary condition
N = 20        # number of subintervals
M = 10        # maximum number of iterations
tol = 0.00001 # tolerance


shoot_nonlinear(a,b,alpha,beta,N,tol,M)

i    x      W1      W2
0.00 1.00 17.0000 -14.0002
1.00 1.10 15.7555 -11.0233
2.00 1.20 14.7734 -8.7113
3.00 1.30 13.9978 -6.8676
4.00 1.40 13.3886 -5.3634
5.00 1.50 12.9167 -4.1112
6.00 1.60 12.5600 -3.0501
7.00 1.70 12.3018 -2.1364
8.00 1.80 12.1289 -1.3384
9.00 1.90 12.0311 -0.6322
10.00 2.00 12.0000 -0.0001
11.00 2.10 12.0291 0.5718
12.00 2.20 12.1127 1.0942
13.00 2.30 12.2465 1.5754
14.00 2.40 12.4267 2.0222
15.00 2.50 12.6500 2.4400
16.00 2.60 12.9138 2.8331
17.00 2.70 13.2159 3.2052
18.00 2.80 13.5543 3.5592
19.00 2.90 13.9272 3.8975
20.00 3.00 14.3333 4.2222

Nova Versão do código com grafico¶

In [8]:
# versao com grafico
from numpy import zeros, abs, linspace
import matplotlib.pyplot as plt

def shoot_nonlinear(a, b, alpha, beta, n, tol, M):
    w1 = zeros(n+1)  
    w2 = zeros(n+1)
    h = (b-a)/n
    k = 1
    TK = (beta - alpha)/(b - a)

    print("i   x        W1        W2")
    while k <= M:
        w1[0] = alpha
        w2[0] = TK
        u1 = 0
        u2 = 1

        for i in range(1, n+1):
            x = a + (i-1)*h
            t = x + 0.5*(h)

            k11 = h*w2[i-1]
            k12 = h*f(x, w1[i-1], w2[i-1])
            k21 = h*(w2[i-1] + (1/2)*k12)
            k22 = h*f(t, w1[i-1] + (1/2)*k11, w2[i-1] + (1/2)*k12)
            k31 = h*(w2[i-1] + (1/2)*k22)
            k32 = h*f(t, w1[i-1] + (1/2)*k21, w2[i-1] + (1/2)*k22)
            t = x + h
            k41 = h*(w2[i-1]+k32)
            k42 = h*f(t, w1[i-1] + k31, w2[i-1] + k32)
            w1[i] = w1[i-1] + (k11 + 2*k21 + 2*k31 + k41)/6
            w2[i] = w2[i-1] + (k12 + 2*k22 + 2*k32 + k42)/6   
            
            kp11 = h*u2
            kp12 = h*(fy(x, w1[i-1], w2[i-1])*u1 + fyp(x, w1[i-1], w2[i-1])*u2)
            t = x + 0.5*(h)
            kp21 = h*(u2 + (1/2)*kp12)
            kp22 = h*((fy(t, w1[i-1], w2[i-1])*(u1 + (1/2)*kp11)) + fyp(x+h/2, w1[i-1], w2[i-1])*(u2 + (1/2)*kp12))
            kp31 = h*(u2 + (1/2)*kp22)
            kp32 = h*((fy(t, w1[i-1], w2[i-1])*(u1 + (1/2)*kp21)) + fyp(x+h/2, w1[i-1], w2[i-1])*(u2 + (1/2)*kp22))
            t = x + h
            kp41 = h*(u2 + kp32)
            kp42 = h*(fy(t, w1[i-1], w2[i-1])*(u1+kp31) + fyp(x + h, w1[i-1], w2[i-1])*(u2 + kp32))
            u1 = u1 + (1/6)*(kp11 + 2*kp21 + 2*kp31 + kp41)
            u2 = u2 + (1/6)*(kp12 + 2*kp22 + 2*kp32 + kp42)

        r = abs(w1[n] - beta)
        if r < tol:
            x_points = [a + i*h for i in range(0, n+1)]
            for i in range(0, n+1):
                print("%.2f %.2f %.4f %.4f" % (i, x_points[i], w1[i], w2[i]))
            
            # Plotando o gráfico
            plt.figure(figsize=(10, 6))
            plt.plot(x_points, w1, 'b-', linewidth=2, label='Solução y(x)')
            plt.xlabel('x')
            plt.ylabel('y(x)')
            plt.title('Solução Numérica pelo Método do Shooting Não Linear')
            plt.grid(True)
            plt.legend()
            plt.show()
            
            return w1, x_points

        TK = TK - (w1[n]-beta)/u1
        k = k+1

    print("Maximum number of iterations exceeded")   
    return None, None

# Você precisará definir as funções f, fy e fyp antes de chamar shoot_nonlinear
# Exemplo:
# def f(x, y, y_prime):
#     return ...  # sua equação diferencial aqui
# 
# def fy(x, y, y_prime):
#     return ...  # derivada parcial de f em relação a y
# 
# def fyp(x, y, y_prime):
#     return ...  # derivada parcial de f em relação a y'

# Exemplo de chamada:
# a, b = 0, 1
# alpha, beta = 1, 2
# n = 10
# tol = 1e-6
# M = 100
# shoot_nonlinear(a, b, alpha, beta, n, tol, M)

Exemplo de Bratu¶

O problema de Bratu é um caso clássico em equações diferenciais não lineares, frequentemente usado como teste para métodos numéricos, como o método do shooting não linear que estamos discutindo. Ele aparece em problemas de combustão, teoria de reações químicas e transferência de calor.

In [9]:
# exemplo de bratu

from numpy import zeros, abs, linspace, exp
import matplotlib.pyplot as plt

def f(x, y, y_prime):
    return -exp(y)  # Equação de Bratu: y'' = -e^y

def fy(x, y, y_prime):
    return -exp(y)

def fyp(x, y, y_prime):
    return 0

# Configuração do problema
a, b = 0, 1
alpha, beta = 0, 0
n = 30
tol = 1e-6
M = 100

# Resolver
solution, x_points = shoot_nonlinear(a, b, alpha, beta, n, tol, M)

# Plotagem
if solution is not None:
    plt.figure(figsize=(10, 6))
    plt.plot(x_points, solution, 'r-o', linewidth=2, markersize=4)
    plt.title('Solução da Equação de Bratu')
    plt.xlabel('x')
    plt.ylabel('y(x)')
    plt.grid(True)
    plt.show()

i   x        W1        W2
0.00 0.00 0.0000 0.5494
1.00 0.03 0.0178 0.5157
2.00 0.07 0.0344 0.4815
3.00 0.10 0.0498 0.4467
4.00 0.13 0.0642 0.4114
5.00 0.17 0.0773 0.3757
6.00 0.20 0.0892 0.3394
7.00 0.23 0.0999 0.3028
8.00 0.27 0.1094 0.2658
9.00 0.30 0.1176 0.2284
10.00 0.33 0.1246 0.1908
11.00 0.37 0.1303 0.1529
12.00 0.40 0.1348 0.1149
13.00 0.43 0.1380 0.0767
14.00 0.47 0.1399 0.0384
15.00 0.50 0.1405 0.0000
16.00 0.53 0.1399 -0.0384
17.00 0.57 0.1380 -0.0767
18.00 0.60 0.1348 -0.1149
19.00 0.63 0.1303 -0.1529
20.00 0.67 0.1246 -0.1908
21.00 0.70 0.1176 -0.2284
22.00 0.73 0.1094 -0.2658
23.00 0.77 0.0999 -0.3028
24.00 0.80 0.0892 -0.3394
25.00 0.83 0.0773 -0.3757
26.00 0.87 0.0642 -0.4114
27.00 0.90 0.0498 -0.4467
28.00 0.93 0.0344 -0.4815
29.00 0.97 0.0178 -0.5157
30.00 1.00 -0.0000 -0.5494
No description has been provided for this image
No description has been provided for this image
In [ ]: