数値計算には気をつけよ

某所で，以下の数学入試問題を Python で解くというのをやっていたけど，ちょっと雑な扱いだったのでちゃんとやっておく。

元記事では f(x) のグラフを描くだけでお茶を濁しているけど，どうせなら sympy を使えば良いかと思います。

Julia でも SymPy で使えます。

(1) は導関数を求めて x > 0 で 負であることを示せば良いかな。グラフを描くなら導関数のグラフを描こう。外観がつかめたら lim x → 0，lim x  → ∞ でどうなるか limit で示す。

using SymPy
@syms a b f g k w x z

f = (x + 1) * log(1 + 1/x)
g = diff(f) # 導関数

limit(g, x, 0) # -∞
limit(g, x, Inf) # 0

(2) でも，「解く方法としてはいたって簡単です。代入していくだけです。脳筋です。」として，図を描こうとしているのだが，

pow(k+1,k+1)/special.factorial(k)/pow(a,k)

をそのまま使ったものだから下のような図になりまして，

「10~20の間で数字が落ちているようです。今回は最大値を求めろなので問題なさそうです。どうやらk=2,3のときに最大値を取るようです。」としているが，何故「数字が落ちている」のか考察しようともしていないのはちょっとねえ。

exp((k+1)*log(k + 1) - k*log(2^8/ 3^4)- logfactorial(k))

でやらないとダメです。比較するとどうしてダメだったのかわかるでしょう。

julia> using SpecialFunctions, Printf

julia> funcb(k) = (k + 1) ^ (k + 1) / (2^8/ 3^4)^k / factorial(k)
funcb (generic function with 1 method)

julia> funcb2(k) = exp((k+1)*log(k + 1) - k*log(2^8/ 3^4)- logfactorial(k))
funcb2 (generic function with 1 method)

julia> for k = 0:20
               @printf("%2d %20.15f %20.15f\n", k, funcb(k), funcb2(k))
       end

 0    1.000000000000000    1.000000000000000
 1    1.265625000000000    1.265625000000000
 2    1.351524353027344    1.351524353027344
 3    1.351524353027344    1.351524353027344
 4    1.305025489273248    1.305025489273248
 5    1.232967201857537    1.232967201857539
 6    1.147688698458480    1.147688698458479
 7    1.056829030163603    1.056829030163603
 8    0.965211505305653    0.965211505305653
 9    0.875875642792911    0.875875642792914
10    0.790691444309025    0.790691444309024
11    0.710746653906890    0.710746653906893
12    0.636600915295705    0.636600915295703
13    0.568457341240705    0.568457341240706
14    0.506280597184257    0.506280597184257
15    0.000000000000000    0.449879200714691
16   -0.001380880825857    0.398963280942955
17   -0.000004461526006    0.353185192585937
18    0.000001002536991    0.312167985688245
19   -0.000000005521767    0.275525189547352
20   -0.000000000066342    0.242874345836610

なお，「1.351524353027344 が最大値Mのようです(既約分数ではありませんが...)」とあるが，これもJulia の rational で計算すると
(3^3 // (2^8 // 3^4)^2 // 2)　# 177147//131072
(4^4 // (2^8 // 3^4)^3 // 2 // 3)　# 177147//131072

funcb (generic function with 1 method)

z(2)

1.351524353027344

z(3)

1.351524353027344

 

(3^3 // (2^8 // 3^4)^2 // 2)

177147//131072

float(ans)

1.3515243530273438

(4^4 // (2^8 // 3^4)^3 // 2 // 3)

177147//131072

747/2

373.5