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

Question

Я никогда не использовал python, но Mathematica не может справиться с уравнением, которое я пытаюсь решить. Я пытаюсь решить для переменной «а» следующих уравнений, где s, c, mu и delta t являются известными параметрами.

Я пытался делать NSolve, решить, и т.д. в Mathematica, но он работает уже в течение часа, не повезло. Поскольку я не знаком с Python, есть ли способ использовать Python для решения этого уравнения для a?

Answer 1

Вы не найдете аналитического решения этих уравнений, потому что они трансцендентны, содержащие a как внутри, так и снаружи тригонометрической функции.

Я думаю, что проблема с численными решениями заключается в том, что диапазон допустимых значений для a ограничен arcsin. Поскольку arcsin определяется только для аргументов между -1 и 1 (если вы хотите, чтобы a был реальным), ваши формулы для alpha и beta требуют, чтобы a > s/2 и a > (s-c)/2.

В Python, вы можете найти нуль вашего третьего уравнения (переписать в виде f(a) = 0) с помощью brentq функции:

import numpy as np 
from scipy.optimize import brentq 

s = 10014.6 
c = 6339.06 
mu = 398600.0 
dt = 780.0 
def f(a): 
    alpha = 2*np.arcsin(np.sqrt(s/(2*a))) 
    beta = 2*np.arcsin(np.sqrt((s-c)/(2*a))) 
    return alpha - beta - (np.sin(alpha)-np.sin(beta)) - np.sqrt(mu/a**3)*dt 

a0 = max(s/2, (s-c)/2) 
a = brentq(f, a0, 10*a0) 

Edit:

Чтобы уточнить путь brentq(f,a,b) работы является что он ищет нуль f с интервалом [a,b]. Здесь мы знаем, что a составляет не менее max(s/2, (s-c)/2). Я просто догадался, что в 10 раз это была правдоподобная верхняя граница, и это работало для данных параметров. В более общем плане вам необходимо убедиться, что f меняет знак между a и b. Вы можете узнать больше о том, как работает функция в SciPy docs.

Answer 2

Я думаю, что его стоит изучить поведение функции, прежде чем пытаться ее решить. Не делая этого, вы не знаете, есть ли уникальное решение, множество решений или нет решения. (Самая большая проблема заключается в многих решениях, где числовые методы могут не дать вам решение, которое вам нужно/ожидайте, - и если вы вслепую его используете, могут возникнуть «плохие вещи»). Вы хорошо изучаете поведение, используя scipy и ipython. Это пример ноутбук, который делает это

# -*- coding: utf-8 -*- 
# <nbformat>3.0</nbformat> 

# <codecell> 

s = 10014.6 
c = 6339.06 
mu = 398600.0 
dt = 780.0 

# <codecell> 

def sin_alpha_2(x): 
    return numpy.sqrt(s/(2*x)) 
def sin_beta_2(x): 
    return numpy.sqrt((s-c)/(2*x)) 
def alpha(x): 
    return 2*numpy.arcsin(numpy.clip(sin_alpha_2(x),-0.99,0.99)) 
def beta(x): 
    return 2*numpy.arcsin(numpy.clip(sin_beta_2(x),-0.99,0.99)) 

# <codecell> 

def fn(x): 
    return alpha(x)-beta(x)-numpy.sin(alpha(x))+numpy.sin(beta(x)) - dt * numpy.sqrt(mu/numpy.power(x,3)) 

# <codecell> 

xx = numpy.arange(1,20000) 
pylab.plot(xx, numpy.clip(fn(xx),-2,2)) 

# <codecell> 

xx=numpy.arange(4000,10000) 
pylab.plot(xx,fn(xx)) 

# <codecell> 

xx=numpy.arange(8000,9000) 
pylab.plot(xx,fn(xx)) 

Это показывает, что мы ожидаем, чтобы найти решение с между 8000 и 9000. d изгиб на кривой около 5000 и более ранний раствор при температуре около 4000 обусловлен отсечением, необходимым для того, чтобы заставить арксин вести себя. Действительно, уравнение не имеет смысла ниже около а = 5000. (точное значение - a0, заданное в решении Rays). Это дает хороший диапазон, который можно использовать с методами в решении Rays.