I'm new to stackoverflow so hopefully this post follows all the rules and is in the correct place. If there's any problems let me know so I don't do it again! Also I hope a little LaTeX notation is ok for the greek letters and maths operators.
Here goes...
I'm attempting to reproduce a graph from a maths paper as part of my PhD project. This involves solving a set of 4 nonlinear equations in 4 unknowns (a, w, \alpha, \beta), corresponding to in the code. These equations feature a 5th variable, \sigma' (dsigma in the code), however we prescribe a value to this before solving the equations. I also need to define the power
P = 2 * a^2 * w.
I would like to plot a graph of P versus \log(\sigma') for \sigma' in the range \exp(-4) to \exp(-0.5). This therefore requires that I loop through values of \sigma' in that range and solve the nonlinear equations to ultimately determine P for each value of \sigma'. The difficulty of this problem is that the equations contain some integral terms as well as derivatives of these integrals. The Equations, Integral Terms
Now I'll explain my implementation. I'll include the code in it's entirety and in order, breaking it up to explain a few bits. The first bit of code is the preamble where I define a few constants that will be used. There is also a comment on how the elements of x represent the unknowns that I mentioned.
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 ]
The next few lines of code define functions that I'll call to compute the integral terms in the equations
# 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), 0, 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, 2)), 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))), 0, 200)[0]
Now I define functions to compute the derivatives of the integral terms. I use a finite difference approximation.
# 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
The system of equations I would like to solve for each \sigma' is then given by:
# 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
Finally, this is my method to solve the equations for values of \sigma' in the range I previously specified and then determine P. It also prints solutions, x, along with those solutions evaluated on system of equations fn_vals(x). Then the plot is drawn up.
# 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()
The problem is that the code produces warning messages and the graph it produces is nothing like the one that I expect to see i.e. the one from the paper. Graph my code produces. Graph I am aiming for (red curve). The x-axis of the curve i'm trying to achieve is actually mislabelled with \log(\sigma) rather than \log(\sigma'), but that's of no consequence.
My questions is therefore how do I fix these error messages and more generally, is there an issue with my overall implementation? Something worth mentioning is that the code slows down considerably as the for loop gets to larger \sigma' values. I see no good reason for this since changing the \sigma' simply changes a single term in one of the equations slightly, and so does not introduce additional complexity as it increases (as far as I can tell, that is!!)
Console output:
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)) + x[2] * pow(1 / np.cosh(y), 2) * np.log(eta_0 + x[2] * 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)) + x[2] * pow(1 / np.cosh(y), 2) * np.log(eta_0 + x[2] * 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) * pow(1 / np.cosh(y / x[1]), 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))/pow(1 - eta_0 - (x[2] * pow(1 / np.cosh(y), 2)), 2) - (4 - 2 * eta_0) / pow(1 - eta_0, 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)) + x[2] * pow(1 / np.cosh(y), 2) * np.log(eta_0 + x[2] * 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