Solving coupled differential equations with an eigenvalue

Question

I have 2 coupled differential equations with an eigenvalue Ei and want to solve them

 ϕ''[r] + (2/r) ϕ'[r] - mb^2 ϕ[r] + (Ei + g*A[r])^2 ϕ[r] == 0
 A''[r] + (2/r) A'[r] - mv^2 A[r] - 2 g (Ei + g*A[r]) (ϕ[r])^2 == 0

where mb, mv and g are constants equal to 1. The boundary conditions of these equations are

ϕ[0] = 1, ϕ'[0] = 0, A[0] = 0, A'[0] = 0

Because of the singularity of r, we assume r = 1*10^-8. What I want to find is the maximum radius (r) when ϕ[rmax] = 0, ϕ'[rmax] = 0, A[rmax] = 0, A'[rmax] = 0. I've with this

mb = mv = g = 1; 
b = ParametricNDSolveValue[{ϕ''[r] + (2/r) ϕ'[r] - mb^2 ϕ[r] + (Ei + g A[r])^2 ϕ[r] == 0,
          A''[r] + (2/r) A'[r] - mv^2 A[r] - 2 g (Ei + g A[r]) (ϕ[r])^2 == 0,
          ϕ[0.00000001] == 1, ϕ'[0.00000001] == 0, A[0.00000001] == 0,
          A'[0.00000001] == 0}, {ϕ, A},
     {r, 0.00000001, 100}, {Ei}]

but it didn't work when I wanted to find Ei like this

val = Map[FindRoot[b[Ei][100], {Ei, #}] &, {1, 2, 3}]

and the maximum radius. How should I solve this issue ? Also I want to plot ϕ and A vs r.

What approach do you want to take for finding Ei? b[Ei] actually returns 2 InterpolatingFunction while you only have 1 unknown! If you replace the {ϕ, A} inside ParametricNDSolve with ϕ or A, you'll get a solution, whether the solution is correct or not is another story. BTW, you may want to try this package: http://library.wolfram.com/infocenter/Demos/4482/ — xzczd, Apr 18 '15 at 04:11
I want to find Ei where the ϕ and the ϕ' when r (maximum radius) are equal to 0. So I think I can use findroot. — M. Fitrah Alfian R. S., Apr 21 '15 at 07:53
As mentioned above, ϕ[Ei] == 0 and ϕ'[Ei] == 0 give two equations, use either of them inside FindRoot or RootSearch will give you a result. — xzczd, Apr 21 '15 at 07:56
How do I have to write the code of FindRoot like what you've said in Mathematica? — M. Fitrah Alfian R. S., Apr 21 '15 at 09:02
I noticed that you've modified your equations for several times, are you sure it's correct this time? Maybe you can show us the original problem? — xzczd, Apr 23 '15 at 03:40
Yes, I think. I just change the sign of some parts of the equations. The original problem is I want to solve those 2 differential equations where there is an eigenvalue that I have to get. The boundary conditions are still the same. Till now, I don't get the eigenvalue because there will always an error when I want to find eigenvalue by using FindRoot. Maybe you can tell me the other way how to find eigenvalue like I want. Thanks — M. Fitrah Alfian R. S., Apr 26 '15 at 15:26
Seems that this hasn't been done: Hello, welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the [faq]! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! — xzczd, Apr 28 '15 at 10:14

score 1 · Accepted Answer · edited Jun 16 '20 at 09:23

OK, let me explain why your attempt doesn't succeed. Here I'll modify your sample a little for clarity (Of course the troublesome part is leaved untouched ):

lb = 10^-8;
mb = mv = g = 1;
eqn = {ϕ''[r] + (2/r) ϕ'[r] - mb^2 ϕ[r] + (Ei + gA[r])^2 ϕ[r] == 0,
       A''[r] + (2/r) A'[r] - mv^2 A[r] - 2 g (Ei + gA[r]) (ϕ[r])^2 == 0};
bc = {ϕ[lb] == 1, ϕ'[lb] == 0, A[lb] == 0, A'[lb] == 0};
b = ParametricNDSolveValue[{eqn, bc}, {ϕ, A}, {r, lb, 100}, {Ei}]

enter image description here

As you see, ParametricNDSolveValue returns a ParametricFunction, and again when you give a numeric value to the ParametricFunction, it will return something: it will return what? Let's have a try:

pb = b[1]

enter image description here

Hmm… the calculation stopped at about t = 8.95, far below 100, but it's not the key point here, the key point is, b[1] doesn't return a InterpolatingFunction, it returns a List of InterpolatingFunction! Then what will something like pb[1] return? A expression that FindRoot doesn't understand:

pb[1]

enter image description here

There're many way to circumvent this problem, let me show several here. Since you only have one unknown Ei, only one equation is needed. I'll use ϕ[Ei] for example:

(* Circumvention 1 *)
fb[x_?NumericQ] := b[x][[1]]
FindRoot[fb[Ei][5], {Ei, 1}]
(* Circumvention 2 *)
fb2[x_?NumericQ] := b[x][5] // Through // First
FindRoot[fb2[Ei], {Ei, 1}]
(* Circumvention 3 *)
newb = ParametricNDSolveValue[{eqn, bc}, ϕ, {r, lb, 100}, {Ei}]
FindRoot[newb[Ei][5], {Ei, 1}]
(* Circumvention 4 *)
newb2 = ParametricNDSolve[{eqn, bc}, {ϕ, A}, {r, lb, 100}, {Ei}]
fϕ = ϕ /. newb2
FindRoot[fϕ[Ei][5], {Ei, 1}]

You'll still see some warnings when you execute the above code, but it's mainly because of the nature of your equation and it's another issue, at least FindRoot works this time!

Well, personally I feel the behavior of ParametricNDSolveValue undesirable, why it doesn't return a List of ParametricFunction?

Finally, my intuition told me that you may be interested in this post.

It's very helpful. Thank you for your help. – M. Fitrah Alfian R. S. Apr 28 '15 at 09:18 — M. Fitrah Alfian R. S., Apr 28 '15 at 09:18

Solving coupled differential equations with an eigenvalue

1 Answers1