Calcul symbolique avec SymPy

pip install sympy

import sympy as sp

# Déclarer des symboles
x, y, z = sp.symbols('x y z')

# Créer des expressions
expr = x**2 + 2*x + 1
print(expr)  # x**2 + 2*x + 1

# Factoriser
print(sp.factor(expr))  # (x + 1)**2

# Développer
print(sp.expand((x + 1)**3))  # x**3 + 3*x**2 + 3*x + 1

# Simplifier
expr2 = sp.sin(x)**2 + sp.cos(x)**2
print(sp.simplify(expr2))  # 1

import sympy as sp

x = sp.Symbol('x')

# Dérivée de x³
print(sp.diff(x**3, x))  # 3*x**2

# Dérivée de sin(x²)
print(sp.diff(sp.sin(x**2), x))  # 2*x*cos(x**2)

# Dérivée de exp(x)*ln(x)
print(sp.diff(sp.exp(x) * sp.ln(x), x))  # exp(x)*log(x) + exp(x)/x

# Dérivée seconde
print(sp.diff(x**4, x, 2))  # 12*x**2

# Dérivée n-ième
n = sp.Symbol('n', positive=True, integer=True)
print(sp.diff(sp.exp(x), x, 3))  # exp(x)

import sympy as sp

x = sp.Symbol('x')

# Notre fonction : sqrt(1 + cos²(x))
f = sp.sqrt(1 + sp.cos(x)**2)

# Calculer f''(x)
f_seconde = sp.diff(f, x, 2)
print("f''(x) =", f_seconde)

# Simplifier
f_seconde_simp = sp.simplify(f_seconde)
print("f''(x) simplifié =", f_seconde_simp)

# Trouver le maximum sur [0, π]
# On cherche les points critiques de |f''|
f_tierce = sp.diff(f_seconde, x)
points_critiques = sp.solve(f_tierce, x)
print("Points critiques de f'' :", points_critiques)

import sympy as sp

x = sp.Symbol('x')

# Primitive de x²
print(sp.integrate(x**2, x))  # x**3/3

# Primitive de 1/(1+x²)
print(sp.integrate(1/(1 + x**2), x))  # atan(x)

# Primitive de exp(-x²) — n'existe pas en forme close !
resultat = sp.integrate(sp.exp(-x**2), x)
print(resultat)  # sqrt(pi)*erf(x)/2

# SymPy retourne la fonction spéciale erf (fonction d'erreur)

import sympy as sp

x = sp.Symbol('x')

# Intégrale de 0 à 1 de x²
print(sp.integrate(x**2, (x, 0, 1)))  # 1/3

# Intégrale de 0 à π de sin(x)
print(sp.integrate(sp.sin(x), (x, 0, sp.pi)))  # 2

# Intégrale de -∞ à +∞ de exp(-x²) — la célèbre gaussienne
print(sp.integrate(sp.exp(-x**2), (x, -sp.oo, sp.oo)))  # sqrt(pi)

# Notre intégrale elliptique
integrale_elliptique = sp.integrate(sp.sqrt(1 + sp.cos(x)**2), (x, 0, sp.pi))
print(integrale_elliptique)  # 2*sqrt(2)*elliptic_e(1/2)

# SymPy reconnaît l'intégrale elliptique et retourne le résultat exact !

import sympy as sp

x = sp.Symbol('x')

# Équation du second degré
print(sp.solve(x**2 - 5*x + 6, x))  # [2, 3]

# Équation du troisième degré
print(sp.solve(x**3 - 6*x**2 + 11*x - 6, x))  # [1, 2, 3]

# Équation avec paramètres
a, b, c = sp.symbols('a b c')
solutions = sp.solve(a*x**2 + b*x + c, x)
print(solutions)  # [(-b - sqrt(b**2 - 4*a*c))/(2*a), (-b + sqrt(b**2 - 4*a*c))/(2*a)]

import sympy as sp

x = sp.Symbol('x')

