问题描述
我是 stackoverflow 的新手，所以希望这篇文章遵循所有规则并且在正确的位置。如果有任何问题，请告诉我，这样我就不会再这样做了！此外，我希望一些 LaTeX 符号适用于希腊字母和数学运算符。
这里...
作为我的博士项目的一部分，我正在尝试从数学论文中复制图表。这涉及在 4 个未知数（a、w、α、β）中求解一组 4 个非线性方程，对应于代码中。这些方程具有第 5 个变量 \sigma'（代码中的 dsigma），但是我们在求解方程之前为此指定了一个值。我还需要定义功率
P = 2 * a^2 * w。
我想在 \exp(-4) 到 \exp(-0.5) 的范围内为 \sigma' 绘制 P 与 \log(\sigma') 的关系图。因此，这要求我遍历该范围内的 \sigma' 值并求解非线性方程，以最终确定 \sigma' 的每个值的 P。这个问题的难点在于方程包含一些积分项以及这些积分的导数。 The Equations、Integral Terms
现在我将解释我的实现。我将完整地按顺序包含代码，将其分解以解释一些内容。代码的第一位是序言，我在其中定义了一些将要使用的常量。还有一个评论是关于 x 的元素如何表示我提到的未知数。
import scipy.integrate as int
import scipy.optimize
import numpy as np
import matplotlib.pyplot as plt
eta_0 = pow(10,-3) # background packing fraction
g = (3 - eta_0) / pow(1 - eta_0,3) + np.log(eta_0) # constant appears in equations
d = 0.001 # finite difference
# x = [ a,w,alpha,beta ]
接下来的几行代码定义了函数，我将调用这些函数来计算方程中的积分项
# fns to define the integral terms
# uses 200 to approximate infinity in the integrals
def Omega(x):
return int.quad(lambda y: pow(1 / np.cosh(y / x[3]),2) * pow(1 / np.cosh(y / x[1]),2),200)[0]
def Xi(x):
return int.quad(lambda y: 2 * ((4 - 2 * eta_0 - 2 * x[2] * pow(1 / np.cosh(y),2))/pow(1 - eta_0 - (x[2] * pow(1 / np.cosh(y),2)),2) - (4 - 2 * eta_0) / pow(1 - eta_0,200)[0]
def Theta(x):
return int.quad(lambda y: (eta_0 * np.log(1 + x[2] / eta_0 * pow(1 / np.cosh(y),2)) + x[2] * pow(1 / np.cosh(y),2) * np.log(eta_0 + x[2] * pow(1 / np.cosh(y),2))),200)[0]
现在我定义函数来计算积分项的导数。我使用有限差分近似。
# Functions computing the finite difference derivatives of Omega,Xi and Theta
def dOmega_dw(x): # finite difference
x2 = np.zeros(4) # initialise x2
for j in range(len(x)):
x2[j] = x[j]
x2[1] = x[1] + d # x2 is x perturbed in the direction of derivative
return (Omega(x2) - Omega(x)) / d
def dOmega_dbeta(x):
x2 = np.zeros(4) # initialise x2
for j in range(len(x)):
x2[j] = x[j]
x2[3] = x[3] + d # x2 is x perturbed in the direction of derivative
return (Omega(x2) - Omega(x)) / d
def dXi_dalpha(x):
x2 = np.zeros(4) # initialise x2
for j in range(len(x)):
x2[j] = x[j]
x2[2] = x[2] + d # x2 is x perturbed in the direction of derivative
return (Xi(x2) - Xi(x)) / d
def dTheta_dalpha(x):
x2 = np.zeros(4) # initialise x2
for j in range(len(x)):
x2[j] = x[j]
x2[2] = x[2] + d # x2 is x perturbed in the direction of derivative
return (Theta(x2) - Theta(x)) / d
我想为每个 \sigma' 求解的方程组由下式给出：
# fn to compute system functions evaluated at x
def fn_vals(x,dsigma):
f = np.zeros(4)
f[0] = 1 - 3 * x[2] * x[1] * (Omega(x) - x[1] * dOmega_dw(x))
f[1] = dsigma + 1 / (2 * pow(x[1],2)) - (x[2] / x[1]) * (2 * Omega(x) - x[1] * dOmega_dw(x))
f[2] = 4 * pow(x[0],2) * x[2] * (Omega(x) - x[3] * dOmega_dbeta(x)) - x[3] * (x[2] * dXi_dalpha(x) - Xi(x)) - 4 * x[3] * (x[2] * dTheta_dalpha(x) - Theta(x))
f[3] = 4 * x[2] * pow(x[0],2) * dOmega_dbeta(x) - Xi(x) - 4 * Theta(x) + 4 * x[2] * (1 + g)
return f
最后，这是我在我之前指定的范围内求解 \sigma' 值的方程，然后确定 P 的方法。它还打印解 x，以及在方程组 fn_vals(x) 上评估的那些解.然后剧情就画好了。
# Main Method
plot_points = 50 # number of sigma' used from linspace
dsigma = np.linspace(np.exp(-4),np.exp(-0.5),plot_points) # linear space of sigma'
x0 = np.array([5,2,.5,2]) # Initial Guess
P = np.zeros(len(dsigma)) # Initialise Power array
for i in range(len(dsigma)):
x = scipy.optimize.fsolve(fn_vals,x0,args=(dsigma[i]))
#x = scipy.optimize.newton_krylov(lambda y: fn_vals(y,dsigma[i]),x0)
print(str(i) + ' x: ' + str(x))
print(str(i) + ' f(x): ' + str(fn_vals(x,dsigma[i])))
P[i] = 2 * x[0] ** 2 * x[1]
plt.plot(np.log(dsigma),P,color='red')
plt.xlabel('logσ\'')
plt.ylabel('P')
plt.ylim((0,100))
plt.show()
问题是代码产生警告信息，它产生的图形与我期望看到的不同，即论文中的图形。 Graph my code produces。 Graph I am aiming for (red curve)。我试图实现的曲线的 x 轴实际上被错误地标记为 \log(\sigma) 而不是 \log(\sigma')，但这无关紧要。
因此，我的问题是如何修复这些错误消息，更一般地说，我的整体实现是否存在问题？值得一提的是，随着 for 循环达到更大的 \sigma' 值，代码会显着减慢。我看不出有什么好的理由，因为改变 \sigma' 只是稍微改变了其中一个方程中的一个项，因此不会随着它的增加而引入额外的复杂性（据我所知，那就是！！）>
控制台输出：
C:\Users\jr00427\Anaconda3\python.exe C:/Users/jr00427/PycharmProjects/pythonProject/fig1.py
0 x: [2.17223218 3.99095975 0.06229124 2.20338412]
0 f(x): [ 8.58901839e-12 -1.41724826e-12 5.19292387e-11 7.05557834e-12]
1 x: [2.37561955 3.08515307 0.1046654 1.67791836]
1 f(x): [ 1.33013378e-10 7.77328202e-12 -5.00405273e-11 2.37623254e-10]
2 x: [2.53243759 2.62483192 0.14401699 1.44845902]
2 f(x): [ 6.08402217e-14 4.68194927e-13 -6.63180622e-12 3.13977733e-11]
3 x: [2.66874353 2.33907902 0.18012946 1.32517916]
3 f(x): [-5.18762810e-12 -3.58851837e-13 -2.00288675e-11 -8.10018719e-13]
C:/Users/jr00427/PycharmProjects/pythonProject/fig1.py:29: RuntimeWarning: invalid value encountered in log
return int.quad(lambda y: (eta_0 * np.log(1 + x[2] / eta_0 * pow(1 / np.cosh(y),2)))
C:/Users/jr00427/PycharmProjects/pythonProject/fig1.py:29: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze
the integrand in order to determine the difficulties. If the position of a
local difficulty can be determined (singularity,discontinuity) one will
probably gain from splitting up the interval and calling the integrator
on the subranges. Perhaps a special-purpose integrator should be used.
return int.quad(lambda y: (eta_0 * np.log(1 + x[2] / eta_0 * pow(1 / np.cosh(y),2)))
C:/Users/jr00427/PycharmProjects/pythonProject/fig1.py:19: IntegrationWarning: The occurrence of roundoff error is detected,which prevents
the requested tolerance from being achieved. The error may be underestimated.
return int.quad(lambda y: pow(1 / np.cosh(y / x[3]),2)
C:/Users/jr00427/PycharmProjects/pythonProject/fig1.py:24: IntegrationWarning: The occurrence of roundoff error is detected,which prevents
the requested tolerance from being achieved. The error may be
underestimated.
return int.quad(lambda y: 2 * ((4 - 2 * eta_0 - 2 * x[2] * pow(1 / np.cosh(y),2))
C:/Users/jr00427/PycharmProjects/pythonProject/fig1.py:29: IntegrationWarning: The occurrence of roundoff error is detected,which prevents
the requested tolerance from being achieved. The error may be underestimated.
return int.quad(lambda y: (eta_0 * np.log(1 + x[2] / eta_0 * pow(1 / np.cosh(y),2)))
4 x: [2.79424801 2.14206677 0.21321325 1.25261518]
4 f(x): [ 1.79067872e-12 5.53390667e-13 -1.31199496e-11 8.22670820e-11]
5 x: [2.91363137 1.99713794 0.24358556 1.20844003]
5 f(x): [-1.94449123e-10 -5.12074549e-11 -5.97554228e-10 -4.86269469e-09]
6 x: [3.02952519 1.88573693 0.27156854 1.18189116]
6 f(x): [ 2.65837463e-10 -9.57736668e-12 5.44213119e-10 -1.18931176e-09]
7 x: [3.14357663 1.79736488 0.29745628 1.16713633]
7 f(x): [ 3.66652264e-11 -8.64475158e-13 3.57718299e-10 -2.08828954e-10]
8 x: [3.25690492 1.72559242 0.32150642 1.16074372]
8 f(x): [ 2.51010324e-12 3.94073663e-13 -5.16942045e-12 1.63424829e-12]
9 x: [3.37032322 1.66624182 0.34394086 1.16056739]
9 f(x): [ 8.99047503e-12 1.85185201e-12 -1.26084920e-10 -1.56481939e-10]
10 x: [3.48445711 1.61647025 0.36494952 1.16520186]
10 f(x): [ 1.51456181e-10 1.84359750e-11 -1.86038251e-09 -3.35343220e-09]
11 x: [3.59981239 1.57427044 0.38469474 1.17369322]
11 f(x): [ 2.68745026e-11 5.78193049e-12 -8.87578899e-10 -8.24673663e-11]
12 x: [3.71681645 1.5381816 0.40331557 1.18537603]
12 f(x): [-2.52327048e-11 -7.29138971e-12 1.18753518e-09 -2.03906225e-10]
13 x: [3.83584463 1.50711355 0.42093148 1.19977648]
13 f(x): [-2.59998689e-11 -1.58172364e-11 2.84988366e-09 -1.40117162e-09]
14 x: [3.95723796 1.48023525 0.43764556 1.21655242]
14 f(x): [ 4.40729675e-11 -1.04632414e-11 8.94391228e-11 -1.26940236e-10]
15 x: [4.08131539 1.45690165 0.45354711 1.2354549 ]
15 f(x): [ 2.13939977e-12 -2.95874436e-13 4.08411083e-11 3.81827903e-12]
16 x: [4.20838276 1.43660421 0.4687139 1.25630271]
16 f(x): [-1.07911458e-11 1.32516220e-12 -3.45652973e-09 1.35364253e-09]
17 x: [4.33873954 1.41893662 0.48321392 1.27896513]
17 f(x): [-1.60516045e-12 2.15716334e-13 8.93907171e-12 -1.02069464e-11]
18 x: [4.47268409 1.4035705 0.49710688 1.30335004]
18 f(x): [ 4.24849045e-12 -1.23717703e-12 1.20773169e-10 -1.39364076e-11]
19 x: [4.61051803 1.39023772 0.51044549 1.32939549]
19 f(x): [ 5.78759263e-13 6.64357458e-13 -1.60298441e-11 3.89155375e-11]
20 x: [4.75254994 1.37871755 0.52327641 1.35706372]
20 f(x): [-9.43467526e-13 -2.15827356e-12 -7.86481991e-13 -2.03732142e-10]
21 x: [4.89909869 1.36882688 0.53564116 1.3863368 ]
21 f(x): [-2.78221890e-13 3.82582854e-13 3.28626015e-11 1.19682042e-10]
22 x: [5.05049648 1.36041294 0.54757682 1.41721355]
22 f(x): [ 9.33837452e-11 1.42422740e-11 -8.25012636e-09 -6.23565910e-09]
23 x: [5.20709183 1.35334762 0.55911661 1.44970722]
23 f(x): [ 7.51620988e-13 3.42614825e-13 2.43605136e-11 -6.86162238e-11]
24 x: [5.36925249 1.34752316 0.5702904 1.48384394]
24 f(x): [ 7.73503483e-12 3.11162207e-12 1.21335253e-09 -6.22270235e-10]
25 x: [5.53736851 1.34284865 0.58112515 1.51966152]
25 f(x): [-2.74983147e-10 4.06052969e-12 2.11696025e-08 1.47733950e-08]
26 x: [5.71185528 1.33924738 0.59164527 1.55720875]
26 f(x): [ 2.51110244e-12 1.37390099e-12 5.69415626e-10 -2.89000823e-10]
27 x: [5.89315698 1.33665465 0.60187293 1.59654494]
27 f(x): [ 3.32889272e-12 5.01754194e-12 9.32416633e-09 -9.03774833e-10]
28 x: [6.08175009 1.33501608 0.61182831 1.63773973]
28 f(x): [ 3.05033776e-12 3.52717855e-12 5.88860027e-09 -4.02114786e-10]
29 x: [6.27814738 1.33428619 0.62152985 1.68087312]
29 f(x): [ 6.07291994e-13 1.09035003e-12 2.63212208e-09 -7.91580135e-11]
30 x: [6.48290222 1.33442729 0.63099445 1.72603569]
30 f(x): [ 1.54798396e-12 3.75588449e-12 5.64094460e-09 -3.88435950e-10]
31 x: [6.6966134 1.33540861 0.64023763 1.77332899]
31 f(x): [ 1.82942550e-12 5.65669733e-12 9.27171229e-09 -2.74568812e-10]
32 x: [6.91993047 1.33720551 0.64927366 1.8228661 ]
32 f(x): [-4.96491737e-13 1.39122047e-12 2.48073828e-09 -2.24294361e-10]
33 x: [7.15355982 1.33979896 0.65811575 1.8747724 ]
33 f(x): [-9.73665593e-13 5.70654635e-14 1.58705671e-09 6.05054673e-10]
34 x: [7.39827148 1.34317504 0.6667761 1.92918642]
34 f(x): [-2.40674147e-12 -2.71782596e-13 3.28928618e-09 4.33761471e-10]
35 x: [7.65490678 1.34732459 0.67526603 1.98626102]
35 f(x): [ 3.39795969e-11 5.65969493e-12 -5.47825794e-08 -3.73836073e-09]
36 x: [7.92438721 1.35224294 0.68359609 2.04616461]
36 f(x): [-3.82072152e-12 -1.14885879e-12 3.55177843e-09 -8.76771544e-10]
37 x: [8.20772441 1.35792973 0.69177608 2.10908273]
37 f(x): [-3.35331762e-12 -4.57078819e-13 6.45956888e-09 -3.00442338e-10]
38 x: [8.50603162 1.36438876 0.69981519 2.17521979]
38 f(x): [ 4.41535697e-12 4.16777723e-13 -3.25887761e-09 8.72297790e-10]
39 x: [8.82053697 1.37162796 0.70772201 2.24480114]
39 f(x): [ 1.52505786e-11 -5.03452835e-12 -1.05404521e-08 1.63349156e-09]
40 x: [9.15259869 1.37965939 0.71550461 2.31807548]
40 f(x): [-1.54432023e-12 4.28101998e-12 -2.21988481e-08 -9.88716664e-10]
41 x: [9.50372283 1.38849928 0.72317057 2.39531758]
41 f(x): [-3.16764392e-11 8.52407034e-12 -8.59034888e-09 -3.73157683e-09]
42 x: [9.87558377 1.3981682 0.73072706 2.47683154]
42 f(x): [ 5.08562081e-11 -2.41295872e-12 -1.26152932e-07 6.87874646e-10]
43 x: [10.27004828 1.40869121 0.73818084 2.56295449]
43 f(x): [-8.52651283e-13 -3.46722651e-13 6.32552410e-09 3.05883319e-10]
44 x: [10.68920363 1.42009813 0.74553833 2.65406095]
44 f(x): [-5.94775340e-11 -3.21842553e-12 3.37865455e-07 5.49956347e-09]
45 x: [11.13539076 1.43242391 0.75280561 2.75056791]
45 f(x): [ 6.83897383e-14 6.65023592e-13 -1.48019073e-08 -3.61112029e-10]
46 x: [11.61124336 1.44570895 0.75998848 2.85294078]
46 f(x): [-4.19131396e-12 -6.77458090e-13 -8.53521698e-11 -3.44803297e-10]
47 x: [12.11973447 1.45999971 0.76709249 2.96170039]
47 f(x): [-9.09827769e-12 4.06008560e-13 3.41526398e-08 -1.99264605e-10]
48 x: [12.6642319 1.47534924 0.77412292 3.07743127]
48 f(x): [ 3.03934655e-12 -6.43929354e-14 -8.62290594e-08 -1.41965728e-09]
C:\Users\jr00427\Anaconda3\lib\site-packages\scipy\optimize\minpack.py:175: RuntimeWarning: The iteration is not making good progress,as measured by the
improvement from the last five Jacobian evaluations.
warnings.warn(msg,RuntimeWarning)
49 x: [8.96606493 1.47419045 0.71120114 2.44107678]
49 f(x): [-0.14923 0.08035658 -0.00437131 -0.00085175]
Process finished with exit code 0
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）