Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
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
2  
@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... –  J. M. Apr 13 '13 at 4:33

2 Answers 2

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$.

share|improve this answer
    
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

share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.