読者です 読者をやめる 読者になる 読者になる

北海道苫小牧市出身の初老PGが書くブログ

永遠のプログラマを夢見る、苫小牧市出身のおじさんのちらしの裏

問題4.7と4.8の解答(不定積分と定積分)

math python scipy sympy

概説微分積分の問題4.7と4.8の解答。手計算に飽きたので現実逃避。

4.7 は sympy を使う。理論的に解けるのはわかってるけども、手計算じゃもはや無理。

7番目の問題にやたら時間がかかったんだけど、辛かったのはどの辺だろ。セオリー的にはt = tan(x/2)と置いて置換積分して有理式に持ち込むんだろうけど、sympyの実装がどうなってるかは見てない。

import sympy as sym
from sympy.abc import x

print sym.integrate(x**2 / sym.sqrt(4 - x**2), x)
print sym.integrate(x**2 / sym.sqrt(x**2 + 4), x)
print sym.integrate(1 / (x**2 + 4)**5, x)
print sym.integrate(1 / (x**2 + 2 * x * 10)**4, x)
print sym.integrate(1 / sym.sin(x)**5, x)
print sym.integrate(sym.log(x)**4, x)
print sym.integrate(sym.sin(x) / (sym.cos(x)**2 + sym.sin(x) + 1), x)
print sym.integrate(x * sym.sqrt((x - 1) * (x + 2)), x)
print sym.integrate(1 / (sym.exp(x) + sym.exp(-x) + 1), x)
print sym.integrate(sym.exp(x) * sym.sin(2 * x), x)
print sym.integrate(x**2 * sym.exp(x) * sym.cos(x), x)
print sym.integrate(x**2 * sym.sin(sym.sqrt(x)), x)
-x*(-x**2 + 4)**(1/2)/2 + 2*asin(x/2)
x*(x**2 + 4)**(1/2)/2 - 2*asinh(x/2)
(105*x**7 + 1540*x**5 + 8176*x**3 + 17856*x)/(98304*x**8 + 1572864*x**6 + 9437184*x**4 + 25165824*x**2 + 25165824) + 35*atan(x/2)/65536
-log(x)/64000000 + log(x + 20)/64000000 - (3*x**5 + 150*x**4 + 2200*x**3 + 6000*x**2 - 24000*x + 160000)/(9600000*x**6 + 576000000*x**5 + 11520000000*x**4 + 76800000000*x**3)
(3*cos(x)**3 - 5*cos(x))/(8*cos(x)**4 - 16*cos(x)**2 + 8) + 3*log(cos(x) - 1)/16 - 3*log(cos(x) + 1)/16
x*log(x)**4 - 4*x*log(x)**3 + 12*x*log(x)**2 - 24*x*log(x) + 24*x
-2*3**(1/2)*I*log(tan(x/2) - 1/2 - 3**(1/2)*I/2)*tan(x/2)/(9*tan(x/2) + 9) - 2*3**(1/2)*I*log(tan(x/2) - 1/2 - 3**(1/2)*I/2)/(9*tan(x/2) + 9) + 2*3**(1/2)*I*log(tan(x/2) - 1/2 + 3**(1/2)*I/2)*tan(x/2)/(9*tan(x/2) + 9) + 2*3**(1/2)*I*log(tan(x/2) - 1/2 + 3**(1/2)*I/2)/(9*tan(x/2) + 9) + 6/(9*tan(x/2) + 9)
Integral(x*((x - 1)*(x + 2))**(1/2), x)
3**(1/2)*I*log(1/2 - 3**(1/2)*I/2 + exp(-x))/3 - 3**(1/2)*I*log(1/2 + 3**(1/2)*I/2 + exp(-x))/3
exp(x)*sin(2*x)/5 - 2*exp(x)*cos(2*x)/5
x**2*exp(x)*sin(x)/2 + x**2*exp(x)*cos(x)/2 - x*exp(x)*sin(x) + exp(x)*sin(x)/2 - exp(x)*cos(x)/2
-2*x**(5/2)*cos(x**(1/2)) + 40*x**(3/2)*cos(x**(1/2)) - 240*x**(1/2)*cos(x**(1/2)) + 10*x**2*sin(x**(1/2)) - 120*x*sin(x**(1/2)) + 240*sin(x**(1/2))

4.8は近似値を出せって話なのでscipyを利用。warning出てるけどlimit増やせば誤差は減る。これで大丈夫かは解答も載ってないのでわからず。

import scipy as sp
import scipy.integrate
import math

def f1(x):
    return math.sin(8. * x) / math.sin(x)
print sp.integrate.quad(f1, 0, math.pi / 2.)

def f2(x):
    return (math.exp(2 * x) - 1) / (math.exp(x) + 1)
print sp.integrate.quad(f2, 1., 2.)

def f3(x):
    return math.sin(x) ** 2 / x ** 2
print sp.integrate.quad(f3, 0., sp.inf, limit=1000)

def f4(x, a):
    return math.exp(-x) * math.sin(a * x)
print sp.integrate.quad(f4, 0., sp.inf, args=(1.))
print sp.integrate.quad(f4, 0., sp.inf, args=(2.))
(1.4476190476190476, 2.3723635847410075e-11)
(3.6707742704716053, 4.075378113279188e-14)
/Users/tarara/work/py-virtualenvs/numpy-test/lib/python2.7/site-packages/scipy/integrate/quadpack.py:288: UserWarning: The maximum number of subdivisions (1000) 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.
  warnings.warn(msg)
(1.5707758031670074, 1.0326045566833031e-05)
(0.5000000000000002, 1.4875911973447031e-08)
(0.39999999999999836, 1.4432856577266494e-08)