# cos(x) = x (équation de point fixe)
solutions = sp.solve(sp.cos(x) - x, x)
print(solutions)  # [] — SymPy ne trouve pas de solution analytique

# C'est normal ! Cette équation n'a pas de solution en forme close
# Il faut une méthode numérique (Newton, dichotomie...)

import sympy as sp

x = sp.Symbol('x')

# Développement de exp(x) autour de 0, ordre 5
print(sp.series(sp.exp(x), x, 0, 6))
# 1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + O(x**6)

# Développement de sin(x)
print(sp.series(sp.sin(x), x, 0, 8))
# x - x**3/6 + x**5/120 - x**7/5040 + O(x**8)

# Développement de 1/(1-x)
print(sp.series(1/(1-x), x, 0, 6))
# 1 + x + x**2 + x**3 + x**4 + x**5 + O(x**6)

# Développement autour d'un autre point
print(sp.series(sp.ln(x), x, 1, 4))  # Autour de x=1
# (x - 1) - (x - 1)**2/2 + (x - 1)**3/3 + O((x - 1)**4, (x, 1))

import sympy as sp
import numpy as np

x = sp.Symbol('x')

# Approximation de sin(x) par son développement
def erreur_taylor_sin(x_val, n_termes):
  """Compare sin(x) exact avec son développement de Taylor."""
  # Développement de Taylor
  taylor = sp.series(sp.sin(x), x, 0, 2*n_termes).removeO()

  # Évaluer numériquement
  exact = np.sin(x_val)
  approx = float(taylor.subs(x, x_val))

  return abs(exact - approx)

# Test pour x = 1
for n in [1, 2, 3, 4, 5]:
  err = erreur_taylor_sin(1.0, n)
  print(f"{n} terme(s) : erreur = {err:.2e}")

import sympy as sp

# Matrice symbolique
a, b, c, d = sp.symbols('a b c d')
M = sp.Matrix([[a, b], [c, d]])

# Déterminant
print("det(M) =", M.det())  # a*d - b*c

# Inverse
print("M^(-1) =", M.inv())

# Valeurs propres
print("Valeurs propres :", M.eigenvals())

# Matrice numérique
A = sp.Matrix([[1, 2], [3, 4]])
print("det(A) =", A.det())  # -2
print("A^(-1) =", A.inv())  # Matrix([[-2, 1], [3/2, -1/2]])

import sympy as sp

t = sp.Symbol('t')
y = sp.Function('y')

# y' = y (solution : exp(t))
eq1 = sp.Eq(y(t).diff(t), y(t))
print(sp.dsolve(eq1))  # Eq(y(t), C1*exp(t))

# y'' + y = 0 (oscillateur harmonique)
eq2 = sp.Eq(y(t).diff(t, 2) + y(t), 0)
print(sp.dsolve(eq2))  # Eq(y(t), C1*sin(t) + C2*cos(t))

# y' = y² (équation de Riccati simple)
eq3 = sp.Eq(y(t).diff(t), y(t)**2)
print(sp.dsolve(eq3))  # Eq(y(t), -1/(C1 + t))

# Avec condition initiale
eq4 = sp.Eq(y(t).diff(t), -2*y(t))
sol = sp.dsolve(eq4, ics={y(0): 5})
print(sol)  # Eq(y(t), 5*exp(-2*t))

import sympy as sp
import numpy as np

x = sp.Symbol('x')

# Expression symbolique
expr = sp.sin(x**2) + sp.exp(-x)

# Convertir en fonction NumPy
f_numpy = sp.lambdify(x, expr, 'numpy')

# Utiliser avec des tableaux NumPy
x_vals = np.linspace(0, 2, 100)
y_vals = f_numpy(x_vals)

# Maintenant on peut utiliser l'intégration numérique
from scipy import integrate
result, error = integrate.quad(f_numpy, 0, 2)
print(f"Intégrale numérique : {result:.10f}")