
# isto é um comentário. abaixo, fazemos a soma de 1 com 1
1+1

2


232323*8987

2087886801


11/4

2.75


11//4

2


11 % 4

3


import sympy as sp
sp.isprime(11)

True


def checaprimo(n):
    if n == 2:
        return True
    if n % 2 == 0 or n <= 1:
        return False
    for divisor in range(3, n, 2):
        if n % divisor == 0:
            return False
    return True
checaprimo(101)

True


import math
def checaprimo2(n):
    if n == 2:
        return True
    if n % 2 == 0 or n <= 1:
        return False
    for divisor in range(3, int(math.sqrt(n)), 2):
        if n % divisor == 0:
            return False
    return True
checaprimo2(101)

True


def minharaizquadrada(x):
    y=1
    # aqui é que definimos o numero de iterações, troque 10 pelo que quiser;
    # quando maior, melhor aproximação
    for k in range(1,10):
        y = (y + x/y)/2
    return y

# calculando sqrt(2)
minharaizquadrada(2)

1.414213562373095


def checaprimo3(n):
    if n == 2:
        return True
    if n % 2 == 0 or n <= 1:
        return False
    for divisor in range(3, int(minharaizquadrada(n)), 2):
        if n % divisor == 0:
            return False
    return True
checaprimo3(101)

True


import sympy as sp


x, y = sp.symbols('x y')


sp.solve(x+1)

[-1]


sp.solve(x**2-9)

[-3, 3]


sp.solve(sp.Eq(x**2,9))

[-3, 3]


sp.solve(x**7-x**6+x**3-9)

[CRootOf(x**7 - x**6 + x**3 - 9, 0),
 CRootOf(x**7 - x**6 + x**3 - 9, 1),
 CRootOf(x**7 - x**6 + x**3 - 9, 2),
 CRootOf(x**7 - x**6 + x**3 - 9, 3),
 CRootOf(x**7 - x**6 + x**3 - 9, 4),
 CRootOf(x**7 - x**6 + x**3 - 9, 5),
 CRootOf(x**7 - x**6 + x**3 - 9, 6)]


sp.nsolve(x**7-x**6+x**3-9,1)


x = sp.symbols('x')
sp.plot(x**7-x**6+x**3-9,(x,-2,2),ylim=[-50,50])

<sympy.plotting.plot.Plot at 0x7f79b89c98e0>


sp.solve((x**2+y**2+x-10,x+y-2),(x,y),dict=True)

[{x: 3/4 - sqrt(57)/4, y: 5/4 + sqrt(57)/4},
 {x: 3/4 + sqrt(57)/4, y: 5/4 - sqrt(57)/4}]


sp.nsolve((x**2+y**2+x-10,x+y-2),(x,y),[1,1],prec=30)


sp.nsolve((x**2+y**2+x-10,x+y-2),(x,y),[1,1],dict=True,prec=30)

[{x: 2.63745860881768742430917120174, y: -0.637458608817687424309171201737}]


x, y = sp.symbols('x y')
F=x**2+y**3-2
F.subs(x,1).subs(y,2)


k=F.subs(x,1).subs(y,2)
k**2


sp.linsolve((x+y-1,x-y-3),(x,y))


sp.solve((x+y-1,x-y-3),(x,y))

{x: 2, y: -1}


sp.nsolve((x+y-1,x-y-3),(x,y),[1,1],prec=5)


sp.nsolve((x**5+y**6+x**3-1,x+y-1),(x,y),[0.5,0],dict=True)

[{x: 0.837615759482675, y: 0.162384240517325}]


sp.nsolve((x**5+y**6+x**3-1,x+y-1),(x,y),[1,2],dict=True)

[{x: -1.12228828894446e-23, y: 1.00000000000000}]


sp.plot_implicit(sp.Or(sp.Eq(x**5+y**6+x**3,1),sp.Eq(x+y,1)))

<sympy.plotting.plot.Plot at 0x7f79d8da7df0>


# a funcao vai pegar uma funcao f, extremos de um intervalo a, b e N será o
# número máximo de iterações que serão realizadas
def bisseccao(f,a,b,N):

    if f(a)*f(b) >= 0:
        print("Falha, escolha outro intervalo.")
        return None
    a_n = a
    b_n = b
    for n in range(1,N+1):
        m_n = (a_n + b_n)/2
        f_m_n = f(m_n)
        if f(a_n)*f_m_n < 0:
            a_n = a_n
            b_n = m_n
        elif f(b_n)*f_m_n < 0:
            a_n = m_n
            b_n = b_n
        elif f_m_n == 0:
            print("Solução encontrada!")
            return m_n
        else:
            print("O método falhou.")
            return None
    return (a_n + b_n)/2

# vamos testar com o nosso polinômio, fazendo 50 iterações no intervalo [0,1]
f = lambda x: x**7 +x - 1
bisseccao(f,0,1,50)

Solução encontrada!

0.7965443541284571


import sympy as sp
x = sp.symbols('x')
sp.plot(x**7+x-1,(x,0,1),ylim=[-1,1])

<sympy.plotting.plot.Plot at 0x7f79b8fe2730>


# importando as bibliotecas
import numpy as np
import matplotlib.pyplot as plt

# definindo a funcao
def f(x):
    return x**2+1 

