I will be very grateful if anybody can tell me why I`m getting this
error:
raise TypeError("cannot create mpf from " + repr(x))
TypeError: cannot create mpf from (nan, nan)
Full code:
from mpmath import *
from scipy.integrate import quad
mp.dps = '10'
def calkalicznik(r):
wynikcalkilicznik=quad(lambda s:((mpf('-2')/power(s,mpf('3')))*(M1*
(mpf('-3')*(M2/power(s,mpf('3'))-M3/power(s,mpf('5'))+M4/power(s,mpf
('6'))+M5/power(s,mpf('8'))-(M6/power(s,mpf('5'))+M7/power(s,mpf
('8'))))-(s*(mpf('-3')*M2/power(s,mpf('4'))+mpf('5')*M3/power(s,mpf
('6'))-mpf('6')*M4/power(s,mpf('7'))-mpf('8')*M5/power(s,mpf('9')))))*
(M2/power(s,mpf('3'))-M3/power(s,mpf('5'))+M4/power(s,mpf('6'))+M5/
power(s,mpf('8'))-mpf('1'))+(D0/kB*T)*((b*(mpf('1')/sigma1a-M8/power
(s,mpf('4'))+M9/power(s,mpf('6')))+a*(mpf('1')/sigma1b-M10/power(s,mpf
('4'))+M11/power(s,mpf('6'))))/(a+b)+(mpf('4')*a*b*(M12/s-mpf('0.25')
*M13/(power(s,mpf('3')))+(mpf('75')/mpf('8'))*M14/(power(s,mpf
('7')))))/power((a+b),mpf('2')))*(mpf('-1')/mpf('6'))*AH*(mpf('-2')
*a*s*b/power(power(s,mpf('2'))-power((a+b),mpf('2')),mpf('2'))+mpf
('-2')*a*s*b/power(power(s,mpf('2'))-power((a-b),mpf('2')),mpf('2'))+
((mpf('2')*s/power(s,mpf('2'))-power((a-b),mpf('2'))-(mpf('2')*s*power
(s,mpf('2'))-power((a+b),mpf('2'))/power(power(s,mpf('2'))-power((a-
b),mpf('2')),mpf('2')))))*power(s,mpf('2'))-power((a-b),mpf('2'))/power
(s,mpf('2'))-power((a+b),mpf('2'))))/((mpf('2')/mpf('15'))*power(s,mpf
('2'))*power(S,mpf('2'))*tals*power((M2/power(s,mpf('3'))-M3/power
(s,mpf('5'))+M4/power(s,mpf('6'))+M5/power(s,mpf('8'))-mpf('1')),mpf
('2'))+D0*((b*(mpf('1')/sigma1a-M8/power(s,mpf('4'))+M9/power(s,mpf
('6')))+a*(mpf('1')/sigma1b-M10/power(s,mpf('4'))+M11/power(s,mpf
('6'))))/(a+b)+(mpf('4')*a*b*(M12/s-mpf('0.25')*M13/(power(s,mpf('3')))
+(mpf('75')/mpf('8'))*M14/(power(s,mpf('7')))))/power((a+b),mpf
('2'))))),inf, r)
return wynikcalkilicznik
def licznikimianownik():
#wynik=1
wynik = quad(lambda r: ((mpf('1')/(power(e,(calkalicznik(r))))) /
(power(r,mpf('2'))*(power(r,mpf('2'))*(M1*power((M2/power(r,mpf('3'))-
M3/power(r,mpf('5'))+M4/power(r,mpf('6'))+M5/power(r,mpf('8'))-mpf
('1')),mpf('2')))+D0*((b*(mpf('1')/sigma1a-M8/power(r,mpf('4'))+M9/
power(r,mpf('6')))+a*(mpf('1')/sigma1b-M10/power(r,mpf('4'))+M11/power
(r,mpf('6'))))/(a+b)+(mpf('4')*a*b*(M12/r-mpf('0.25')*M13/(power(r,mpf
('3')))+(mpf('75')/mpf('8'))*M14/(power(r,mpf('7')))))/power((a+b),mpf
('2')))))), (a+b), inf)
print 'LICZNIKMIANOWNIK:'
print wynik
return wynik
#dokladnosc
print e
#Stale
AH=power('10','-20')
u=power('10','-3')
T=mpf('293')
kB=mpf('1.38')*power('10','-23')
#Zmienne
a=mpf('1')*power('10','-6')
b=mpf('1')*power('10','-6')
D=mpf('1.7')
c=mpf('0.25')
daa=power('10','-4')
dpa=a
dab=daa
dpb=b
eps=mpf('0.1')
rop=mpf('1000')
lambd=b/a
v=u/rop
S=power((eps/v),mpf('0.5'))
tals=mpf('0.17')/S
#rownania
D0=(kB*T/mpf('6')*pi*u)*(mpf('1')/a+mpf('1')/b)
temp1=(daa/dpa)
temp2=D-mpf('3')
EPSpa=mpf('1')-c*power(temp1,temp2)
temp1=(dab/dpb)
EPSpb=mpf('1')-c*power(temp1,temp2)
temp3=mpf('8')/(mpf('1')-EPSpa)-mpf('3')
kappaba=power(dpa,'2')/mpf('72') * (mpf('3') + (mpf('3')/(mpf('1')-
EPSpa)) - (mpf('3') * power(temp3, (mpf('0.5')))))
#print kappaba
temp3=mpf('8')/(mpf('1')-EPSpb)-mpf('3')
kappabb=power(dpb,'2')/mpf('72') * (mpf('3') + (mpf('3')/(mpf('1')-
EPSpb)) -(mpf('3') * power(temp3, (mpf
('0.5')))))
#print kappabb
print 'zetaa'
zetaa=daa/(mpf('2')*power(kappaba,(mpf('1')/mpf('2'))))
print zetaa
zetab=dab/(mpf('2')*power(kappabb,mpf('0.5')))
print zetab
print 'Q1a'
Q1a=((mpf('30')+mpf('15')*power(zetaa,mpf('2'))+power(zetaa,mpf('4')))
*tanh(zetaa)-mpf('30')*zetaa-mpf('5')*power(zetaa,mpf('2')))/((mpf
('30')+mpf('10')*power(zetaa,mpf('2'))+power(zetaa,mpf('4')))*tanh
(zetaa)-mpf('30')*zetaa)
print Q1a
Q1b=((mpf('30')+mpf('15')*power(zetab,mpf('2'))+power(zetab,mpf('4')))
*tanh(zetab)-mpf('30')*zetab-mpf('5')*power(zetab,mpf('2')))/((mpf
('30')+mpf('10')*power(zetab,mpf('2'))+power(zetab,mpf('4')))*tanh
(zetab)-mpf('30')*zetab)
print 'Q1b'
print Q1b
Q2a=((mpf('3')*power(zetaa,mpf('2'))+power(zetaa,mpf('4')))*tanh
(zetaa)-mpf('3')*power(zetaa,mpf('2')))/((mpf('30')+mpf('10')*power
(zetaa,mpf('2'))+power(zetaa,mpf('4')))*tanh(zetaa)-mpf('30')*zetaa)
print 'Q2a'
print Q2a
Q2b=((mpf('3')*power(zetab,mpf('2'))+power(zetab,mpf('4')))*tanh
(zetab)-mpf('3')*power(zetab,mpf('2')))/((mpf('30')+mpf('10')*power
(zetab,mpf('2'))+power(zetab,mpf('4')))*tanh(zetab)-mpf('30')*zetab)
print 'Q2b'
print Q2b
sigma1a=(mpf('2')*power(zetaa,mpf('2'))*(zetaa-tanh(zetaa)))/(mpf('2')
*power(zetaa,mpf('3'))+mpf('3')*(zetaa-tanh(zetaa)))
print 'sigma1a'
print sigma1a
sigma1b=(mpf('2')*power(zetab,mpf('2'))*(zetab-tanh(zetab)))/(mpf('2')
*power(zetab,mpf('3'))+mpf('3')*(zetab-tanh(zetab)))
print 'sigma1b'
print sigma1b
sigma2a=(mpf('6')*(power(zetaa,mpf('2'))-mpf('2'))*(zetaa-tanh(zetaa))-
mpf('4')*power(zetaa,mpf('2')))/(mpf('2')*power(zetaa,mpf('3'))+mpf
('3')*(zetaa-tanh(zetaa)))
print 'sigma2a'
print sigma2a
sigma2b=(mpf('6')*(power(zetab,mpf('2'))-mpf('2'))*(zetab-tanh(zetab))-
mpf('4')*power(zetab,mpf('2')))/(mpf('2')*power(zetab,mpf('3'))+mpf
('3')*(zetab-tanh(zetab)))
print 'sigma2b'
print sigma2b
sigma3a=(power(zetaa,mpf('2'))*((mpf('3')+power(zetaa,mpf('2')))*
(zetaa-tanh(zetaa))-power(zetaa,mpf('3'))))/((mpf('30')+mpf('10')*power
(zetaa,mpf('2'))+power(zetaa,mpf('4')))*(zetaa-tanh(zetaa))-mpf('10')
*power(zetaa,mpf('3'))-power(zetaa,mpf('5')))
print 'sigma3a'
print sigma3a
sigma3b=(power(zetab,mpf('2'))*((mpf('3')+power(zetab,mpf('2')))*
(zetab-tanh(zetab))-power(zetab,mpf('3'))))/((mpf('30')+mpf('10')*power
(zetab,mpf('2'))+power(zetab,mpf('4')))*(zetab-tanh(zetab))-mpf('10')
*power(zetab,mpf('3'))-power(zetab,mpf('5')))
print 'sigma3b'
print sigma3b
sigma4a=((mpf('30')+mpf('15')*power(zetaa,mpf('2'))+power(zetaa,mpf
('4')))*(zetaa-tanh(zetaa))-mpf('10')*power(zetaa,mpf('3'))-power
(zetaa,mpf('5')))/((mpf('30')+mpf('10')*power(zetaa,mpf('2'))+power
(zetaa,mpf('4')))*(zetaa-tanh(zetaa))-mpf('10')*power(zetaa,mpf('3'))-
power(zetaa,mpf('5')))
print 'sigma4a'
print sigma4a
sigma4b=((mpf('30')+mpf('15')*power(zetab,mpf('2'))+power(zetab,mpf
('4')))*(zetab-tanh(zetab))-mpf('10')*power(zetab,mpf('3'))-power
(zetab,mpf('5')))/((mpf('30')+mpf('10')*power(zetab,mpf('2'))+power
(zetab,mpf('4')))*(zetab-tanh(zetab))-mpf('10')*power(zetab,mpf('3'))-
power(zetab,mpf('5')))
print 'sigma4b'
print sigma4b
sigma5a=(mpf('0.2')*power(zetaa,mpf('2'))*((mpf('15')+mpf('6')*power
(zetaa,mpf('2')))*(zetaa-tanh(zetaa))-mpf('5')*(power(zetaa,mpf
('3')))))/((mpf('315')/mpf('4')+(mpf('126')/mpf('4'))*power(zetaa,mpf
('2'))+(mpf('1')/mpf('5'))*power(zetaa,mpf('4')))*(zetaa-tanh(zetaa))-
(mpf('105')/mpf('4')*(power(zetaa,mpf('3')))))
print 'sigma5a'
print sigma5a
sigma5b=(mpf('0.2')*power(zetab,mpf('2'))*((mpf('15')+mpf('6')*power
(zetab,mpf('2')))*(zetab-tanh(zetab))-mpf('5')*(power(zetab,mpf
('3')))))/((mpf('315')/mpf('4')+(mpf('126')/mpf('4'))*power(zetab,mpf
('2'))+(mpf('1')/mpf('5'))*power(zetab,mpf('4')))*(zetab-tanh(zetab))-
(mpf('105')/mpf('4')*(power(zetab,mpf('3')))))
print 'sigma5b'
print sigma5b
sigma6a=(((mpf('15')+mpf('6')*power(zetaa,mpf('2')))*(zetaa-tanh
(zetaa))-mpf('5')*(power(zetaa,mpf('3')))))/(power(zetaa,mpf('2'))*
(zetaa-tanh(zetaa)))
print 'sigma6a'
print sigma6a
sigma6b=(((mpf('15')+mpf('6')*power(zetab,mpf('2')))*(zetab-tanh
(zetab))-mpf('5')*(power(zetab,mpf('3')))))/(power(zetab,mpf('2'))*
(zetab-tanh(zetab)))
print 'sigma6b'
print sigma6b
sigma7a=(mpf('20')*((mpf('3')+(mpf('12')/mpf('5'))*power(zetaa,mpf
('2'))+(mpf('1')/mpf('2'))*power(zetaa,mpf('4')))*(zetaa-tanh(zetaa))-
power(zetaa,mpf('3'))-(mpf('4')/mpf('10'))*power(zetaa,mpf('5'))))/
((power(zetaa,mpf('2')))*(mpf('2')*power(zetaa,mpf('3'))+mpf('3')*
(zetaa-tanh(zetaa))))
print 'sigma7a'
print sigma7a
sigma7b=(mpf('20')*((mpf('3')+(mpf('12')/mpf('5'))*power(zetab,mpf
('2'))+(mpf('1')/mpf('2'))*power(zetab,mpf('4')))*(zetab-tanh(zetab))-
power(zetab,mpf('3'))-(mpf('4')/mpf('10'))*power(zetab,mpf('5'))))/
((power(zetab,mpf('2')))*(mpf('2')*power(zetab,mpf('3'))+mpf('3')*
(zetab-tanh(zetab))))
print 'sigma7b'
print sigma7b
Sa=sigma2a/sigma1a
Sb=sigma2b/sigma1b
#stale dodatkowe
print 'start od M1'
M1=mpf('2')*power(S,mpf('2'))*tals/mpf('15')
print M1
M2=(mpf('5')/mpf('2'))*(Q2a*power(a,mpf('3'))+Q2b*power(b,mpf('3')))
print M2
M3=(mpf('3')/mpf('2'))*(Q1a*power(a,mpf('5'))+Q1b*power(b,mpf('5')))+
(mpf('5')/mpf('2'))*power(a,mpf('2'))*power(b,mpf('2'))*(a*Q2a*Sb
+b*Q2b*Sb)
print M3
M4=(mpf('25')*(Q2a*power(a,mpf('3'))*Q2b*power(b,mpf('3'))))
print M4
M5=(mpf('25')*(Q2a*power(a,mpf('3'))*Q2b*power(b,mpf('3')))*(Sa*power
(a,mpf('2'))+Sb*power(b,mpf('2'))))
print M5
M6=Q1a*power(a,mpf('5'))+Q1b*power(b,mpf('5'))+(mpf('5')/mpf('3'))
*power(a,mpf('2'))*power(b,mpf('2'))*(a*Q2a*Sb+b*Q2b*Sa)
print M6
M7=(mpf('25')/mpf('3'))*Q2a*power(a,mpf('3'))*Q2b*power(b,mpf('3'))*
(Sa*power(a,mpf('2'))+Sb*power(b,mpf('2')))
print M7
M8=(mpf('15')/mpf('4'))*a*sigma3b*power(b,mpf('3'))
print M8
M9=((mpf('1')/mpf('4'))*sigma2b*Sb+(mpf('9')/mpf('2'))*sigma4b-(mpf
('63')/mpf('10'))*sigma5b-(mpf('9')/mpf('20'))*sigma7b)*a*power(b,mpf
('5'))+(mpf('15')/mpf('2'))*sigma3b*Sa*power(a,mpf('3'))*power(b,mpf
('3'))
print M9
M10=(mpf('15')/mpf('4'))*sigma3a*(mpf('1')/lambd)*power(b,mpf('4'))
print M10
M11=((mpf('1')/mpf('4'))*sigma2a*Sa+(mpf('9')/mpf('2'))*sigma4a-(mpf
('63')/mpf('10'))*sigma5a-(mpf('9')/mpf('20'))*sigma7a)*(power
(lambd,mpf('-1')))*power(b,mpf('6'))+(mpf('15')/mpf('2'))
*sigma3a*Sb*power(lambd,mpf('-3'))*power(b,mpf('6'))
print M11
M12=(mpf('3')/mpf('4'))*(a+b)
print M12
M13=(Sa*power(a,mpf('2'))+Sb*power(b,mpf('2')))*(a+b)
print M13
M14=sigma3a*sigma3b*power(a,mpf('3'))*power(b,mpf('3'))*(a+b)
print M14
print 'poczatek obliczen'
#s=mpf('0.0000000009')
#CC1=M2/power(s,mpf('3'))-M3/power(s,mpf('5'))+M4/power(s,mpf('6'))+M5/
power(s,mpf('8'))-(M6/power(s,mpf('5'))+M2/power(s,mpf('8')))
#print CC1
licznikimianownik()
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to sympy@googlegroups.com
To unsubscribe from this group, send email to sympy+unsubscr...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/sympy?hl=en
-~----------~----~----~----~------~----~------~--~---
