In [1]:
using BasisMatrices
using Plots
using NLsolve
pyplot()
Out[1]:
最初に, 解析的に解が求まるケースを考える.
$MC = c $, $D(p) = p^{-\eta}$ のとき, $S(p) = \eta (p - c) p^{-\eta}$
$\alpha = 1.0, c = 0.2, \eta = 1.5$ とする.
Settings
In [67]:
alpha = 1.0
eta = 1.5
c = 0.5
degree = 25
lower_b = 0.1
upper_b = 5.0
true_S(p) = eta .* (p .- c) .* p .^ (-eta)
marginal_demand(p) = -eta ./ p .^(eta+1)
marginal_cost(q) = q .* 0 .+ c
compute_residual(p, S_p) = (p .+ S_p ./ marginal_demand(p)) .- marginal_cost(S_p)
basis = Basis(LinParams(degree, lower_b, upper_b))
chev_nodes, = nodes(basis)
println(compute_residual([2, 1], [0.1, 0]))
Calculate coefs using NLSolve
broyden method: https://github.com/matthieugomez/NLsolve.jl.git
In [68]:
function objective!(coefs, residuals)
S_p = funeval(coefs, basis, chev_nodes)
residuals .= compute_residual(chev_nodes, S_p)
end
# initial coefficients
init_coefs = zeros(degree)
init_coefs .= 0.1
result = nlsolve(objective!, init_coefs, method = :broyden, xtol = 0.0, ftol = 1.0e-14)
Out[68]:
Chev nodesでの残差は0に近づいているか
In [69]:
coefs = result.zero
residuals_start = compute_residual(chev_nodes, funeval(init_coefs, basis, chev_nodes))
residuals_end = compute_residual(chev_nodes, funeval(coefs, basis, chev_nodes))
for i in 1:degree
@printf "%0.4f, %0.4f \n" residuals_start[i] residuals_end[i]
end
In [70]:
xgrid = collect(linspace(lower_b, upper_b, 101))
f_vals = funeval(coefs, basis, xgrid)
true_vals = true_S(xgrid)
plot(xgrid, f_vals, xlim=(0, 3), ylim=(-2, 3.5), label="Approx S(p)", xlabel = "p", ylabel = "S(p)")
plot!(xgrid, true_vals, xlim=(0, 3), ylim=(-2, 3.5), label="True S(p)", xlabel = "p", ylabel = "S(p)")
Out[70]:
次に問題のケース
In [63]:
alpha = 1.0
eta = 1.5
lower_b = 0.1
upper_b = 3.0
initial_c = 0.5
marginal_demand(p) = -eta ./ p .^(eta+1)
marginal_cost(q) = alpha .* q .^0.5 + q .^2
compute_residual(p, S_p) = (p .+ S_p ./ marginal_demand(p)) .- marginal_cost(S_p)
# 25
degree = 25
basis_25 = Basis(ChebParams(degree, lower_b, upper_b))
chev_nodes_25, = nodes(basis_25)
init_coefs = zeros(degree)
init_coefs .= initial_c
function objective_25!(coefs, residuals)
S_p = funeval(coefs, basis_25, chev_nodes_25)
S_p = abs(S_p)
residuals .= compute_residual(chev_nodes_25, S_p)
end
result_25 = nlsolve(objective_25!, init_coefs, method = :broyden, xtol = 0.0, ftol = 1.0e-12)
coefs_25 = result_25.zero
println(result_25)
# 10
degree = 10
basis_10 = Basis(ChebParams(degree, lower_b, upper_b))
chev_nodes_10, = nodes(basis_10)
init_coefs = zeros(degree)
init_coefs .= initial_c
function objective_10!(coefs, residuals)
S_p = funeval(coefs, basis_10, chev_nodes_10)
S_p = abs(S_p)
residuals .= compute_residual(chev_nodes_10, S_p)
end
result_10 = nlsolve(objective_10!, init_coefs, method = :broyden, xtol = 0.0, ftol = 1.0e-12)
coefs_10 = result_10.zero
println(result_10)
零点で0になっているか?
In [64]:
xgrid = collect(linspace(lower_b, upper_b, 301))
f_vals_25 = compute_residual(xgrid, abs(funeval(coefs_25, basis_25, xgrid)))
plot(xgrid, f_vals_25, xlim=(0, 3), label="Node=25", xlabel = "p", ylabel = "S(p)")
plot!(chev_nodes_25, zeros(25), linestyle = :dot, marker='o')
Out[64]:
In [65]:
xgrid = collect(linspace(lower_b, upper_b, 501))
f_vals_25 = funeval(coefs_25, basis_25, xgrid)
f_vals_10 = funeval(coefs_10, basis_10, xgrid)
plot(xgrid, f_vals_25, xlim=(0, 3), label="Node=25", xlabel = "p", ylabel = "S(p)")
plot!(xgrid, f_vals_10, label="Node=10")
Out[65]:
In [66]:
plot(xgrid, f_vals_25.-f_vals_10, xlim=(0, 3), label="diffs", xlabel = "p", ylabel = "S(p)")
Out[66]:
In [ ]: