La résolution des équations non linéaires en python

Il y a deux façons de le faire.

1. Utiliser un non-solveur linéaire
2. Linéariser le problème et de le résoudre dans la méthode des moindres carrés sens

Installation

Donc, si je comprends bien votre question, vous savez F, a, b, et c à 4 points de vue différents, et que vous souhaitez inverser pour les paramètres du modèle X, Y et Z. Nous avons 3 inconnues et 4 les données observées points, donc le problème est surdéterminé. Donc, nous allons être dans la résolution de la méthode des moindres carrés sens.

Il est plus courant d'utiliser l'inverse de la terminologie dans ce cas, nous allons donc retourner votre équation de autour de. Au lieu de:

F_i = X^2 + a_i Y^2 + b_i X Y cosZ + c_i X Y sinZ

Écrivons:

F_i = a^2 + X_i b^2 + Y_i a b cos(c) + Z_i a b sin(c)

Où nous savons F, X, Y, et Z à 4 points différents (par exemple F_0, F_1, ... F_i).

Nous sommes juste en changeant les noms des variables, pas l'équation elle-même. (C'est plus pour mon aise de penser que quoi que ce soit d'autre.)

Linéaire De La Solution

Il est possible de linéariser cette équation. Vous pouvez facilement résoudre pour a^2, b^2, a b cos(c), et a b sin(c). Pour rendre cela un peu plus facile, nous allons reclasser les choses encore une fois:

d = a^2
e = b^2
f = a b cos(c)
g = a b sin(c)

Maintenant l'équation est beaucoup plus simple: F_i = d + e X_i + f Y_i + g Z_i. Il est facile de faire un linéaire des moindres carrés inversion pour d, e, f, et g. Nous pouvons alors obtenir a, b, et c à partir de:

a = sqrt(d)
b = sqrt(e)
c = arctan(g/f)

