Как решить два уравнения числовые по Python?
Возможно, это самый простой способ показать вам код и попытаться объяснить его позже:
sigma_aussen = 7.1
roh_aussen = 38
lambda_schreib = 532*10**-9
lambda_rek = 432*10**-9
sigma_aussen=radians(sigma_aussen)
roh_aussen=radians(roh_aussen)
def BraggMatch(sigma_aussen, roh_aussen, n_brech, lambda_schreib, lambda_rek):
gw2, sig_innen, ref_innen, sig_din, ref_din, p2, sig, ref =var('gw2 sig_innen ref_innen sig_din ref_din p2 sig ref')
sigma = asin(1/n_brech*sin(sigma_aussen))
roh = asin(1/n_brech*sin(roh_aussen))
gw = (sigma+roh)/2
sigma_din = (roh-sigma)/2
roh_din = (sigma-roh)/2
gw2 = (asin(1/n_brech*sin(sig))+asin(1/n_brech*sin(ref))/2)
print("gw2: ", gw2)
sig_innen = asin(1/n_brech*sin(sig))
ref_innen = asin(1/n_brech*sin(ref))
print ("sig_innen: ", sig_innen)
print ("ref_innen: ", ref_innen)
sig_din = gw2-sig_innen
ref_din = gw2-ref_innen
print ("ref_din: ", ref_din)
p = lambda_schreib/(n_brech*(sin(sigma_din)-sin(roh_din)))
print ("p: ",p)
p2 = lambda_rek/n_brech*(sin(sig_din)-sin(ref_din))
print ("p2: ", p2)
Winkel=nsolve([gw-gw2, p-p2],[sig,ref], [0,0])
return Winkel
Winkel = BraggMatch(sigma_aussen, roh_aussen, n_brech, lambda_schreib, lambda_rek)
Вот мое объяснение: BraggMatch - это метод, который возвращает два угла в радианах.
два уравнения: gw-gw2=0 и p-p2=0. gw и p - две известные переменные. неизвестны только две переменные, sig и ref. что две неизвестные переменные должны быть решены численно и возвращены в метод BraggMatch.
С программным обеспечением Maple, решение двух уравнений не проблема. Может быть, это может помочь, показывая вам решение: sig=0.064 и ref=0.734
Это ошибка: "ValueError: Не удалось найти root в пределах заданного допуска. (0.0448851 > 2.1684e-19) Попробуйте другую начальную точку или настройте аргументы".
Вы делаете меня очень благодарным, помогая мне решить эту проблему. Это ядро моей магистерской диссертации.
Во-первых, спасибо за вашу быструю реакцию.
Но я не думаю, что это правильный путь. Вы правы, я забыл n_brech=1.5
, Сегодня я узнаю, что gw-gw2
может решить аналитическую. Так что мне нужно только решение для p-p2
, Решения ref
а также sig
могут быть между PI и -PI, путем изменения sigma_aussen, roh_aussen, lambda_rek
Вот почему я начинаю искать решение в 0.
Здесь вы видите мой новый код:
n_brech = 1.5
sigma_aussen = 7.1
roh_aussnen = 38
lambda_schreib = 532*10**-9
lambda_rek = 432*10**-9
sigma_aussen=radians(sigma_aussen)
roh_aussen=radians(roh_aussen)
def BraggMatch(sigma_aussen, roh_aussen, n_brech, lambda_schreib, lambda_rek):
gw2, sig_innen, ref_innen, sig_din, ref_din, p2, sig, ref =var('gw2 sig_innen ref_innen sig_din ref_din p2 sig ref')
sigma = asin(1/n_brech*sin(sigma_aussen))
roh = asin(1/n_brech*sin(roh_aussen))
gw = (sigma+roh)/2
sigma_din = (roh-sigma)/2
roh_din = (sigma-roh)/2
gw2 = (asin(1/n_brech*sin(sig))+asin(1/n_brech*sin(ref))/2)
ref_list=solve(gw-gw2,ref)
ref=ref_list[0]
print("ref: ", ref)
gw2 = (asin(1/n_brech*sin(sig))+asin(1/n_brech*sin(ref))/2)
sig_innen = asin(1/n_brech*sin(sig))
ref_innen = asin(1/n_brech*sin(ref))
sig_din = gw2-sig_innen
ref_din = gw2-ref_innen
p = lambda_schreib/(n_brech*(sin(sigma_din)-sin(roh_din)))
p2 = lambda_rek/(n_brech*(sin(sig_din)-sin(ref_din)))
print ("p2: ", p2)
print ("p-p2: ", p-p2)
sig=fsolve(p-p2,0)
Winkel=[sig, ref]
return Winkel
Winkel=BraggMatch(sigma_aussen, roh_aussen, n_brech, lambda_schreib, lambda_rek)
Если вы хорошо разбираетесь в Maple, на изображении, которое я загружаю, вы видите проблему, которую мне нужно решить в Python, потому что здесь я могу визуализировать ее с помощью VPython. Код Maple возвращает углы в градусах, что не проблема.
Теперь я попытался написать код, как в Maple. Я получаю некоторые ошибки сейчас:
ZeroDivisionError
Вот код:
n_brech = 1.5
sigma_aussen = 7.1
roh_aussen = 38
lambda_schreib = 532*10**-9
lambda_rek = 432*10**-9
sigma_aussen=radians(sigma_aussen)
roh_aussen=radians(roh_aussen)
def BraggMatch(sigma_aussen, roh_aussen, n_brech, lambda_schreib, lambda_rek):
sig, ref =var('sig ref')
def sigma(n_brech, sigma_aussen):
return (asin(1/n_brech*sin(sigma_aussen)))
def roh(n_brech, roh_aussen):
return (asin(1/n_brech*sin(roh_aussen)))
def gw(n_brech, sigma_aussen, roh_aussen):
return ((sigma(n_brech, sigma_aussen)+roh(n_brech, roh_aussen))/2)
def sigma_din(n_brech, sigma_aussen, roh_aussen):
return (gw(n_brech, sigma_aussen, roh_aussen)-sigma(n_brech, sigma_aussen))
def roh_din(n_brech, sigma_aussen, roh_aussen):
return (gw(n_brech, sigma_aussen, roh_aussen)-roh(n_brech, roh_aussen))
def p(n_brech, sigma_aussen, roh_aussen, lambda_schreib):
return (lambda_schreib/(n_brech*(sin(sigma_din(n_brech, sigma_aussen, roh_aussen))-sin(roh_din(n_brech, sigma_aussen, roh_aussen)))))
return (nsolve([gw(n_brech, sigma_aussen, roh_aussen)-gw(n_brech, sig, ref), p(n_brech, sigma_aussen, roh_aussen, lambda_rek)-p(n_brech, sig, ref, lambda_schreib)], [sig, ref], [0, 0]))
Winkel=BraggMatch(sigma_aussen, roh_aussen, n_brech, lambda_schreib, lambda_rek)
Я не понимаю, почему это не работает...
1 ответ
Эта проблема также может быть решена численно, например, с scipy.optimize.minimize()
:
from scipy.optimize import minimize;
from math import radians, asin, sin;
sigma_aussen = 7.1
roh_aussen = 38
lambda_schreib = 532*10**-9
lambda_rek = 432*10**-9
sigma_aussen=radians(sigma_aussen)
roh_aussen=radians(roh_aussen)
n_brech = 1
def BraggMatch(sig,ref,sigma_aussen, roh_aussen, n_brech, lambda_schreib, lambda_rek):
sigma = asin(1/n_brech*sin(sigma_aussen))
roh = asin(1/n_brech*sin(roh_aussen))
gw = (sigma+roh)/2
sigma_din = (roh-sigma)/2
roh_din = (sigma-roh)/2
gw2 = (asin(1/n_brech*sin(sig))+asin(1/n_brech*sin(ref))/2)
sig_innen = asin(1/n_brech*sin(sig))
ref_innen = asin(1/n_brech*sin(ref))
sig_din = gw2-sig_innen
ref_din = gw2-ref_innen
p = lambda_schreib/(n_brech*(sin(sigma_din)-sin(roh_din)))
p2 = lambda_rek/n_brech*(sin(sig_din)-sin(ref_din))
return (gw-gw2)**2 + (p-p2)**2
def f(x,sigma_aussen, roh_aussen, n_brech, lambda_schreib, lambda_re):
return BraggMatch(x[0],x[1],sigma_aussen, roh_aussen, n_brech, lambda_schreib, lambda_re)
print(minimize(f,[0,0],args=(sigma_aussen, roh_aussen, n_brech, lambda_schreib, lambda_rek)).x)
Однако я должен был догадаться, что n_brech = 1
потому что его значение не было показано в вопросе. scipy.optimize.minimize()
действительно, находит решение ([ 0.3148574 0.1574287]
) но это не то же самое решение, что и Maple.
Не зная источника уравнений и того, что они означают, я не могу быть уверен, что они сформулированы правильно. Я мог бы, однако, проверить несколько решений. Начиная оптимизацию с [0.064, 0.734]
вместо того, чтобы в [0, 0]
дал результат [ 0.0340574 0.7190287]
,
Таким образом, кажется, что есть несколько решений для этих уравнений и (по крайней мере, для этого значения n_brech
), оптимизация найдет разные решения, если она будет запущена в разных точках.
(Далее РЕДАКТИРОВАТЬ)
Чтобы исследовать несколько решений, я использовал следующий код:
from scipy.optimize import minimize;
from math import radians, asin, sin, log;
sigma_aussen = 7.1
roh_aussen = 38
lambda_schreib = 532*10**-9
lambda_rek = 432*10**-9
sigma_aussen=radians(sigma_aussen)
roh_aussen=radians(roh_aussen)
n_brech = 1.5
def BraggMatch(sig,ref,sigma_aussen, roh_aussen, n_brech, lambda_schreib, lambda_rek):
sigma = asin(1/n_brech*sin(sigma_aussen))
roh = asin(1/n_brech*sin(roh_aussen))
gw = (sigma+roh)/2
sigma_din = (roh-sigma)/2
roh_din = (sigma-roh)/2
gw2 = (asin(1/n_brech*sin(sig))+asin(1/n_brech*sin(ref))/2)
sig_innen = asin(1/n_brech*sin(sig))
ref_innen = asin(1/n_brech*sin(ref))
sig_din = gw2-sig_innen
ref_din = gw2-ref_innen
p = lambda_schreib/(n_brech*(sin(sigma_din)-sin(roh_din)))
p2 = lambda_rek/n_brech*(sin(sig_din)-sin(ref_din))
return log((gw-gw2)**2 + (p-p2)**2)
import matplotlib
import itertools
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
matplotlib.rcParams['xtick.direction'] = 'out'
matplotlib.rcParams['ytick.direction'] = 'out'
delta = 0.001
sig = np.arange(0, 1, delta)
ref = np.arange(0, 1, delta)
XY = [(x,y) for x in sig for y in ref]
Z = [ BraggMatch(x,y,sigma_aussen, roh_aussen, n_brech, lambda_schreib, lambda_rek) for x,y in XY]
X = [x for x,y in XY]
Y = [y for x,y in XY]
Z = np.reshape(Z,(len(sig),-1))
plt.figure()
CS = plt.contour(sig, ref, Z)
plt.clabel(CS, inline=1, fontsize=10)
plt.title('Contour plot')
plt.xlabel('sig')
plt.ylabel('ref')
plt.show()
... который производит следующий вывод:
... из которого видно, что вдоль кривой в интересующем диапазоне много локальных минимумов. Таким образом, я не удивлен, что разные методы (и разное программное обеспечение) находят разные минимумы.