Решение трансцендентного уравнения с помощью Mathematica и Python - PullRequest
Новая
       14

Решение трансцендентного уравнения с помощью Mathematica и Python

0 голосов
/ 28 апреля 2018

У меня есть следующая проблема с поиском корней нелинейного уравнения. Уравнение следующее:

tanh [5 * log [(2 / t) ^ (0.00990099) (1 + x) ^ (0.990099) (1-x) ^ (- 1)]] -x = 0

Решение с помощью NSolve, для {t, 0, 100} возвращает следующее с Mathematica :

plot

Это то, что я ожидал, построив полученные корни в зависимости от параметра времени в этом диапазоне. Теперь я попытался воспроизвести этот результат с помощью Python, используя scipy.optimize.root, но кажется, что мой код возвращает в качестве решения любое значение, которое я использую в качестве начального условия, следовательно, это ничто иное, как карта идентичности. Это также можно увидеть на картинке ниже, где я использовал начальное условие 0.7:

plot1

Я предоставил код ниже: 

import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import root

#Setting up the function
def delta(v,t):
    epsilon = 10**(-20)
    return np.tanh( 5*np.log( (2/(1.0*t+epsilon))**(0.00990099)*(1+v+epsilon)**(0.990099)*(1-v+epsilon)**(-1)))-v
#Setting up time paramerer
time = np.linspace(0, 101)
res = [root(delta, 0.7, args=(t, )).x[0] for t in time]
print res
plt.plot(time, res)
plt.savefig("plot.png")

Я не совсем уверен, правильно ли я использую scipy.optimize.root, поскольку функция выглядит нормально, насколько я ожидаю от ее поведения. Возможно, ошибка в том, как я передаю args?

1 Ответ

0 голосов
/ 29 апреля 2018

Методы поиска корней, которые начинаются с интервала скобок [a, b] (тот, где f (a) и f (b) имеют противоположные знаки), как правило, более устойчивы, чем методы, начинающиеся с одной точки x0 Вылет из. Причина в том, что у первых есть определенное поле для работы, и оно может уточнять его итеративно. Классический метод является классическим примером этого, но он медленный. SciPy реализует более сложные методы, такие как brentq . Здесь он работает нормально, с скобкой [-0,1, 0,1] (этого должно быть достаточно, если смотреть на график Mathematica).

Кроме того, t = 0 проблематично в уравнении, так как тогда оно даже не определено. Вместо этого положите небольшое положительное число, например 0,01. 

time = np.linspace(0.01, 101, 500)
res = [brentq(delta, -0.1, 0.1, args=(t, )) for t in time]

plot

...