# definindo o intervalo de x e quantos pontos serão criados.
# abaixo definimos que o gráfico será criado com base em 50 pontos
# com coordenada x igualmente distribuída entre -3 e 3.
# a variável x armazenará uma lista com esses 50 pontos.
# esse processo é conhecido por "discretização".
x = np.linspace(-3, 3, 50)

# definicao de y=f(x). na verdade, o que esse comando faz é criar uma
# lista aplicando a funcao f(x) que definimos acima nos pontos da lista
# que contém os valores de x (os tais 50 valores entre -3 e 3)
y = f(x)

# plotando a funcao y=f(x).
plt.plot(x, y)

[<matplotlib.lines.Line2D at 0x7f79b913e9d0>]


x = np.linspace(-3, 3, 4)
y=f(x)
plt.plot(x, y)

[<matplotlib.lines.Line2D at 0x7f79b929fbb0>]


# vamos chamar de "a" um primeiro plot vazio, só com os eixos.
a = plt.axes()

# agora damos nomes para os eixos do gráfico vazio.
a.set_xlabel('x');
a.set_ylabel('y');

# delimitando a área da janela de visualização para ser com
# x entre -2 e 2 e y entre -4 e 4
a.set(xlim=(-2,2), ylim=(-4,4))


# e finalmente adicionamos o gráfico.
# note que agora vamos até aumentar o valor de x, mas isso
# não tem impacto na janela de visualização, que já definimos
# anteriormente. colocamos 100 pontos para o gráfico ficar suave.
x = np.linspace(-10, 10, 100)
y=f(x)

# agora o plot tem que ser dado com o prefixo a, para usar as configurações
# anteriores.
a.plot(x, y)

[<matplotlib.lines.Line2D at 0x7f79b9047e50>]


t = np.linspace(0, 2*np.pi, 200)
x = np.cos(t)
y = np.sin(2*t)
plt.plot(x, y)

[<matplotlib.lines.Line2D at 0x7f79da22b130>]


import numpy as np
import matplotlib.cm as cm
import matplotlib.pyplot as plt

# o comando abaixo é muito parecido com o linspace que usamos antes.
# a diferença é que aqui, ao invés de definirmos o número de pontos
# que iremos usar para plotar, vamos definir ponto inicial, ponto
# final e o incremento, neste caso delta=0.025. 
# é só uma outra forma de fazer a discretização.
delta = 0.025
x = np.arange(-3.0, 3.0, delta)
y = np.arange(-2.0, 2.0, delta)


# aqui o primeiro truque: transformamos os dois "intervalos" de x e y
# numa malha bidimensional, fazendo uma espécie de produto cartesiano
# dos pontos, e renomeamos a primeira coordenada desses pontos como x
# e a segunda como y. esse é só um artifício técnico.
x, y = np.meshgrid(x, y)

# aqui definimos a funcao/equacao da qual vamos plotar as curvas de nivel.
# note que ela usa as variáveis com letras maiúsculas, pois o que vamos
# armazenar no vetor Z é o valor dessa funcao em cada um dos pontos da
# malha criada.
z = x**2-2*y**2-1

fig, ax = plt.subplots()

# a sintaxe do contor é intervalo X, intervalo Y, funcao e quantidade
# de contornos. abaixo, plotamos 20 curvas de nível.
contornos = ax.contour(x, y, z, 20)

ax.clabel(contornos, inline=True, fontsize=10)
ax.set_title('Curvas de nível')

plt.contour(x, y, z)

<matplotlib.contour.QuadContourSet at 0x7f79c97dd9a0>


import numpy as np
import matplotlib.pyplot as plt

# definindo a função
def f(x, y):
    return np.sqrt(9-x**2+y**2) 

# discretizando x e y
x = np.linspace(-3, 3, 50)
y = np.linspace(-3, 3, 50)

# artifício do produto cartesiano
x, y = np.meshgrid(x, y)

# define z como z=f(x,y)
z = f(x, y)

# indica que iremos fazer um gráfico 3d e nomeia os eixos
ax = plt.axes(projection='3d')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

# finalmente, adiciona a função (no caso, estamos plotando
# em versão wireframe, ou seja, você vai ver a malha sobre a
# superfície.

ax.plot_wireframe(x, y, z)

<mpl_toolkits.mplot3d.art3d.Line3DCollection at 0x7f79b90965e0>


# discretizando as duas variáveis de ângulo
angle = np.linspace(0, 2 * np.pi, 32)
theta, phi = np.meshgrid(angle, angle)

# pegando valores comuns para os raios interno e exteno
r, R = 0.3, 1

# parametrização
x = (R + r * np.cos(phi)) * np.cos(theta)
y = (R + r * np.cos(phi)) * np.sin(theta)
z = r * np.sin(phi)

# criando o ambiente 3d
fig = plt.figure()
ax = plt.axes(projection='3d')

# limitando os eixos
ax.set_xlim3d(-1, 1)
ax.set_ylim3d(-1, 1)
ax.set_zlim3d(-1, 1)

# exibindo a figura
ax.plot_surface(x, y, z, rstride=1, cstride=1)
plt.show()


# criando o ambiente 3d
ax = plt.axes(projection='3d')

# discretizando t
t = np.linspace(0, 2*np.pi, 200)

