# How can I use FindRoot on an expression from NDSolve?

I have a second order ODE that I can only solve numerically using NDSolve, but I then need to use the solution in FindRoot and am running into errors. A simplified but analogous problem is the following: Find the solution to $\phi^{\prime\prime}(u) = -\omega^2 \phi(u)$ on the interval $0<u<1$ with the boundary conditions $\phi(0)=\phi(1)=0$, which will only be true for certain values of $\omega$ ($\omega=n\pi$ in this case). The code I am using is

eqnp = p''[ u] + ω^2 p[u];

psol[ω_?NumericQ] = NDSolve[{eqnp == 0, p[0] == 0, p[1] == 0}, p, u];

FindRoot[psol[ω], {ω, 3}]


but I am receiving errors. How can I feed a solution from NDSolve into FindRoot? I know that in the specific case above everything can be done analytically, but my actual ODE must be solved numerically, so I need to figure out how to solve this simplified problem completely numerically.

EDIT: I have updated the post with the full problem. The ODE is both homogeneous and linear

λ = 0.00001; k = 1.0;
f[u_] := 1/(2λ u) (1-Sqrt[1-4λ+4λ u^2]);
A = Sqrt[1/2(1+Sqrt[1-4λ])];
K2[u_] := u^(-3/2)f[u](1 + 2λ u^2f'[u]);
K1[u_] := K2[u](4ω^2)/(A^2f[u]^2) - 4k^2u^(-1/2)(1 - λ(6u^2f'[u] + 4 u^3f''[u]));
eomu = 4u^3K2[u]ϕ''[u] + (6u^2K2[u] + 4u^3K2'[u])ϕ'[u] + K1[u]ϕ[u];


I define the boundary conditions at $u=1$ by finding the first few coefficients in the series solution (using normal Froebenius method). I then want to find the solution and use Dirichlet boundary conditions at $u=0$ to determine the eigenfrequencies. The first coefficient is $a1$ and the indicial exponent is $-\frac{I\omega}{2A}$ (the other solution to the indicial equation is unphysical and is discarded)

a1 = (k^2(1 + Sqrt[1 - 4λ])(1 + 8λ) + 2(-1 + 2λ)ω^2 + (I Sqrt[1 + Sqrt[1 - 4λ]] (12λ ω))/Sqrt[2])/(2Sqrt[2] Sqrt[1 + Sqrt[1 - 4λ]] (Sqrt[1 + Sqrt[1 - 4λ]]/Sqrt[2] - I ω))


(I have separate code to find an arbitrary number of coefficients, but to enter them here would be too much) Then define the solution to use as boundary conditions

ϕhnum[u_] := (1 - u)^(-I ω/(2A))(1 + a1(1 - u))
numh = 0.999999;
numb = 0.000001;


Then I use the code suggested by Jens

ϕsol[ω_?NumericQ] := ϕ /. First@NDSolve[{eomu == 0, ϕ[numh] == ϕhnum[numh], ϕ'[
numh] == ϕhnum'[numh]}, ϕ, {u, numb, numh}];

FindRoot[ϕsol[ω][numb], {ω, 1.0 - 1.0 I}]


When I do this I get the following error

FindRoot::lstol: The line search decreased the step size to within tolerance specified
by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the
merit function. You may need more than MachinePrecision digits of working precision to
meet these tolerances. >>


I know that for the first solution I expect something like $\omega = 1.95- 1.27 I$. Hopefully this clarifies the problem, and thanks again for any help.

-
You have to give a non trivial example because this one has the trivial solution ($u=0$) when $\omega\ne\kappa\pi$. –  Spawn1701D Apr 13 '13 at 3:58
@Spawn1701D The fact that his example has a trivial solution doesn't mean that the example is trivial. It has a non-trivial solution too (any linear homogeneous equation has a trivial solution, but nobody cares). How to obtain the latter using NDSolve is a valid question, which I answered. –  Jens Apr 13 '13 at 4:11
@Jens yes you answered but you used the knowledge of the solution itself, i.e. that $p'[0]==0$. I believe the problem he wants to attack: a) doesn't know the general solution b) his problem probably is a nonlinear or/and nonhomogeneous one. No one questions the validity of your answer and probably with some slight modifications it can be adapted to his problem I was just curious of his proper problem. –  Spawn1701D Apr 13 '13 at 4:26
@Spawn1701D The initial condition is $p'[0]=1$, and that's a standard choice which doesn't require knowing the solution. It's just selecting the non-trivial solution consistent with the boundary condition at $0$. The other linearly independent solution is $p[0]=1$, $p'[0]=0$; this is always true for such equations. Basically what I'm doing is a shooting method that converts the BVP into an initial value problem. I'm sure that's what Paul is doing, and that's why FindRoot comes in. –  Jens Apr 13 '13 at 4:31
In any case, OP might consider posting his actual problem instead... –  Ｊ. Ｍ. Apr 13 '13 at 4:33

There are several things that need to be corrected. First of all, to get a non-trivial solution from NDSolve you have to

1. specify the start and end of the interval for u
2. Specify two initial conditions for a second-order differential equation. Initial and final conditions (corresponding to your boundary-value problem) won't work.

Since you want Dirichlet boundary conditions, you can choose the second initial condition to be a non-vanishing derivative. With this, you'll get a solution that doesn't satisfy the Dirichlet boundary condition p[1] == 0.

Then the job of FindRoot will be to adjust $\omega$ until the latter condition is saitsfied:

eqnp = p''[u] + ω^2 p[u];

psol[ω_?NumericQ] := p /. First@NDSolve[{
eqnp == 0,
p[0] == 0,
p'[0] == 1
},
p, {u, 0, 1}];

FindRoot[psol[ω][1], {ω, 3}]

(* ==> {ω -> 3.14159} *)


To make FindRoot work with the solution, I also had to convert the output of NDSolve to an actual function first. The output is a rule, of which I take the First part which then replaces the generic p. So in FindRoot, the argument psol[ω][1] is the solution function for a given $\omega$, evaluated at the right-hand boundary u=1. When that's zero, you've found the correct value of $\omega$.

-
Thank you for your help with this problem, and for the clear explanation. I have updated my original post with the full problem. When I used the provided code, I still find an error, and I am not sure where it is from. Perhaps it has to do with the boundary conditions depending on $\omega$? –  Paul Apr 13 '13 at 9:27
I can't execute your edited new code, but the fact that you get the quoted error message means that you probably did everything right up to the point of finding $\omega$. The issue you're seeing now is unrelated to the original one. There's a lack of precision in calculating the NDSolve solution, or maybe the zero is not where you think it is. –  Jens Apr 13 '13 at 18:07
Thanks for the help. I was thinking that the issue is now different so it is nice to confirm that. However I am pretty sure that the $\omega$ is near where I am searching, since I have confirmed this value in an independent way. I will see what I can find about my current error message (which is the one given at the bottom of the original post) –  Paul Apr 13 '13 at 19:57

Here's a more concise version of Jens's solution, using the new NDSolve-related features introduced in version 9:

fun =
ParametricNDSolveValue[{p''[u] + ω^2 p[u] == 0, p[0] == 0, p'[0] == 1},
p[1], {u, 0, 1}, ω]

FindRoot[fun[ω], {ω, 1}]

(* ==> {ω -> 3.14159} *)


This will find ω so that p[1] == 0.

Taking the fixed equation from your last update, and applying the same method works:

parfun = ParametricNDSolveValue[{eomu ==
0, \[Phi][numh] == \[Phi]hnum[numh], \[Phi]'[numh] == \[Phi]hnum'[
numh]}, \[Phi], {u, numb, numh}, \[Omega]]


Attempting FindRoot gives the error you encounter and gives no correct solution. There's a good chance that there's no solution: there's no ω so that parfun[ω][numb] == 0

-
I was working in version 8 when I wrote my answer so I couldn't go this route, but I agree this is a nice new feature (+1). –  Jens Apr 13 '13 at 18:10
Maybe you have time to try the OP's new original code. I can't do it right now. –  Jens Apr 13 '13 at 18:13