python - 如何用 Python 求解两个数值方程？

Question

也许这是向您展示代码并尝试稍后解释的最简单方法：

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：在给定的容差内找不到根。（0.0448851 > 2.1684e-19）尝试另一个起点或调整参数。”

你帮助我解决这个问题让我非常感激。这是我硕士论文的核心。

---

首先，感谢您的快速反应。

但我不认为这是正确的方法。你是对的，我忘记了n_brech=1.5。今天，我发现，gw-gw2可以解决解析。所以我只需要一个解决方案p-p2。ref和的解决方案sig可以在 PI 和sigma_aussen, roh_aussen, lambda_rek-PI 之间，通过 的变化，这就是我开始在 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)

我不明白为什么它不起作用...

Answer 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]从 at而不是 at开始优化[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()

...产生以下输出：

...从中可以看出，在感兴趣的范围内，沿曲线有许多局部最小值。因此，我对不同的方法（和不同的软件）找到不同的最小值并不感到惊讶。