A 2nd order nonlinear ODE with one boundary and two algebrac equation constraints

Question

How to solve the following nonlinear ODE with two algebraic equations and one boundary condition?

$$y''(x)=\dfrac{2\left((x+15)y'(x)-y(x)\right)\left(y'(x)^2+1\right)}{\left(y(x)^2+x(x+30)+236\right)^2}$$

The boundary condition: $$y(-14)=0$$

The algebraic equation constraint:

$$\left\{ \begin{array}{ll} y(x_0)=\sqrt{1-x_0^2} &\\[15pt] y'(x_0)=\dfrac{-x_0}{\sqrt{1-x_0^2}}& \text{where: }-1\lt x_0\lt 0 \\ \end{array} \right.$$

Answer 1

Make a function that depends on $x_0$ to integrate from $x_0$ to $-14$ using the given formulas as initial conditions and returns the value $y(-14)$. Use the secant method, or some bracketed method to be sure to stay in the interval, to find the value of $x_0$ that gives $y(-14)=0$.

---

The following code gives the tangency to the unit circle at $x_0=-0.0716992989368$, $y(-14)=3.17161661517·10^{-15}$, $y'(-14)=0.0699203372635$. One could replace the fixed step classical Runge-Kutta with one of scipy-odeint or scipy-ode, but that should only change little in the result.

import numpy as np

def odefunc(x,y):
    u,v = y
    # y''(x)=[ 2((x+15)y'(x)-y(x))(y'(x)^2+1) ] / [ (y(x)^2+x(x+30)+236)^2 ]
    return [ v, ( 2*((x+15)*v-u)*(v**2+1) ) / ( (u**2+x*(x+30)+236)**2 ) ];


def odeint(f,t0,y0,tf,h):
    '''Classical RK4 with fixed step size, modify h to fit
    the full interval'''
    N = np.ceil( (tf-t0)/h )
    h = (tf-t0)/N

    t = t0
    y = np.array(y0)
    for k in range(int(N)):
        k1 = h*np.array(f(t      ,y       ))
        k2 = h*np.array(f(t+0.5*h,y+0.5*k1))
        k3 = h*np.array(f(t+0.5*h,y+0.5*k2))
        k4 = h*np.array(f(t+    h,y+    k3))
        y = y + (k1+2*(k2+k3)+k4)/6
        t = t + h
    return t, y

def objective(x0):
    '''integrate a ray tangential to the unit circle 
    to x=-14, return y value'''
    y0 = [ np.sqrt(1-x0**2), -x0/np.sqrt(1-x0**2) ]
    t, y = odeint(odefunc, x0, y0, -14.0, -0.01)
    print x0,t,y
    return y[0]

def illinois(f,a,b):
    '''regula falsi resp. false postion method with
    the Illinois anti-stalling variation'''
    fa = f(a)
    fb = f(b)
    while(np.abs(b-a)>1e-10):
        c = (a*fb-b*fa)/(fb-fa)
        fc = f(c)
        if fa*fc < 0:
            fa *= 0.5
        else:
            a = b; fa = fb
        b = c; fb = fc
    return b, fb

# function value table, inconsequential for solution of the problem
# for k in range(21):
#     objective(-0.99+(0.99*k)/20)
#     
# print "--------------"

x,y = illinois(objective, -0.5,0)
print x,y