Approximating the computation of arctan via a boundary value problem

Question

Suppose

$$ f(x) = \arctan(x) $$

We can write

$$ f^{(2)}(x) = 2x\left( f^{(1)}(x) \right)^2 $$

Subject to some boundary condition $f(a) = c, f(b) = d$. I'm not an expert in this boundary value problems, but I'd like to know what kind of iteration I can use to solve such kind of problem.

I've tried to look at the "shooting method" (I think that's how it's called) but I can't use it because it applies only to linear problems.

Is there anything else I can use?

Even a reference is fine, I'd prefer something which is easy to implement, so I can test it quickly.

Answer 1

user8469759: "I'm not sure I understand how to transform the boundary condition into initial conditions in my case."

Exactly that is the purpose of the shooting method, to transform the second boundary condition into a slope at the first boundary. That is to solve the inverse problem to the equation $\phi(v)=d$ where $\phi(v)=f(b)$ for the solution $f$ to the initial value problem with $f(a)=c$, $f'(a)=v$.

If there were a known inverse function to $\phi$ there would be no problem to solve with the shooting method or any other. In general $\phi$ is a smooth non-linear function, and neither the existence nor the uniqueness of a solution to $\phi(v)=d$ are guaranteed.

As always with non-linear root finders you need some, most likely heuristic, globalization strategy to try a sufficient number of initial slopes to get a global picture of where to look for the correct solution. This can be as simple as picking two random slopes $v_1,v_2$ to start the secant method for $\phi$.

def shoot(x0,y0,xf,yf, eps):
    def derivs(y,x): return [ y[1], 2*x*min(1e10,y[1]**2) ]
    def fb(v0): return odeint(derivs, [y0, v0], [x0, xf], atol=0.01*eps, rtol=0.1*eps)[-1,0]

    v = (yf-y0)/(xf-x0)
    v1 = 0.735*v; fb1 = fb(v1)-yf;
    v2 = 1.234*v; fb2 = fb(v2)-yf;
    while abs(v2-v1)>eps*(1+abs(v2)):
        v3 = v2-fb2*(v2-v1)/(fb2-fb1);
        v1, fb1, v2, fb2 = v2, fb2, v3, fb(v3)-yf
        print "v=%.13f fb(v)=%.13g"%(v2,fb2)
    return v2

x0=0; y0=0; xf=1; yf=0.8; eps=1e-9
v0 = shoot(x0, y0, xf, yf, eps)
A = (v0/(1+v0*x0**2))**0.5
print "exact solution with computed parameters: f(xf)=%.13f"%(A*m.atanh(A*xf))

with the result

v=0.5925363255730 fb(v)=-0.01503501130161
v=0.5954088487471 fb(v)=-0.00959298035706
v=0.6004724102490 fb(v)=7.126956548664e-05
v=0.6004350687242 fb(v)=-3.369897517702e-07
v=0.6004352444583 fb(v)=-1.181577058418e-11
v=0.6004352444645 fb(v)=1.387778780781e-14
exact solution with computed parameters: f(xf)=0.7999999986750

The exact solution is: \begin{align} v'=2xv^2 &\implies 1/v_0-1/v(x)=x^2-x_0^2 \\ f'=v=\dfrac1{1/v_0+x_0^2-x^2} &\implies f(x)=A\cdot \operatorname{Artanh}(Ax) \text{ where } A=\sqrt{\frac{v_0}{1+v_0x_0^2}} \end{align} so that with the computed values $A=1/\sqrt{x_0^2+1/v_0}=\sqrt{v_0}=\sqrt{0.6004352444645}=0.77487756740305624$ and indeed $$ A\cdot{\rm Artanh}(A\cdot 1)=0.79999999867501892 $$ gives the second boundary condition within the error level.

You might have remarked that the solution is not the arcus tangent, as that has the differential equation $f''=-2xf'^2$.