-
-
Notifications
You must be signed in to change notification settings - Fork 298
/
compute_fp.jl
73 lines (56 loc) · 1.9 KB
/
compute_fp.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
#=
Compute the fixed point of a given operator T, starting from
specified initial condition v.
@author : Spencer Lyon <spencer.lyon@nyu.edu>
@date: 2014-07-05
References
----------
http://quant-econ.net/jl/dp_intro.html
=#
"""
Repeatedly apply a function to search for a fixed point
Approximates `T^∞ v`, where `T` is an operator (function) and `v` is an initial
guess for the fixed point. Will terminate either when `T^{k+1}(v) - T^k v <
err_tol` or `max_iter` iterations has been exceeded.
Provided that `T` is a contraction mapping or similar, the return value will
be an approximation to the fixed point of `T`.
##### Arguments
* `T`: A function representing the operator `T`
* `v::TV`: The initial condition. An object of type `TV`
* `;err_tol(1e-3)`: Stopping tolerance for iterations
* `;max_iter(50)`: Maximum number of iterations
* `;verbose(true)`: Whether or not to print status updates to the user
* `;print_skip(10)` : if `verbose` is true, how many iterations to apply between
print messages
##### Returns
---
* '::TV': The fixed point of the operator `T`. Has type `TV`
##### Example
```julia
using QuantEcon
T(x, μ) = 4.0 * μ * x * (1.0 - x)
x_star = compute_fixed_point(x->T(x, 0.3), 0.4) # (4μ - 1)/(4μ)
```
"""
function compute_fixed_point{TV}(T::Function, v::TV; err_tol=1e-3,
max_iter=50, verbose=true, print_skip=10)
iterate = 0
err = err_tol + 1
while iterate < max_iter && err > err_tol
new_v = T(v)::TV
iterate += 1
err = Base.maxabs(new_v - v)
if verbose
if iterate % print_skip == 0
println("Compute iterate $iterate with error $err")
end
end
v = new_v
end
if iterate < max_iter && verbose
println("Converged in $iterate steps")
elseif iterate == max_iter
warn("max_iter exceeded in compute_fixed_point")
end
return v
end