# definindo as coordenadas da parametrizacao
x = np.cos(t)
y = np.sin(2*t)
z = 2*t

# plotando
ax.plot(x, y, z)

[<mpl_toolkits.mplot3d.art3d.Line3D at 0x7f79da4e3340>]


import sympy as sp
x, y, z = sp.symbols('x y z')

eq1=sp.Eq(x+y-2*z,1)
eq2=sp.Eq(x+y-z,3)
eq3=sp.Eq(3*x+y-z,3)

sp.linsolve((eq1,eq2,eq3),(x,y,z))


sp.solve((eq1,eq2,eq3),(x,y,z))

{x: 0, y: 5, z: 2}


import sympy as sp

# define as variaveis e constantes que vamos usar
x, y, z = sp.symbols('x y z')
a, b, c, d = sp.symbols('a b c d')

# equacao do plano na forma geral
eq = sp.Eq(a*x+b*y+c*z,d)

# definindo os pontos
P=sp.Array([0,0,1])
Q=sp.Array([3,1,0])
R=sp.Array([0,2,2])

# substituindo x,y,z na equacao do plano pelos pontos
eq1=eq.subs([(x,P[0]),(y,P[1]),(z,P[2])])
eq2=eq.subs([(x,Q[0]),(y,Q[1]),(z,Q[2])])
eq3=eq.subs([(x,R[0]),(y,R[1]),(z,R[2])])

# as equacoes acima resultam num sistema nas variaveis a, b, c, d.
# resolvendo o sistema:
sol=sp.linsolve((eq1,eq2,eq3),(a,b,c,d))

#armazenando as solucoes
(A,B,C,D)=tuple(*sol)

# finalmente, exibindo a equação do plano, substituindo a,b,c,d pelos
# valores que encontramos
eq.subs([(a,A),(b,B),(c,C),(d,D)])


# importando a biblioteca para calculos simbolicos
import sympy as sp

# vamos usar a biblioteca NumPy para calcular o produto vetorial e tambem
# o produto interno
import numpy as np

# define as variaveis e constantes que vamos usar
x, y, z = sp.symbols('x y z')
a, b, c, d = sp.symbols('a b c d')

# equacao do plano na forma geral
eq = sp.Eq(a*x+b*y+c*z,d)

# definindo os pontos
P=sp.Array([0,0,1])
Q=sp.Array([3,1,0])
R=sp.Array([0,2,2])

# criando os vetores
u=Q-P
v=R-P

# obtendo o vetor normal
n=np.cross(u,v)
(a,b,c)=n

# termo nao-homogeneo
d=np.dot(n,P)

# exibindo a equação do plano
sp.Eq(a*x+b*y+c*z,d)


import sympy as sp

A=sp.Matrix([[1,2],[-3,1]])
B=sp.Matrix([[2],[4]])
C=sp.Matrix([[2,5],[4,10]])

A

B

C


sp.eye(2)


sp.zeros(2,2)


sp.ones(2,2)

A*C

A*B

B*A

---------------------------------------------------------------------------
ShapeError                                Traceback (most recent call last)
/var/folders/90/lb7fswt93477cpj1pl8z494m0000gn/T/ipykernel_41171/3451148411.py in <module>
----> 1 B*A

~/opt/anaconda3/lib/python3.9/site-packages/sympy/core/decorators.py in binary_op_wrapper(self, other)
    134                     if f is not None:
    135                         return f(self)
--> 136             return func(self, other)
    137         return binary_op_wrapper
    138     return priority_decorator

