The object of this is to solve the differential equation

{\frac{d^{3}}{d{x}^{3}}}y(x)-\mu\,\left(1-\left(y(x)\right)^{2}\right){\frac{d^{2}}{d{x}^{2}}}y(x)+2\,\mu\, y(x)\left({\frac{d}{dx}}y(x)\right)^{2}+{\frac{d}{dx}}y(x)=0

for the following boundary conditions and parameters:

  • \mu=\frac{1}{2}
  • y\left(0\right)=0
  • {\frac{d}{dx}}y(0)=\frac{1}{2}
  • y\left(2\right)=1

Conventional wisdom would indicate that, because of the high order of the derivatives, this problem cannot be solved using a scalar implementation of simple shooting. However, this is not the case.

To see why this is so, let us begin by implementing the following notation:

\frac{d^{3}}{d{x}^{3}}y(x)=y_{4}\left(x\right)

\frac{d^{2}}{d{x}^{2}}y(x)=y_{3}\left(x\right)

\frac{d}{d{x}}y(x)=y_{2}\left(x\right)

y(x)=y_{1}\left(x\right)

If we substitute these into the differential equation and solve for the third derivative, we have
y_{4}\left(x\right)=\mu\,{\it y_{3}}(x)-\mu\, y_{3}(x)\left(y_{1}(x)\right)^{2}-2\,\mu\, y_{1}(x)\left(y_{2}(x)\right)^{2}-{\it y_{2}}(x)

Now we can construct a series of first-order differential equations for our Runge-Kutta integration scheme as follows:

{\frac{d}{d{x}}}y_{1}(x)=y_{2}\left(x\right)
{\frac{d}{d{x}}}y_{2}(x)=y_{3}\left(x\right)
{\frac{d}{d{x}}}y_{3}(x)=\mu\,{\it y_{3}}(x)-\mu\, y_{3}(x)\left(y_{1}(x)\right)^{2}-2\,\mu\, y_{1}(x)\left(y_{2}(x)\right)^{2}-{\it y_{2}}(x)

Now let us consider things at our first boundary point, namely

y_{4}\left(0\right)=\mu\,{\it y_{3}}(0)-\mu\, y_{3}(0)\left(y_{1}(0)\right)^{2}-2\,\mu\, y_{1}(0)\left(y_{2}(0)\right)^{2}-{\it y_{2}}(0)

Making the appropriate substitutions yields
y_{4}\left(0\right)=1/2\,{\it y_{3}}(0)-1/2

We thus see that, if we use {\it y_{3}}(0) as our independent variable for shooting purposes, we also assume a value of {\it y_{4}}(0) as well. Put another way, since we are given the dependent variable and its first derivative of the ODE at the first boundary, if we guess the second derivative the third derivative automatically follows in spite of the non-linearity of the problem. So we can use a scalar shooting scheme for this problem as well. The dependent variable for root-finding purposes is

F=y_{1}\left(2\right)-y\left(2\right)

which goes to zero as the far boundary condition is met.

To implement this we used a same simple shooting routine with a regula falsi root-finding technique. One thing that was varied was the number of integration steps in the interval; we wanted to see how this affected the convergence.

The results of this are summarized below.

The number of shooting iterations is fairly constant with the variation in integration steps; however, the final value of {\it y_{3}}(0) seems to be refining itself until around 100 integration steps, at which point the accumulation of numerical error begins to affect the precision of the result.

A plot of the results is shown below.

FORTRAN 77 code is below.  Some of the “includes”, subroutines and functions are as follows:

  • matsize.for is a brief snippet of code which includes a parameter statement which sizes the variable size arrays.
  • tstamp calls an OPEN WATCOM routine that returns the date and time and places it in the output.
  • scrsho is a “line plotter” routine, very old school.
