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

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

問題3.16-3.18の解答(SymPy)

概説微分積分の問題3.16から3.18までの解答。

ソフト使って解く問題だったのでSymPyを使ってみた。機能を探すのにかなり手間取ったけど、数式を代数的に扱えて面白い。∞はooだったりする。SymPyとnumpyとscipyとの関連性はまだよくわからず。

ただしPlotはpygletがまともに動かなくて駄目だった。この辺はMatplotlibに任せるべきだったかもしれないけど、力尽きた。

import sympy as sym
from sympy.abc import x, y
from sympy.geometry import Point, Line
from sympy.utilities.autowrap import autowrap

print "[3.16 (1)]"
f1 = (x ** 2 - 1) * (x + 2) * (x - 3) * sym.sin(2 * sym.pi * x)
df1 = sym.diff(f1, x)
print df1.subs(x, 2.)
print Line(Point(2., f1.subs(x, 2)), slope=df1.subs(x, 2.)).equation(), "= 0"

print "[3.16 (2)]"
f2 = (sym.sqrt(1 + x ** 2) + sym.sqrt(1 - x ** 2)) / \
     (sym.sqrt(1 + x ** 2) - sym.sqrt(1 - x ** 2))
df2 = sym.diff(f2, x)
print df2

print "[3.17 (1)]"
f3 = sym.sin(x) - 3 * sym.sin(2 * x) + 2 * sym.cos(3 * x)
# sym.Plot(f3, [x, 0, 2 * sym.pi, 100])
df3 = sym.diff(f3, x)
df3_callable = autowrap(df3)
print sym.mpmath.findroot(df3_callable, (0, 1.5), solver='anderson')
print sym.mpmath.findroot(df3_callable, (1.5, 3), solver='anderson')
print sym.mpmath.findroot(df3_callable, (3, 4), solver='anderson')
print sym.mpmath.findroot(df3_callable, (4, 4.7), solver='anderson')
print sym.mpmath.findroot(df3_callable, (4.7, 5), solver='anderson')
print sym.mpmath.findroot(df3_callable, (5, 6.5), solver='anderson')

print "[3.17 (2)]"
f4 = 3 * sym.exp(- x ** 2) \
     - 2 * sym.exp(- 2 * (x - 1) ** 2) - 2 * sym.exp(- 2 * (x + 1) ** 2) \
     + sym.exp(- 3 * (x - 2) ** 2) + sym.exp(- 3 * (x + 2) ** 2)
# sym.Plot(f4, [x, -3, 3])
df4 = sym.diff(f4, x)
df4_callable = autowrap(df4)
print sym.mpmath.findroot(df4_callable, (-3, -2), solver='anderson')
print sym.mpmath.findroot(df4_callable, (-2, -1), solver='anderson')
print sym.mpmath.findroot(df4_callable, (-1, 1), solver='anderson')
print sym.mpmath.findroot(df4_callable, (1, 2), solver='anderson')
print sym.mpmath.findroot(df4_callable, (2, 3), solver='anderson')

print "[3.18 (1)]"
print sym.limit(x ** (1. / x), x, sym.oo)
print "[3.18 (2)]"
print sym.limit((x - sym.sin(x)) / (sym.atan(x) - x) , x, 0)
[3.16 (1)]
-24.0*pi
24.0*pi*x + 1.0*y - 48.0*pi = 0
[3.16 (2)]
(-x/(x**2 + 1)**(1/2) - x/(-x**2 + 1)**(1/2))*((-x**2 + 1)**(1/2) + (x**2 + 1)**(1/2))/(-(-x**2 + 1)**(1/2) + (x**2 + 1)**(1/2))**2 + (x/(x**2 + 1)**(1/2) - x/(-x**2 + 1)**(1/2))/(-(-x**2 + 1)**(1/2) + (x**2 + 1)**(1/2))
[3.17 (1)]
0.921193032136374
2.17913969064355
3.51092140472367
4.645447617455
4.71238898038469
6.02205784978526
[3.17 (2)]
-2.10620408177109
-1.15895914580887
0.0
1.15895914580887
2.10620408177109
[3.18 (1)]
1
[3.18 (2)]
-1/2