~/opt/anaconda3/lib/python3.9/site-packages/sympy/matrices/common.py in __mul__(self, other)
   2761         """
   2762 
-> 2763         return self.multiply(other)
   2764 
   2765     def multiply(self, other, dotprodsimp=None):

~/opt/anaconda3/lib/python3.9/site-packages/sympy/matrices/common.py in multiply(self, other, dotprodsimp)
   2783              getattr(other, 'is_MatrixLike', True))):
   2784             if self.shape[1] != other.shape[0]:
-> 2785                 raise ShapeError("Matrix size mismatch: %s * %s." % (
   2786                     self.shape, other.shape))
   2787 

ShapeError: Matrix size mismatch: (2, 1) * (2, 2).


A*C-C*A


A**(-1)


A.det()


A.eigenvals()

{1 - sqrt(6)*I: 1, 1 + sqrt(6)*I: 1}


A.eigenvects()

[(1 - sqrt(6)*I,
  1,
  [Matrix([
   [sqrt(6)*I/3],
   [          1]])]),
 (1 + sqrt(6)*I,
  1,
  [Matrix([
   [-sqrt(6)*I/3],
   [           1]])])]


sp.init_printing(use_unicode=True)
A.eigenvects()


M=sp.Matrix([[3,0,0],[0,3,0],[0,1,1]])
M.eigenvects()


A.nullspace()


C.nullspace()


C.eigenvects()


A.columnspace()


C.columnspace()


A.diagonalize()


Q,D = A.diagonalize()


sp.simplify(Q*D*(Q**(-1)))


a=1
b=2
c=-3
d=1
e=2
f=-4
x, y = sp.symbols('x y')
sp.Eq(a*x**2+b*x*y+c*y**2+d*x+e*y+f,0)


import numpy as np
import matplotlib.cm as cm
import matplotlib.pyplot as plt
delta = 0.01
x = np.arange(-4.0, 4.0, delta)
y = np.arange(-4.0, 4.0, delta)
x, y = np.meshgrid(x, y)
z = a*x**2+b*x*y+c*y**2+d*x+e*y+f

fig = plt.figure()
ax = fig.add_subplot()
ax.spines['left'].set_position('center')
ax.spines['bottom'].set_position('center')
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')

plt.contour(x, y, z,0)

<matplotlib.contour.QuadContourSet at 0x7f79b9335280>


A=sp.Matrix([[a,sp.Rational(b, 2)],[sp.Rational(b, 2),c]])


Q,D = A.diagonalize()

Q,D


n1=sp.simplify(sp.sqrt(Q[0]**2+Q[2]**2))
n2=sp.simplify(sp.sqrt(Q[1]**2+Q[3]**2))
Q2=sp.simplify(sp.Matrix([[Q[0]/n1,Q[1]/n2],[Q[2]/n1,Q[3]/n2]]))
Q=Q2
Q


sp.simplify(Q**(-1)-Q.transpose())


sp.simplify(Q**(-1)*A*Q)


# aplicando a mudança de coordenadas; z2 será a nova equacao

x, y = sp.symbols('x y')
z2=sp.simplify(
    (sp.Matrix([[x,y]])*(Q**(-1)*A*Q)*sp.Matrix([[x],[y]]))+(sp.Matrix([[d,e]]))*Q*(sp.Matrix([[x],[y]]))+sp.Matrix([[f]])
)[0]


sp.init_printing(use_unicode=False)


# exibindo a nova equacao, e como vamos usá-la para plotar, passando tudo
# pra ponto flutuante. você vai precisar copiar o resultado desse comando e
# colar na próxima linha de código, no z = ....
z3=z2.evalf()
repr(z3)

'1.23606797749979*x**2 + 1.43275483056186*x - 3.23606797749979*y**2 + 1.71674505838791*y - 4.0'


import numpy as np
import matplotlib.cm as cm
import matplotlib.pyplot as plt
delta = 0.01
x = np.arange(-4.0, 4.0, delta)
y = np.arange(-4.0, 4.0, delta)
x, y = np.meshgrid(x, y)
z3 = 1.23606797749979*x**2 + 1.43275483056186*x - 3.23606797749979*y**2 + 1.71674505838791*y - 4.0

fig = plt.figure()
ax = fig.add_subplot()
ax.spines['left'].set_position('center')
ax.spines['bottom'].set_position('center')
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')

plt.contour(x, y, z3,0)

<matplotlib.contour.QuadContourSet at 0x7f79a850eaf0>


lam = D[0]
mu = D[3]


x, y = sp.symbols('x y')
termolinear=sp.simplify((sp.Matrix([[d,e]]))*Q*(sp.Matrix([[x],[y]])))[0]
alfa=sp.diff(termolinear,x)
beta=sp.diff(termolinear,y)


# definindo quem é f-til
ftil=f-(alfa**2)/(4*lam)-beta**2/(4*mu)

# resetando as variáveis e escrevendo a equacao final
x, y = sp.symbols('x y')
print((lam*x**2 + mu*y**2 +ftil).evalf())

# plotando
import numpy as np
import matplotlib.cm as cm
import matplotlib.pyplot as plt
delta = 0.01
x = np.arange(-4.0, 4.0, delta)
y = np.arange(-4.0, 4.0, delta)
x, y = np.meshgrid(x, y)
z4 = 1.23606797749979*x**2 - 3.23606797749979*y**2 - 4.1875

fig = plt.figure()
ax = fig.add_subplot()
ax.spines['left'].set_position('center')
ax.spines['bottom'].set_position('center')
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')

plt.contour(x, y, z4,0)

1.23606797749979*x**2 - 3.23606797749979*y**2 - 4.1875

<matplotlib.contour.QuadContourSet at 0x7f79a84cd580>


# plotando
import numpy as np
import matplotlib.cm as cm
import matplotlib.pyplot as plt
delta = 0.01
x = np.arange(-4.0, 4.0, delta)
y = np.arange(-4.0, 4.0, delta)
x, y = np.meshgrid(x, y)

z = a*x**2+b*x*y+c*y**2+d*x+e*y+f
z3 = 1.23606797749979*x**2 + 1.43275483056186*x - 3.23606797749979*y**2 + 1.71674505838791*y - 4.0
z4 = 1.23606797749979*x**2 - 3.23606797749979*y**2 - 4.1875

fig = plt.figure()
ax = fig.add_subplot()
ax.spines['left'].set_position('center')
ax.spines['bottom'].set_position('center')
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')

plt.contour(x, y, z,0,colors=['red'])
plt.contour(x, y, z3,0,colors=['blue'])
plt.contour(x, y, z4,0,colors=['black'])
plt.show()


x, y = sp.symbols('x y')
sp.Eq(lam*x**2+mu*y**2+ftil,0)


import sympy as sp
x = sp.symbols('x')
sp.limit(x**2,x,2)


sp.limit(sp.sin(x)/x,x,0)


import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-0.5, 0.5, 100)
y = np.sin(x)
plt.plot(x, y)
plt.plot(x, x)

[<matplotlib.lines.Line2D at 0x7f79dab1fdc0>]


import sympy as sp
x = sp.symbols('x')
sp.limit(1/(x-2),x,2)


sp.limit(1/x,x,sp.oo)


sp.limit(1/x, x, 0, '+')


sp.limit(1/x, x, 0, '-')


import sympy as sp
x, h = sp.symbols('x h', real=True)
f=x**2
fh=f.subs(x,x+h)
sp.limit( (fh-f)/h , h, 0)


import sympy as sp
x = sp.symbols('x', real=True)
g=sp.Piecewise((x**2+1, x<1),(3*x, x>1))


sp.limit(g,x,1,'-')


sp.limit(g,x,1,'+')


import sympy as sp
x = sp.symbols('x')
f=x**2+x
sp.diff(f,x)


g=sp.exp(x**3+x)+1/x
sp.diff(g,x)


y = sp.symbols('y')
h=sp.cos(x*y)+(3*x+y)/(y**2+1)


# derivada em x
sp.diff(h,x)


# derivada em y
sp.diff(h,x)


# derivada mista, em x depois em y
sp.diff(h,x,y)


sp.diff(sp.diff(h,x),y)


sp.diff(h,x,x,x)


sp.diff(h,x).subs(x,0).subs(y,0)


sp.diff(g,x).subs(x,0).subs(y,0)


import sympy as sp
x = sp.symbols('x')
f=x*sp.cos(x)


# calculando a expansao formal de f em serie de potencias
# com respeito aa variavel x, no ponto x=0, e truncando
# em ordem 10 usando o fps
fn=sp.fps(f,x,0).truncate(10)
display(fn)


# o comando series tambem funciona: a sintaxe é 
# series(f(x),x,x0,k), onde x0 
fm=sp.series(f,x,0,20).as_expr()
display(fm)


import sympy as sp
x = sp.symbols('x')
g = sp.exp(x)
gk = sp.series(g,x,0,10).removeO()
gk.subs(x,2)


20947//2835


sp.series(g,x,0,20).as_expr().removeO()


import sympy as sp
x = sp.symbols('x')
f=x**2
sp.integrate(f,x)


import sympy as sp
x = sp.symbols('x')
g=x*sp.exp(x)
sp.integrate(g,x)


import sympy as sp
x = sp.symbols('x')
h=sp.exp(x)*sp.cos(x)
sp.integrate(h,x)


p=x**3+x
sp.integrate(p,(x,0,1))


# vamos tentar enganar o python?
q=1/(x-2)
sp.integrate(q,(x,0,3))


sp.integrate(sp.exp(-x**2), (x, 0, sp.oo))


sp.integrate(1/(x**(6)), (x, 1,sp.oo))


x, y, a, b, c, d = sp.symbols("x y a b c d")
f = x*y
sp.integrate(f, (y, a, b), (x, c, d))


x, y, a, b, c, d = sp.symbols("x y a b c d")
f = x**2+y
sp.integrate(f, (y, 0, x+1), (x, 0, 1))


import numpy as np
import matplotlib.pyplot as plt

# originamente, a função foi definida com a definicao lambda,
# que de fato é mais simples nesse caso - caso queira saber
# mais sobre ela, leia aqui:
# https://stackabuse.com/lambda-functions-in-python/
# f = lambda x : x**2

# definindo a função
def f(x):
    return x**2

# intervalo
a = 0; b = 10;

# numero de retangulos
N = 10

n = 10 # Use n*N+1 points to plot the function smoothly

# discretizando variáveis
x = np.linspace(a,b,N+1)
y = f(x)
X = np.linspace(a,b,n*N+1)
Y = f(X)

# iniciando plot
plt.figure(figsize=(15,5))
plt.subplot(1,3,1)
plt.plot(X,Y,'b')
x_left = x[:-1]
y_left = y[:-1]
plt.plot(x_left,y_left,'b.',markersize=10)
plt.bar(x_left,y_left,width=(b-a)/N,alpha=0.2,align='edge',edgecolor='b')
plt.title('Soma de Riemann com base na esquerda'.format(N))

plt.show()


import sympy as sp

# definindo a variavel independente
t = sp.symbols('t')

# definindo símbolos para as funcoes que estarão envolvidas
# nas equações diferenciais.
# note que precisamos já declarar esses símbolos como sendo
# da classe "função" com a adição cls=sp.Function
f = sp.symbols('f', cls=sp.Function)

# defina a equacao diferencial. para representar f'(t)
# deve-se usar f(t).diff(t).
# para indicar f''(t) usamos f(t).diff(t,t).
# a equacao sera na forma F=0. defina eq=F.
eq=f(t).diff(t) + f(t)

# resolvendo a equacao
sp.dsolve(eq, f(t))


import sympy as sp
t = sp.symbols('t')
g = sp.symbols('g', cls=sp.Function)
eq=g(t).diff(t,t) -2*g(t).diff(t) +g(t)
sp.dsolve(eq, g(t))


import sympy as sp
t = sp.symbols('t')
q = sp.symbols('q', cls=sp.Function)
eq=q(t).diff(t,t) -6*q(t).diff(t) +2*q(t)-t*sp.cos(t)
sp.dsolve(eq, q(t))


import sympy as sp
t = sp.symbols('t')
q = sp.symbols('q', cls=sp.Function)
eq=q(t).diff(t,t) -6*q(t).diff(t) +2*q(t)-t*sp.cos(t)

# encontrando a solucao e armazenando no nome "solucao"
solucao=sp.dsolve(eq, q(t), ics={q(0): 1, q(t).diff(t).subs(t, 0): 1})

sp.plot(solucao.rhs,(t,-2,2),ylim=[-5,5])

<sympy.plotting.plot.Plot at 0x7f79b994c370>


import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def lorenz(state, t, sigma, beta, rho):
    x, y, z = state
    dx = sigma * (y - x)
    dy = x * (rho - z) - y
    dz = x * y - beta * z
    return [dx, dy, dz]

# valores dos parametros
sigma = 10.0
beta = 8.0 / 3.0
rho = 28.0

# parametros
p = (sigma, beta, rho)

# condicao inicial
ic = [1.0, 1.0, 1.0]

# tempo inicial, tempo final, tamanho da partição.
# se você colocar uma partição muito grande, então serão
# usados menos pontos para avaliar a função, e com isso
# terá a impressão de que são vários segmentos de reta.
# experimente um pouco com isso.
t = np.arange(0.0, 60.0, 0.001)

# integrando
result = odeint(lorenz, ic, t, p)

# plotando
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot(result[:, 0], result[:, 1], result[:, 2])

fig2 = plt.figure()
plt.plot(t,result[:, 0])
plt.plot(t,result[:, 1])
plt.plot(t,result[:, 2])

[<matplotlib.lines.Line2D at 0x7f79dab948e0>]


import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# define o sistema de van der pol
def sistema(variaveis, t):
    x, y = variaveis
    dx = 3*(1-y**2)*x-y
    dy = x
    return [dx, dy]

# define a condicao inicial, onde 
# ic=[ x(0), y(0)]
ic = [0.1,0.2]

# define o tempo - será usado pelo integrador e também
# pelo plot
t = np.arange(0, 40, 0.0001)

# integrando
result = odeint(sistema, ic, t)

# primeiro plot
fig = plt.figure()
ax = plt.axes()
# plotando a curva-solucao
ax.plot(result[:, 0], result[:, 1])

# agora plotando x(t) e y(t) separadamente num segundo plot
fig2=plt.figure()
plt.plot(t,result[:, 0])
plt.plot(t,result[:, 1])

[<matplotlib.lines.Line2D at 0x7f79db019c10>]


from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# note que agora que queremos usar o solve_ivp precisamos
# colocar o t como primeira variável, ao contrário de quando
# queríamos usar o odeint.
def sistema(t, variaveis):
    x, y = variaveis
    dx = 3*(1-y**2)*x-y
    dy = x
    return [dx, dy]

# condicao inicial
ic=[0.1,0.2]

# basicamente a sintaxe eh solve_ivp(sistema [t0,t1], ic). adicionamos
# a opcao dense_output para ficar mais facil recuperar as solucoes para
# que sejam plotadas.

solucao = solve_ivp(sistema, [-40, 40], ic,dense_output=True)

# agora discretizamos o tempo para plotar. é preciso ser compatível com
# o intervalo do tempo que foi passado para o comando solve_ivp
t = np.arange(-40, 40, 0.0001)

# a propriedade sol(t) da solucao carrega em suas colunas os valores
# de x(t) e y(t). para usar no plot, precisamos que isso esteja nas linhas,
# por isso usamos o .T no final, para calcular a transposta. isso vai plotar
# as solucoes x(t), y(t)
plt.plot(t,solucao.sol(t).T)

# agora criamos outra janela gráfica e plotamos a curva parametrizada.
# usamos [:,0] para acessar a primeira linha da matriz solucao.sol(t)
# e [:,1] para acessar a segunda linha. isso vai produzir um plot de
# curva parametrizada, como já sabemos fazer.
fig=plt.figure()
plt.plot(solucao.sol(t).T[:, 0],solucao.sol(t).T[:, 1])

[<matplotlib.lines.Line2D at 0x7f79db1dcca0>]


import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
from scipy.integrate import solve_ivp

fig = plt.figure(num=1)
ax=fig.add_subplot(111)

def sistema(t, variaveis):
    x, y = variaveis
    dx = (1-y**2)*x-y
    dy = x
    return [dx, dy]

# trajetórias em tempo positivo
for P in range(-5,5,1):
  for Q in range(-5,5,1):
    solucao = solve_ivp(sistema, [0,10], [P,Q],dense_output=True)
    t = np.linspace(0, 10, 500)
    plt.plot(solucao.sol(t).T[:, 0],solucao.sol(t).T[:, 1])

# trajetórias em tempo negativo
for P in range(-5,5,1):
  for Q in range(-5,5,1):
    solucao = solve_ivp(sistema, [0,-10], [P,Q],dense_output=True)
    t = np.linspace(0, -10, 500)
    plt.plot(solucao.sol(t).T[:, 0],solucao.sol(t).T[:, 1])

# limita a janela de visualização e mostra o plot
plt.xlim(-5,5)
plt.ylim(-5,5)
plt.show()


import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
from scipy.integrate import solve_ivp

fig = plt.figure(num=1)
ax=fig.add_subplot(111)

# um sistema que só tem órbitas periódicas
def sistema(t, variaveis):
    x, y = variaveis
    dx = -y
    dy = x
    return [dx, dy]

# trajetórias em tempo positivo
for P in range(-5,5,1):
  for Q in range(-5,5,1):
    solucao = solve_ivp(sistema, [0,10], [P,Q],dense_output=True)
    t = np.linspace(0, 10, 500)
    plt.plot(solucao.sol(t).T[:, 0],solucao.sol(t).T[:, 1])

# trajetórias em tempo negativo
for P in range(-5,5,1):
  for Q in range(-5,5,1):
    solucao = solve_ivp(sistema, [0,-10], [P,Q],dense_output=True)
    t = np.linspace(0, -10, 500)
    plt.plot(solucao.sol(t).T[:, 0],solucao.sol(t).T[:, 1])

# limita a janela de visualização e mostra o plot
plt.xlim(-5,5)
plt.ylim(-5,5)
plt.show()


import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

fig = plt.figure(num=1)
ax=fig.add_subplot(111)

# um sistema que só tem órbitas periódicas
def sistema(t, variaveis):
    x, y = variaveis
    dx = y
    dy = x
    return [dx, dy]

# trajetórias em tempo positivo
for P in range(-5,5,1):
  for Q in range(-5,5,1):
    solucao = solve_ivp(sistema, [0,10], [P,Q],dense_output=True)
    t = np.linspace(0, 10, 500)
    plt.plot(solucao.sol(t).T[:, 0],solucao.sol(t).T[:, 1])

# trajetórias em tempo negativo
for P in range(-5,5,1):
  for Q in range(-5,5,1):
    solucao = solve_ivp(sistema, [0,-10], [P,Q],dense_output=True)
    t = np.linspace(0, -10, 500)
    plt.plot(solucao.sol(t).T[:, 0],solucao.sol(t).T[:, 1])

# limita a janela de visualização e mostra o plot
plt.xlim(-5,5)
plt.ylim(-5,5)
plt.show()


import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

fig = plt.figure(num=1)
ax=fig.add_subplot(111)

# um sistema hamiltoniano com órbitas periódicas e uma órbita homoclínica
def sistema(t, variaveis):
    x, y = variaveis
    dx = -y
    dy = -x-x**2
    return [dx, dy]

# trajetórias em tempo positivo
for P in np.linspace(-2, 2, 20):
  for Q in np.linspace(-2, 2, 20):
    solucao = solve_ivp(sistema, [0,10], [P,Q],dense_output=True)
    t = np.linspace(0, 10, 500)
    plt.plot(solucao.sol(t).T[:, 0],solucao.sol(t).T[:, 1])

# trajetórias em tempo negativo
for P in np.linspace(-2, 2, 20):
  for Q in np.linspace(-2, 2, 20):
    solucao = solve_ivp(sistema, [0,-10], [P,Q],dense_output=True)
    t = np.linspace(0, -10, 500)
    plt.plot(solucao.sol(t).T[:, 0],solucao.sol(t).T[:, 1])

# limita a janela de visualização e mostra o plot
plt.xlim(-2,2)
plt.ylim(-2,2)
plt.show()


import matplotlib.pyplot as plt
import numpy as np

# criamos uma malha em R2, com x e y indo de -5 a 5, contendo
# um total de 10 pontos em cada coordenada
x,y = np.meshgrid( np.linspace(-5,5,10),np.linspace(-5,5,10) )

# calculamos o campo vetorial X(x,y)=(y,x) na malha
u = y
v = x

#N = np.sqrt(u**2+v**2)
#U2, V2 = u/N, v/N

# plotando
plt.quiver( x,y,u, v)

<matplotlib.quiver.Quiver at 0x7f79ba40d880>


import matplotlib.pyplot as plt
import numpy as np

# criamos uma malha em R2, com x e y indo de -5 a 5, contendo
# um total de 10 pontos em cada coordenada
x,y = np.meshgrid( np.linspace(-5,5,10),np.linspace(-5,5,10) )

# calculamos o campo vetorial X(x,y)=(y,x) na malha
u = y
v = x

# normalizando o campo
N = np.sqrt(u**2+v**2)
U, V = u/N, v/N

# plotando
plt.quiver( x,y,U, V)

<matplotlib.quiver.Quiver at 0x7f79cb89b940>


import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

fig = plt.figure(num=1)
ax=fig.add_subplot(111)

# um sistema que só tem órbitas periódicas
def sistema(t, variaveis):
    x, y = variaveis
    dx = y
    dy = x
    return [dx, dy]

# calculando e plotando trajetórias em tempo positivo
for P in range(-5,5,1):
  for Q in range(-5,5,1):
    solucao = solve_ivp(sistema, [0,10], [P,Q],dense_output=True)
    t = np.linspace(0, 10, 500)
    plt.plot(solucao.sol(t).T[:, 0],solucao.sol(t).T[:, 1])

# calculando e plotando trajetórias em tempo negativo
for P in range(-5,5,1):
  for Q in range(-5,5,1):
    solucao = solve_ivp(sistema, [0,-10], [P,Q],dense_output=True)
    t = np.linspace(0, -10, 500)
    plt.plot(solucao.sol(t).T[:, 0],solucao.sol(t).T[:, 1])

# grid para os vetores
x,y = np.meshgrid( np.linspace(-5,5,20),np.linspace(-5,5,20) )

# calculando o campo vetorial X(x,y)=(y,x) na malha
u = y
v = x

# normalizando o campo (opcional)
N = np.sqrt(u**2+v**2)
u, v = u/N, v/N

# plotando o campo de direcoes
plt.quiver( x,y,u, v)

# limita a janela de visualização e mostra o plot
plt.xlim(-5,5)
plt.ylim(-5,5)
plt.show()


%matplotlib inline
%matplotlib tk

import numpy as np
from numpy import pi
import matplotlib.pyplot as plt
import matplotlib.animation as animation


fig = plt.figure()
fig.set_dpi(100)
ax1 = fig.add_subplot(1,1,1)

#Diffusion constant
k = 2

#Scaling factor (for visualisation purposes)
scale = 5

#Length of the rod (0,L) on the x axis
L = pi

#Initial contitions u(0,t) = u(L,t) = 0. Temperature at x=0 and x=L is fixed
x0 = np.linspace(0,L+1,10000)
t0 = 0
temp0 = 5  #Temperature of the rod at rest (before heating)

#Increment
dt = 0.01

# Solução da equação do calor
def u(x,t):
    return temp0 + scale*np.exp(-k*t)*np.sin(x)

#Gradient of u
def grad_u(x,t):
                               #du/dx              #du/dt
    return scale*np.array([np.exp(-k*t)*np.cos(x),-k*np.exp(-k*t)*np.sin(x)])

a = []
t = []

for i in range(500):
    value = u(x0,t0) + grad_u(x0,t0)[1]*dt
    t.append(t0)
    t0 = t0 + dt
    a.append(value)

k = 0
def animate(i):         #The plot shows the temperature evolving with time
    global k            #at each point x in the rod
    x = a[k]            #The ends of the rod are kept at temperature temp0
    k += 1              #The rod is heated in one spot, then it cools down
    ax1.clear()
    plt.plot(x0,x,color='red',label='Temperatura em cada ponto x')
    plt.plot(0,0,color='red',label='Tempo decorrido '+str(round(t[k],2)))
    plt.grid(True)
    plt.ylim([temp0-2,2.5*scale])
    plt.xlim([0,L])
    plt.title('Evolução da equação do calor')
    plt.legend()
 
anim = animation.FuncAnimation(fig,animate,frames=360,interval=10)
plt.show()


%matplotlib inline
%matplotlib tk

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation

print("2D heat equation solver")

plate_length = 50
max_iter_time = 750

alpha = 2
delta_x = 1

delta_t = (delta_x ** 2)/(4 * alpha)
gamma = (alpha * delta_t) / (delta_x ** 2)

# Initialize solution: the grid of u(k, i, j)
u = np.empty((max_iter_time, plate_length, plate_length))

# Initial condition everywhere inside the grid
u_initial = 0

# Boundary conditions
u_top = 100.0
u_left = 0.0
u_bottom = 0.0
u_right = 0.0

# Set the initial condition
u.fill(u_initial)

# Set the boundary conditions
u[:, (plate_length-1):, :] = u_top
u[:, :, :1] = u_left
u[:, :1, 1:] = u_bottom
u[:, :, (plate_length-1):] = u_right

def calculate(u):
    for k in range(0, max_iter_time-1, 1):
        for i in range(1, plate_length-1, delta_x):
            for j in range(1, plate_length-1, delta_x):
                u[k + 1, i, j] = gamma * (u[k][i+1][j] + u[k][i-1][j] + u[k][i][j+1] + u[k][i][j-1] - 4*u[k][i][j]) + u[k][i][j]

    return u

def plotheatmap(u_k, k):
    plt.clf()

    plt.title(f"Gráfico da temperatura no instante t = {k*delta_t:.3f} unit time")
    plt.xlabel("x")
    plt.ylabel("y")

    # This is to plot u_k (u at time-step k)
    plt.pcolormesh(u_k, cmap=plt.cm.jet, vmin=0, vmax=100)
    plt.colorbar()

    return plt

# Do the calculation here
u = calculate(u)

def animate(k):
    plotheatmap(u[k], k)

anim = animation.FuncAnimation(plt.figure(), animate, interval=1, frames=max_iter_time, repeat=False)

2D heat equation solver

Tutorial rápido de Python para Matemáticos¶

O que tem aqui?¶

Referências bibliográficas¶

Mensagem legal e copyright¶

Python é uma boa calculadora!¶

Sem preguiça, vamos programar um pouco!¶

Bônus: programando nossa própria raiz quadrada¶

Resolvendo equações¶

Encontrar soluções de equações? Ideias simples podem funcionar.¶

Gráficos¶

Resolvendo sistemas lineares¶

Um exercício legal para testar nossos conhecimentos¶

Matrizes¶

Um exercício interessante: reconhecimento de cônicas¶

Limites¶

Derivadas¶

Expansão em séries¶

Integrais¶

Somas de Riemann¶

Equações diferenciais¶

Equações lineares de primeira e segunda ordens¶

Sistemas de equações diferenciais¶

Retratos de fase e campos de direções¶

EDPs¶