численно решить нелинейную систему уравнений на ограниченном интервале определения, используя python - PullRequest
Новая
       27

численно решить нелинейную систему уравнений на ограниченном интервале определения, используя python

0 голосов
/ 02 мая 2018

Допустим, следующая проблема: Мне нужно, чтобы функции, одна из которых параметризована t, и обе определены в независимой переменной z: f (t, z, * ags), g (z, * args).

Я хочу найти пару значений (t0, z0), чтобы

функции и их производные по z одинаковы f (t0, z0, * args) = g (z0, * args) а также df / dz (t0, z0, * args) = dg / dz (z0, * args).

Я знаю, что решение существует, и я получил разумную отправную точку (tS, zS). Однако, по крайней мере, одна из функций определена только в указанном интервале [zL .. zH] (что я знаю).

Мой вопрос сейчас таков: как лучше всего численно решить систему уравнений в python?

Я попробовал scipys fsolve, но, похоже, он не работает, я думаю, потому что он не может обработать ограниченный интервал определения. Я попробовал пакет diff_evolution, просто чтобы свести к минимуму составную функцию, но это выглядит полным перебором.

У меня есть выражения для всех функций и их производных (хотя они и сложные).

Конечно, должен существовать простой корень поиска Python, который способен решать систему из двух нелинейных уравнений, которые определены только на ограниченном интервале?!

По какой-то причине все, что я нахожу, может либо решить систему уравнений, не принимая во внимание пределы, либо может взять пределы, но решить только одно уравнение за один раз ...

Был бы очень признателен, если бы кто-то мог указать мне правильное направление, куда смотреть!

1 Ответ

0 голосов
/ 08 мая 2018

Этот вопрос требует допущений, что 1) существует несколько решений и 2) хотя бы одно решение находится в границах указанного интервала. Если эти предположения не выполняются, ваша проблема становится оптимизацией - минимизируя разницу между вашими функциями.

Я сгенерировал следующий пример, который соответствует предположениям: 

import numpy as np
from scipy.optimize import fsolve

def f(t,z):
    return t**2 + z + z**2 + np.sin(z)
def g(z):
    return z**2
def dfdz(t,z):
    return 1 + 2*z + np.cos(z)
def dgdz(z):
    return 2*z

def solve(x):
    t,z = x

    #residual array
    r = np.zeros(2)
    #match equations
    r[0] = (f(t,z) - g(z))**2
    #match derivatives
    r[1] = dfdz(t,z) - dgdz(z)

    return r

x0 = [10,-10]
x = fsolve(solve,x0)
t,z = x
print(t,z)

При этих начальных условиях ответом является [3.06998012761 -9.42477800292]. Требование ограниченных переменных возможно с помощью fsolve, добавляя остаток вне диапазона. Это довольно хакерское решение, которое не очень надежно. Мы можем связать z с диапазоном [-5,0] (где есть другой правильный ответ) с помощью этих модификаций: 

def solve(x,zL,zH):
    t,z = x

    #penalize deviation from range
    if z < zL:
        e = z - zL
    elif z > zH:
        e = z - zH
    else:
        e = 0

    #residual array
    r = np.zeros(2)
    #match equations
    r[0] = (f(t,z) - g(z))**2 + (e)**2
    #match derivatives
    r[1] = dfdz(t,z) - dgdz(z)

    return r


x0 = [10,-10]
x = fsolve(solve,x0,args=(-5,0))
t,z = x
print(t,z)

Новый ответ [1.77245385 -3.14159263].

Как упоминалось ранее, ограниченные переменные чаще рассматриваются в задачах оптимизации, поэтому пакет оптимизации будет обрабатывать их более интуитивно. В python есть стандартный scipy.optimize.minimize или гораздо более мощные пакеты, такие как Pyomo или GEKKO . Вот эквивалентная проблема в GEKKO: 

from gekko import GEKKO

#initialize a GEKKO model
m = GEKKO()

#add GEKKO variables
t = m.Var(value = -10)
z = m.Var(value = -10, lb=-5, ub=0)

#define the constraints
m.Equation(t**2 + z + z**2 + m.sin(z) == z**2)
m.Equation(1 + 2*z + m.cos(z) == 2*z)

#solve a system of non-dynamic equations
m.options.IMODE = 1 
m.solve(disp=False)

print(t.value,z.value)

[- 1.772454] [-3.142608]

...