fprime = function (f,a,h=0.0001){(f(a+h)-f(a-h))/(2*h)}
bisection = function(f,a,b,tol=0.0001){
if (f(a)*f(b) > 0){
return ("Boundary Conditions Not Met")
}
else{
middle = a
while (abs(f(middle))>tol){
middle = (a+b)/2
if (f(middle)*f(a)>0) (a = middle)
else (b = middle)
x=middle
y=f(middle)
## if you want to "see" what happens at every step, take off the # of the next line ##
#cat(sprintf("x-Val: %.4f ; f(x-val): %.4f\n",x,y))
}
return (middle)
}
}
## Solve this problem ##
x = seq(0,20)
profit = function(t){(0.65-0.01*t)*(200+5*t)-0.45*t}
plot(x,profit(x),"l")
dProfit = function(t){fprime(profit,t)}
ans = bisection(dProfit,0,20)
## Show what seq does ##
r = seq(.008,.012,.001)
## Show how to do a for loop ##
for (r in seq(.008,.012,.001)){print(r)}
## [1] 0.008
## [1] 0.009
## [1] 0.01
## [1] 0.011
## [1] 0.012
## Solve all of the problems ##
ans.time = 0
ans.prof=0
r = seq(.008,.012,.001)
for (i in 1:length(r)){
profit = function(t){(0.65-r[i]*t)*(200+5*t)-0.45*t}
dProfit = function(t){fprime(profit,t)}
ans.time[i] = bisection(dProfit,0,20)
ans.prof[i] = profit(ans.time[i])
}
print(ans.time) ## See book##
## [1] 15.000000 11.110840 8.000488 5.454102 3.332520
print(ans.prof)
## [1] 139.0000 135.5556 133.2000 131.6364 130.6667
result = data.frame(price=r,time=ans.time,profit=ans.prof)
print(result) ## Table 1.1 in book
## price time profit
## 1 0.008 15.000000 139.0000
## 2 0.009 11.110840 135.5556
## 3 0.010 8.000488 133.2000
## 4 0.011 5.454102 131.6364
## 5 0.012 3.332520 130.6667
## Investigate other times/situations not in the book!! ##
ans.time = 0
ans.prof=0
r = seq(.008,.02,.001)
for (i in 1:length(r)){
profit = function(t){(0.65-r[i]*t)*(200+5*t)-0.45*t}
dProfit = function(t){fprime(profit,t)}
ans.time[i] = bisection(dProfit,-100,20)
ans.prof[i] = profit(ans.time[i])
}
print(ans.time) ## See book##
## [1] 15.0012207031 11.1120605469 7.9992675781 5.4541015625 3.3337402344
## [6] 1.5383911133 -0.0006103516 -1.3336181641 -2.5000000000 -3.5290527344
## [11] -4.4445800781 -5.2630615234 -6.0000610352
print(ans.prof)
## [1] 139.0000 135.5556 133.2000 131.6364 130.6667 130.1538 130.0000 130.1333
## [9] 130.5000 131.0588 131.7778 132.6316 133.6000
result = data.frame(price=r,time=ans.time,profit=ans.prof)
print(result) ## Table 1.1 in book
## price time profit
## 1 0.008 15.0012207031 139.0000
## 2 0.009 11.1120605469 135.5556
## 3 0.010 7.9992675781 133.2000
## 4 0.011 5.4541015625 131.6364
## 5 0.012 3.3337402344 130.6667
## 6 0.013 1.5383911133 130.1538
## 7 0.014 -0.0006103516 130.0000
## 8 0.015 -1.3336181641 130.1333
## 9 0.016 -2.5000000000 130.5000
## 10 0.017 -3.5290527344 131.0588
## 11 0.018 -4.4445800781 131.7778
## 12 0.019 -5.2630615234 132.6316
## 13 0.020 -6.0000610352 133.6000