1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
|
from sympy import *
from mpmath import glaisher, pi
β = Rational(1/8)
δ = 15
α = 0
σ = Rational(5, 6)
Δ = β * δ
s = 2**Rational(1, 12) * exp(-β) * glaisher**Rational(3, 2)
def f(B, θc, θ):
return (θc * Ei(-1 / (B * θc)) * exp(1 / (B * θc)) + (θ - θc) * Ei(1/(B * (θ - θc))) * exp(-1/(B * (θ - θc)))) / pi
def FH(θ0, θ):
return (θ**2 + θ0**2)**σ - (θ0**2)**σ
def FR(AL, AH, As, B, θc, θ0, θ):
return AL * f(B, θc, -θ) + AH * FH(θ0, θ) + sum([A * θ**(2 * i + 2) for i,A in enumerate(As)])
def F(AL, AH, As, B, θc, θ0, θ):
return FR(AL, AH, As, B, θc, θ0, θ) + AL * f(B, θc, θ)
def dfc(AL, B, m):
if m == 0:
return AL * θc * Ei(-1 / (B * θc)) * exp(1 / (B * θc)) / pi
elif m == 1:
return 0
else:
return AL * factorial(m) * factorial(m - 2) * B**(m-1) / pi
def dFc(AL, AH, As, B, θc, θ0, m):
θ = Symbol['θ']
return dfc(AL, B, m) + diff(FR(AL, AH, As, B, θc, θ0, θ), θ).subs(θ, θc)
|
