if( !require('phaseR') ){ install.packages('phaseR') } library(phaseR) if( !require('FNN') ){ install.packages("FNN", dependencies=TRUE, repos='http://cran.rstudio.com/') } library(FNN) source('euler.R') my_yprime = function(t, y, parameters){ x = y[1] y = y[2] dy = numeric(2) dy[1] = 1 dy[2] = -x*y + 0.1*y^3 list( dy ) } #postscript("../../WriteUp/Graphics/Chapter2/chap_2_sect_7_prob_18_plot.eps", onefile=FALSE, horizontal=FALSE) diff_eq_params = c() logistic.flowField <- flowField(my_yprime, x.lim=c(0, 5), y.lim=c(-5, 5), parameters=diff_eq_params, points=21, add=FALSE, main='Problem 18') logistic.trajectory <- trajectory(my_yprime, y0 = c( 0, 2.5 ), t.end = 3.0, parameters = diff_eq_params, col='black') # diverges logistic.trajectory <- trajectory(my_yprime, y0 = c( 0, 2.0 ), t.end = 3.0, parameters = diff_eq_params, col='black') # converges #dev.off() # Thus we know that 2.0 (converges) < alpha_0 < 2.5 (diverges) # Use a bisection like algorithm to reduce the size of the interval. # bkt = c( 2.0, 2.5 ) for( iter in 1:10 ){ alpha_test = mean(bkt) dy.dx = function(x, y){ -x*y + 0.1*y^3 } h = 0.01 res = euler(dy.dx, h=h, start=0, y0=alpha_test, end=3.0) ys = res\$ys ys = ys[is.finite(ys)] npts = length(ys) if( ys[npts] < 0.5 ){ # converges bkt[1] = alpha_test }else{ # diverges bkt[2] = alpha_test } print( bkt ) }