First of all a suggestion: You don't need to use mpf('integer'), just
integer is fine, because they are exact. For floats it makes sense to
use mpf('float').
Your integrals seem to be very hard to calculate, scipy gives a
warning and seems to return (nan, nan), so it could not calculate it.
I can't try your script, because it's broken due to forced 80
characters per line, just attach it to prevent this.
There seems to be a typo:
On Jun 27, 7:00 am, qMic <michal.zubkow...@gmail.com> wrote:
> 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)
^
Maybe it lacks a '-' here? This would justify why it doesn't work.
>
> 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()
Vinzent
--~--~---------~--~----~------------~-------~--~----~
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
-~----------~----~----~----~------~----~------~--~---