c Solution of Third-Order Ordinary Differential Equation
 c Using Simple Shooting Method and Regula Falsi root-finding technique
 c Solution based on Carnahan, Luther and Wilkes (1969)
 include 'matsize.for'
 parameter(neqs=3,nx=10)
 character*40 vinput(nmax)
 character*16 voutpt(nmax)
 c Plot title
 character*60 titlep
 dimension y(neqs),dy(neqs),xyplot(0:nx,2)
 c Define third derivative function
 y4(y1,y2,y3,xmu)=xmu*y3-xmu*y3*y1**2-2.*xmu*y1*y2**2-y2
 c Define second boundary condition function
 bfin(y1,y1fin)=y1-y1fin
 c Enter Boundary Conditions
 data y1strt/0.0/,y2strt/0.5/,y1fin/1.0/xstrt/0.0/,xend/2.0/
 c Enter initial upper and lower bounds for y(3) (Second Derivative)
 data y3left/2.0/,y3rite/0.0/
 c Enter number of shots
 data n/50/
 c Enter friction coefficient
 data xmu/0.5/
 c Enter nomenclature describing initial data as data statements
 data vinput(1)/'Left Initial Bracket'/
 data vinput(2)/'Right Initial Bracket'/
 data vinput(4)/'Delta x increment'/
 data vinput(5)/'First Boundary Condition'/
 data vinput(6)/'Second Boundary Condition'/
 data vinput(7)/'mu'/
 data vinput(8)/'Number of Shots'/
 c Enter nomenclature describing output for individual case
 data vinput(11)/'Shot Number'/
 data vinput(12)/'Value of Guess for Second Derivative'/
 data voutpt(1)/'x'/
 data voutpt(2)/'y(x)'/
 data voutpt(3)/'y`(x)'/
 data voutpt(4)/'y``(x)'/
 data voutpt(5)/'y```(x)'/
 data titlep/'Plot of Final Boundary Condition vs. x'/
 c Description of variables below the comment:
 write(*,*)'Math 5610 Spring 2012 Final Exam'
 write(*,*)'Simple Shooting Method for Third-Order ODE'
 call tstamp
 c Compute value of integration step size for dx
 dx = (xend-xstrt)/float(nx)
 write(*,*)
 c Interval of uncertainty (left)
 1 write(*,*)'Initial Parameters for Problem:'
 write(*,200)vinput(1),y3left
 c Interval of uncertainty (right)
 write(*,200)vinput(2),y3rite
 write(*,200)vinput(4),dx
 write(*,200)vinput(5),xstrt
 write(*,200)vinput(6),xend
 write(*,200)vinput(7),xmu
 write(*,200)vinput(8),n
 do 21 iter=1,n,1
 c ..... Set and print initial conditions
 c ..... Since regula falsi isn't "self starting" it is necessary
 c to compute the left and right values for y(3) (regula
 c falsi independent variable) and compute the corresponding
 c result for bfnact (regula falsi dependent variable)
 c
 c This is done in the first two iterations
 c Iteration 1: Right Value
 c Iteration 2: Left Value
 c Susequent iterations adjust boundaries according to the
 c method
 x = xstrt
 y(1)=y1strt
 y(2)=y2strt
 if(iter.eq.1)then
 y(3) = y3rite
 write(*,*)'Iteration for Right Estimate'
 elseif(iter.eq.2)then
 y(3) = y3left
 write(*,*)'Iteration for Left Estimate'
 else
 y(3) = (y3left*bfinr-y3rite*bfinl)/(bfinr-bfinl)
 endif
 bfnold =bfnact
 y3zero = y(3)
 c Set print index to zero; print index changed at the first, middle
 c and last iteration and printing/plotting takes place
 write(*,*)
 write(*,210)iter
 write(*,200)vinput(1),y3left
 write(*,200)vinput(12),y3zero
 write(*,200)vinput(2),y3rite
 write(*,*)
 write(*,212)(voutpt(np),np=1,5,1)
 write(*,*)
 bfnact=bfin(y(1),y1fin)
 write(*,202)x,y(1),y(2),y(3),bfnact
 c ..... Set Selected variable for plot routine
 xyplot(0,1)=x
 xyplot(0,2)=bfnact
 do 17 icycle=1,nx,1
 c ..... Call Runge-Kutta Subroutine .....
 8 if(irunge(neqs,y,dy,x,dx).ne.1)goto 10
 dy(1)=y(2)
 dy(2)=y(3)
 dy(3)=y4(y(1),y(2),y(3),xmu)
 goto 8
 c ..... Print Solutions, Plot bfnact vs. x values .....
 10 bfnact=bfin(y(1),y1fin)
 write(*,202)x,y(1),y(2),y(3),bfnact
 c ..... Set Selected variable for plot routine
 nplot = icycle
 xyplot(icycle,1)=x
 xyplot(icycle,2)=bfnact
 if(icycle.eq.nx)call scrsho(xyplot,nx,titlep)
 17 continue
 c ..... Ending routine if difference between current and past
 c values of the boundary condition function are small enough
 if(iter.eq.1)then
 bfinr=bfnact
 elseif(iter.eq.2)then
 bfinl=bfnact
 else
 err = abs((bfnold-bfnact)/xmu)
 if(err.lt.1.0e-06)goto 23
 if(bfnact*bfinl.gt.0)then
 y3left = y3zero
 bfinl = bfnact
 else
 y3rite = y3zero
 bfinr = bfnact
 endif
 endif
 21 continue
 23 continue
 200 format(A40,2h ,g15.7)
 202 format(1x,f7.4,2f16.7,2f16.8)
 210 format(35hResults for Regula Falsi Iteration ,i2)
 212 format(1x,5a16)
 stop
 end
