4

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.

Michael
  • 63
  • 1
  • 5
  • 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 03: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 04: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 04: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 04:31
  • In any case, OP might consider posting his actual problem instead... – J. M.'s missing motivation Apr 13 '13 at 04:33
  • @Jens I don't know what you mean by "such" equations, in general there is no consistent way you can turn a BC to an equivalent IC( without knowing the general solution of course). If by such equations you means equations that have cosines or sines as solutions then in a sense you know the general solution. In general a BC problem might not have a solution and I believe you know that the shooting method employed for solving such a problem numerically starts by a guess of the derivative at the initial point and not knowing the value of the derivative that fits the BC . – Spawn1701D Apr 13 '13 at 04:40
  • @Spawn1701D I'm just taking the question as stated: "an analogous problem is..." so we have a linear second order BVP. Then of course you can't use the initial derivative as the parameter. It's simply an eigenvalue problem, and you shoot for the value of $\omega$ that satisfies the BC. My guess is that you're reading too much into the question. – Jens Apr 13 '13 at 04:56
  • @Jens thats why I ask Paul about more info about his problem. Your answer, according to his given information, is perfectly sound. – Spawn1701D Apr 13 '13 at 05:04
  • I have updated the original post with the full problem. Hopefully this helps to clarify things. As mentioned in the post I still seem to have a problem. – Michael Apr 13 '13 at 09:53
  • Do not use K as a function name. It is a symbol used internally by Mathematica (note that it's coloured black in the front end), and will it cause problems. Better: don't use any names starting with capital letters, then you can avoid conflict easily. – Szabolcs Apr 13 '13 at 15:18
  • The edit contains various syntax errors. Definitions like f[u] should be f[u_], spaces are missing between names. – Jens Apr 13 '13 at 18:12
  • Sorry, I believe I have fixed all the syntax errors in my post. These were errors in my posting of the solution though, in my mathematica file these problems were not there and I still have the same error message. – Michael Apr 13 '13 at 19:46

2 Answers2

5

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

Jens
  • 97,245
  • 7
  • 213
  • 499
  • 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 09: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
3

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

Szabolcs
  • 234,956
  • 30
  • 623
  • 1,263
  • 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