Bien, nous allons écrire ça sous forme de matrice. Nous allons traduire 4 observations de (le code, nous allons écrire vont prendre n'importe quel nombre d'observations, mais nous allons garder le béton à l'instant):

F_i = d + e X_i + f Y_i + g Z_i

Dans:

|F_0|   |1, X_0, Y_0, Z_0|   |d|
|F_1| = |1, X_1, Y_1, Z_1| * |e|
|F_2|   |1, X_2, Y_2, Z_2|   |f|
|F_3|   |1, X_3, Y_3, Z_3|   |g|

Ou: F = G * m (je suis un geophysist, nous utilisons donc des G pour les "Fonctions de Green" et m pour "les Paramètres du Modèle". Habituellement, nous utiliserions d pour les "données" au lieu de F, en tant que bien.)

En python, ce serait traduit:

def invert(f, x, y, z):
    G = np.vstack([np.ones_like(x), x, y, z]).T
    m, _, _, _ = np.linalg.lstsq(G, f)

    d, e, f, g = m
    a = np.sqrt(d)
    b = np.sqrt(e)
    c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
    return a, b, c

Non-linéaire de la Solution

Vous pouvez également résoudre ce à l'aide de scipy.optimize, comme @Joe suggéré. Le plus accessible de la fonction dans scipy.optimize est scipy.optimize.curve_fit qui utilise un Levenberg-Marquardt méthode par défaut.

De Levenberg-Marquardt est une "escalade" de l'algorithme de bien, il va de descente, dans ce cas, mais le terme est utilisé de toute façon). Dans un sens, vous faire une estimation initiale des paramètres du modèle (tous ceux que, par défaut, dans scipy.optimize) et de suivre la pente de observed - predicted dans votre espace de paramètre descente vers le bas.

Mise en garde: la Cueillette de la droite non-linéaire de la méthode d'inversion, estimation initiale, et le réglage des paramètres de la méthode est très bien un "dark art". Vous ne l'apprendre par le faire, et il y a beaucoup de situations où les choses ne fonctionnent pas correctement. De Levenberg-Marquardt est une bonne méthode si votre paramètre de l'espace est assez lisse (ce qui devrait être). Il y a beaucoup d'autres (y compris les algorithmes génétiques, réseaux de neurones, etc, en plus de méthodes plus courantes comme le recuit simulé) qui sont mieux dans d'autres situations. Je ne vais pas plonger dans la partie ici.

Il y est une commune de la chasse aux sorcières que certains d'optimisation des boîtes à outils d'essayer de corriger pour que scipy.optimize n'essayez pas de poignée. Si vos paramètres du modèle ont une amplitude différente (par exemple,a=1, b=1000, c=1e-8), vous aurez besoin de redimensionner les choses de sorte qu'ils sont de même importance. Sinon scipy.optimizes '"escalade" des algorithmes (comme LM) de ne pas calculer avec précision de l'estimation du gradient local, et donnera sauvagement des résultats inexacts. Pour l'instant, je suis en supposant que a, b, et c relativement similaire des grandeurs. Aussi, sachez que presque tous les non-linéaires méthodes, vous devrez faire une estimation initiale, et sont sensibles à cette supposition. Je pars it out ci-dessous (il suffit de passer dans le p0 kwarg à curve_fit) car la valeur par défaut a, b, c = 1, 1, 1 est une assez précise deviner pour a, b, c = 3, 2, 1.

Avec les mises en garde de la route, curve_fit s'attend à être passé d'une fonction, d'un ensemble de points où les observations ont été faites (comme un seul ndim x npoints tableau), et les valeurs observées.

Donc, si nous écrire la fonction comme ceci:

def func(x, y, z, a, b, c):
    f = (a**2
         + x * b**2
         + y * a * b * np.cos(c)
         + z * a * b * np.sin(c))
    return f

Nous aurons besoin pour l'envelopper d'accepter légèrement différents arguments avant de les passer à curve_fit.

En un mot:

def nonlinear_invert(f, x, y, z):
    def wrapped_func(observation_points, a, b, c):
        x, y, z = observation_points
        return func(x, y, z, a, b, c)

    xdata = np.vstack([x, y, z])
    model, cov = opt.curve_fit(wrapped_func, xdata, f)
    return model

Stand-alone Exemple des deux méthodes:

Pour vous donner la pleine mise en œuvre, voici un exemple qui

1. génère aléatoirement distribué des points pour évaluer la fonction,
2. évalue la fonction sur ces points (en utilisant l'ensemble des paramètres du modèle),
3. ajoute du bruit pour les résultats,
4. et puis inverse pour les paramètres du modèle en utilisant à la fois linéaire et non-linéaire des méthodes décrites ci-dessus.

---

import numpy as np
import scipy.optimize as opt
def main():
nobservations = 4
a, b, c = 3.0, 2.0, 1.0
f, x, y, z = generate_data(nobservations, a, b, c)
print 'Linear results (should be {}, {}, {}):'.format(a, b, c)
print linear_invert(f, x, y, z)
print 'Non-linear results (should be {}, {}, {}):'.format(a, b, c)
print nonlinear_invert(f, x, y, z)
def generate_data(nobservations, a, b, c, noise_level=0.01):
x, y, z = np.random.random((3, nobservations))
noise = noise_level * np.random.normal(0, noise_level, nobservations)
f = func(x, y, z, a, b, c) + noise
return f, x, y, z
def func(x, y, z, a, b, c):
f = (a**2
+ x * b**2
+ y * a * b * np.cos(c)
+ z * a * b * np.sin(c))
return f
def linear_invert(f, x, y, z):
G = np.vstack([np.ones_like(x), x, y, z]).T
m, _, _, _ = np.linalg.lstsq(G, f)
d, e, f, g = m
a = np.sqrt(d)
b = np.sqrt(e)
c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
return a, b, c
def nonlinear_invert(f, x, y, z):
# "curve_fit" expects the function to take a slightly different form...
def wrapped_func(observation_points, a, b, c):
x, y, z = observation_points
return func(x, y, z, a, b, c)
xdata = np.vstack([x, y, z])
model, cov = opt.curve_fit(wrapped_func, xdata, f)
return model
main()