# This file is part of the Test Set for BVP solvers # # http://www.dm.uniba.it/~bvpsolvers/ # # Problem Flat_Moon # ODE of dimension 7 # # DISCLAIMER: see # http://www.dm.uniba.it/~bvpsolvers/testsetbvpsolvers # # Authors: # # Jeff R. Cash # (Department of Mathematics, Imperial College, London, England.) # Davy Hollevoet # (Vakgroep Toegepaste Wiskunde en Informatica, Universiteit Gent, Belgium.) # Francesca Mazzia # (Dipartimento di Matematica, Universita' di Bari, Italy) # Abdelhameed Nagy Abdo # (Dipartimento di Matematica, Universita' di Bari, Italy) # (Dept. of Mathematics, Faculty of Sciences, Benha University,Egypt) Flat_Moon<- function(){ prob<- function(){ fullnm <- 'Problem Flat_Moon' problm = 'Flat_Moon' typebvp = 'BVP' neqn = 7 nlbc = 4 aleft = 0 aright = 1 numjac = TRUE numbcjac = TRUE linear = FALSE Rpar<-c(0.5) Ipar<-0 return(list(fullnm=fullnm,problm=problm,typebvp=typebvp,neqn=neqn, nlbc=nlbc,aleft=aleft,aright=aright,Rpar=Rpar,Ipar=Ipar, numjac=numjac,numbcjac=numbcjac,linear=linear)) } init <- function(neqn) { givmsh = TRUE givey = TRUE nmsh = 41 yval <- rep(0,neqn) yval[neqn]=700 xguess = seq(0,1,by=1/(nmsh-1)) yguess <- NULL for (i in 1:neqn) yguess <- rbind(yguess,rep(yval[i],nmsh)) return(list(givmsh=givmsh,givey=givey, nmsh=nmsh,xguess=xguess,yguess=yguess)) } feval = function(x,y,eps,Rpar,Ipar){ g<-1.62e0 A<-3e0*1.62e0 TF<-y[7] f <- c(y[3]*TF, y[4]*TF, A*cos(y[5])*TF, (A*sin(y[5])-g )*TF, (-y[6]*cos(y[5]))*TF, (((y[6])^2)*sin(y[5]))*TF, 0 ) return(list(f)) } jeval = function(x,y,eps,Rpar,Ipar){ dfy <- matrix( nrow=6,ncol=6,byrow = TRUE, data = c( 0, 1.0e0, 0, 0, 0, 0, y[4]/eps, -y[3]/eps , -y[2]/eps, y[1]/eps, 0, 0, 0, 0, 0,1.0e0, 0, 0, 0, 0, 0, 0, 1.0e0, 0, 0, 0, 0, 0, 0, 1.0e0, -y[2]/eps, -y[1]/eps, -y[6]/eps, 0, 0, -y[3]/eps )) return((dfy)) } bcfun <- function(i, y, eps,Rpar,Ipar) { VC=1627e0 h=185.2e0 if (i == 1) return(y[1]) if (i == 2) return(y[2]) if (i == 3) return(y[3]) if (i == 4) return(y[4]) if (i == 5) return(y[2]-h) if (i == 6) return(y[3]-VC) if (i == 7) return(y[4]) } dbcfun <- function(i, y, eps,Rpar,Ipar) { if (i == 1) return(c(1, 0)) if (i == 2) return(c(1, 0)) if (i == 1) return(c(1, 0)) if (i == 2) return(c(1, 0)) } setoutput<- function(neqn,plotsol=NULL){ solref = FALSE if (is.null(plotsol)){ nindsol = neqn indsol = 1;neqn } else{ nindsol = length(plotsol) indsol = plotsol } } esolu <- function(X,parms,Rpar,Ipar){ lambda=parms Exact =matrix( c( exp(-X/(lambda)), (-exp(-X/(lambda)))/(lambda)), ncol=2,nrow=length(X),byrow=FALSE) return(Exact=Exact) } return(list(prob=prob,init=init,feval=feval,jeval=jeval,bcfun=bcfun,dbcfun=dbcfun,setoutput=setoutput,esolu=esolu)) }