polinomios - Resolviendo ecuaciones polinómicas multivariadas simultáneas con python

python resolver ecuaciones diferenciales

Question

editar: la referencia de la que obtuve mis ecuaciones contenía un par de errores. Lo he arreglado aquí. ¡Las soluciones podrían tener sentido ahora!

Cuando un fluido de dos capas fluye sobre la topografía, existen varias soluciones diferentes dependiendo del tamaño relativo de la velocidad del flujo y la velocidad de la onda en el fluido.

Estos se denominan "supercríticos", "subcríticos" y "críticos" (los dos primeros me refiero aquí como "extracríticos").

Las siguientes ecuaciones definen las líneas límite entre el comportamiento crítico y extra-crítico en el espacio de parámetros (h, U0):

Quiero eliminar d_1c (es decir, no me importa de qué se trata) y encontrar soluciones a estas ecuaciones en (h, U_0) .

Factores simplificadores:

• Solo necesito respuestas para d_0 dado
• No necesito soluciones exactas , solo un esquema de las curvas de solución, por lo que esto se puede resolver de forma analítica o numérica.
• Solo quiero trazar la región (h, U0) = (0,0) a (0.5, 1).

Me gustaría resolver esto usando módulos disponibles en la distribución Enthought (numpy, scipy, sympy), pero realmente no sé por dónde empezar. Es la eliminación de la variable d1c lo que realmente me confunde.

Aquí están las ecuaciones en python:

def eq1(h, U0, d1c, d0=0.1):
f = (U0) ** 2 * ((d0 ** 2 / d1c ** 3) + (1 - d0) ** 2 / (1 - d1c - d0) ** 3) - 1
return f
def eq2(h, U0, d1c, d0=0.1):
f = 0.5 * (U0) ** 2 * ((d0 ** 2 / d1c ** 2) - (1 - d0) ** 2 / (1 - d1c - d0) ** 2) + d1c + (h - d_0)
return f

Estoy esperando una solución que tenga varias ramas de solución (no siempre físicas, pero no te preocupes) y se ve más o menos así:

¿Cómo hago para implementar esto?

Answer 1

Déjame lidiar con la eliminación de d1c primero. Imagine que logró masajear la primera ecuación para tener la forma d1c = f(U, h, d0) . Luego, sustituirías esto en la segunda ecuación y tendrías cierta relación entre U , h y d0 . Con d0 fijo, esto define una sola ecuación para dos variables, U y h , desde la cual puede, en principio, encontrar U para cualquier h dada. Lo cual parece ser lo que estás llamando soluciones, basado en tu último boceto. Las malas noticias son que no es fácil obtener d1c de cualquiera de sus ecuaciones. La buena noticia es que no es necesario.

fsolve puede tomar un sistema de ecuaciones, digamos dos ecuaciones que dependen de dos variables y le dan la solución. En este caso: arregle h ( d0 ya está fijo), y aliméntelo para fsolve el sistema que tiene, tratándolo como variables U0 y d1c . Registre el valor de U0 , repita para el siguiente valor de h , y así sucesivamente.

Tenga en cuenta que, contrariamente a la sugerencia de @duffymo, recomiendo usar fsolve , o al menos comenzar con él, y buscar otros solucionadores solo si se queda sin energía.

Una posible advertencia es que está esperando tener más de una solución para U0 dado que h : fsolve necesita una suposición inicial, y no hay una forma fácil de decirle que converja a una de las ramas de la solución. Si eso resulta ser un problema, mira el solucionador de brentq .

Una forma diferente sería observar que puede eliminar fácilmente U0 de su sistema. De esta forma obtendrá una ecuación única para h y d1c , la resolverá para d1c para cada valor de h , y luego usará cualquiera de sus ecuaciones originales para calcular U0 dado d1c y h .

Un ejemplo de uso de fsolve :

>>> from scipy.optimize import fsolve
>>> def f(x, p):
... return x**2 -p
...
>>> fsolve(f, 0.5, args=(2,))
array([ 1.41421356])
>>> 

Aquí args=(2,) es la sintaxis para decir fsolve que lo que realmente queremos resolver si f(x,2)=0 , y 0.5 es la conjetura inicial para el valor de x .

Answer 2

Semi-formalmente, el problema que está tratando de resolver es el siguiente: dado d0, resuelva la fórmula lógica "existe d1c tal que eq1 (h, U0, d1c, d0) = eq2 (h, U0, d1c, d0) = 0 "para h y U0.

Existe un algoritmo para reducir la fórmula a una ecuación polinómica "P (h, U0) = 0", se llama "eliminación del cuantificador" y generalmente se basa en otro algoritmo, "descomposición algebraica cilíndrica". Lamentablemente, esto no se implementa en sympy (aún).

Sin embargo, como U0 puede eliminarse fácilmente, hay cosas que puedes hacer con Sympy para encontrar tu respuesta. Empezar con

h, U0, d1c, d0 = symbols(''h, U0, d1c, d0'')
f1 = (U0) ** 2 * ((d0 ** 2 / d1c ** 3) + (1 - d0) ** 2 / (1 - d1c - d0 * h) ** 3) - 1
f2 = U0**2 / 2 * ((d0 ** 2 / d1c ** 2) + (1 - d0) ** 2 / (1 - d1c - d0 * h)) + d1c + d0 * (h - 1)

Ahora, elimine U0 de f1 e inserte el valor en f2 (lo hago "a mano" en lugar de con solve () para obtener una expresión más bonita):

U2_val = ((f1 + 1)/U0**2)**-1
f3 = f2.subs(U0**2, U2_val)

f3 solo depende de hyd1c. Además, dado que es una fracción racional, solo nos importa cuando su numerador va a 0, por lo que obtenemos una ecuación polinómica única en 2 variables:

p3 = fraction(cancel(f3))

Ahora, para un d0 dado, debería poder invertir p3.subs (d0, .1) numéricamente para obtener h (d1c), conectarlo de nuevo a U0 y hacer un diagrama paramétrico de (h, U0) como una función de d1c.

Answer 3

Utiliza solucionadores no lineales como Newton Raphson o BFGS para resolver ecuaciones simultáneas no lineales. Son sensibles a las condiciones iniciales y al acondicionamiento de las matrículas, por lo que se necesitará cierto cuidado.