function irunge(n,y,f,x,h)
 c Fourth-order Runge-Kutta integration routine
 c Taken from Carnahan, Luther and Wilkes (1969)
 c
 c The function irunge employs the fourth-order Runge-Kutta method
 c with Kutta's coefficients to integrate a system of n simultaneous
 c first-order ordinary differential equations f(j)=dy(j)/dx,
 c (j=1,2,...,n), across one step of length h in the independent
 c variable x, subject to initial conditions y(j), (j=1,2,...,n).
 c Each f(j), the derivative of y(j), must be computed four times
 c per integration step by the calling program. The function must
 c be called five times per step (pass(1)...pass(5)) so that the
 c independent variable value (x) and the solution values (y(1)...y(n))
 c can be updated using the Runge-Kutta algorithm. M is the pass
 c counter. Irunge returns as its value 1 to signal that all derivatives
 c (the f(j)) be evaluated or 0 to signal that the integration process
 c for the current step is finished. Savey(j) is used to save the
 c initial value of y(j) and phi(j) is the increment function for the
 c j(th) equation.
 c
 c The size of the savey and phi arrays depends upon the value of the parameter
 c which is included.
 include 'matsize.for'
 dimension phi(nmax),savey(nmax),y(n),f(n)
 data m/0/
 m=m+1
 goto(1,2,3,4,5), m
 c ..... Pass 1 .....
 1 irunge=1
 return
 c ..... Pass 2 .....
 2 do 22 j=1,n,1
 savey(j)=y(j)
 phi(j)=f(j)
 22 y(j)=savey(j)+h*f(j)/2.0
 x=x+h/2.0
 irunge=1
 return
 c ..... Pass 3 .....
 3 do 33 j=1,n,1
 phi(j)=phi(j)+2.0*f(j)
 33 y(j) = savey(j)+0.5*h*f(j)
 irunge=1
 return
 c ..... Pass 4 .....
 4 do 44 j=1,n,1
 phi(j)=phi(j)+2.0*f(j)
 44 y(j)=savey(j)+h*f(j)
 x=x+0.5*h
 irunge=1
 return
 c ..... Pass 5 .....
 5 do 55 j=1,n,1
 55 y(j)=savey(j)+(phi(j)+f(j))*h/6.0
 m=0
 irunge=0
 return
 end