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] = y^2 - x^2 list( dy ) } #postscript("../../WriteUp/Graphics/Chapter2/chap_2_sect_7_prob_19_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 19') logistic.trajectory <- trajectory(my_yprime, y0 = c( 0, 0.75 ), t.end = 3.0, parameters = diff_eq_params, col='black') # diverges logistic.trajectory <- trajectory(my_yprime, y0 = c( 0, 0.5 ), t.end = 5.0, parameters = diff_eq_params, col='black') # converges #dev.off() # Thus we know that 0.5 (converges) < alpha_0 < 0.75 (diverges) # Use a bisection like algorithm to reduce the size of the interval. # bkt = c( 0.5, 0.75 ) for( iter in 1:10 ){ alpha_test = mean(bkt) dy.dx = function(x, y){ y^2 - x^2 } h = 0.01 res = euler(dy.dx, h=h, start=0, y0=alpha_test, end=5.0) npts = length(res\$ys) if( res\$ys[npts] < 0 ){ # converges bkt[1] = alpha_test }else{ # diverges bkt[2] = alpha_test } print( bkt ) }