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 ϕ′′(u)=−ω2ϕ(u)\phi^{\prime\prime}(u) = -\omega^2 \phi(u) on the interval 0>

I know that for the first solution I expect something like ω=1.95−1.27I\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=0u=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]==0p'[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]=1p'[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 00. The other linearly independent solution is p[0]=1p[0]=1, p′[0]=0p'[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…
– J. M.♦
Apr 13 ’13 at 4:33


2 Answers


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

specify the start and end of the interval for u
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?
– Michael
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)
– Michael
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