# Solving system of equations with (nested) Bessel Functions

I'm trying to find solutions to a system of equations which include Bessel functions. It feels like the problem for mathematica is that the arguments of the Bessel functions are functions themselves but I don't see, why mathematica can't get a numeric solution for this.

My Code is

 ν = 10; c = 1/2; α = 0.7 ; β = 0.55  ;
vA[x_] := -c + x;
vB[x_] := c + x;
F[x_] := 1/2 + x/2;

σA[α_,xaA_,xbA_]:=αF[xaA] + (1-α)F[xbA];
σB[α_,xaB_,xbB_]:=α(1 - F[xaB]) + (1-α)(1 - F[xbB]);
τA[β_, xaA_, xbA_] := (1 - β) F[xaA] + β F[xbA];
τB[β_, xaB_, xbB_] := (1 - β) (1 - F[xaB]) + β (1 - F[xbB]);

pAA[ν_Integer, α_, xaA_, xbA_, xaB_, xbB_] :=
Exp[-ν (σA[α, xaA, xbA] + σB[α,
xaB, xbB] )] Sum[(ν σA[α, xaA, xbA])^
k  (ν σB[α, xaB, xbB])^(k +
1)/(k! (k + 1)!), {k, 0, ∞}];

pAB[ν_Integer, α_, xaA_, xbA_, xaB_, xbB_] :=
Exp[-ν (σA[α, xaA, xbA] + σB[α,
xaB, xbB] )] Sum[(ν σA[α, xaA, xbA])^
k  (ν σB[α, xaB, xbB])^k/(k! k!), {k,
0, ∞}];

pBA[ν_Integer, β_, xaA_, xbA_, xaB_, xbB_] :=
Exp[-ν (τA[β, xaA, xbA] + τB[β, xaB,
xbB] )] Sum[(ν τA[β, xaA, xbA])^
k (ν τB[β, xaB, xbB])^(k + 1)/(k! (k + 1)!), {k,
0, ∞}];

pBB[ν_Integer, β_, xaA_, xbA_, xaB_, xbB_] :=
Exp[-ν (τA[β, xaA, xbA] + τB[β, xaB,
xbB] )] Sum[(ν τA[β, xaA, xbA])^
k (ν τB[β, xaB, xbB])^k/(k! k!), {k,
0, ∞}];

w = NSolve[
{
-1 < xbA < 1 &&
-1 < xaA < 1 &&
-1 < xbB < 1 &&
-1 < xaB < 1 &&

vA[xaA]/vB[xaA] α/(1 - β) ==
vA[xbA]/vB[xbA] (1 - α)/β  &&
- vA[xbA]/vB[xbA] ==
pBA[ν, β, xaA, xbA, xaB, xbB]/
pAA[ν, α, xaA, xbA, xaB,
xbB] * β/(1 - α) &&

vA[xaB]/vB[xaB] α/(1 - β) ==
vA[xbB]/vB[xbB] (1 - α)/β  &&
- vA[xbB]/vB[xbB] ==
pBB[ν, β, xaA, xbA, xaB, xbB]/
pAB[ν, α, xaA, xbA, xaB,
xbB] * β/(1 - α)
},
{xaA, xbA, xaB, xbB}, Reals]


The Bessel functions are the infinite sums in the pAA etc. functions. When trying to solve this mathematica's output is just the input to NSolve with the equations inserted to it. I was able to solve a similar system of two instead on four equations (which I can post here if helpful) and I don't understand why this system should be that difficult to solve.

Any help is highly appreciated!!

PS: I have tried to use = instead of := for the first functions because I thought is was an issue with re-evaluation but that this not help either. I have also given start values to the x which didn't help either. I hope the code is in a copyable fashion. I tried and it worked for me but if there are any changes I should make before you can help me, please let me now.)

• 1. A wide range of Bessel functions and related are implemented natively. Using those may simplify your code. 2) Solve and NSolve work best with polynomial equations; try FindRoot in your case instead. – MarcoB May 28 '16 at 13:58

In addition to using FindRoot as pointed out by Mariusz Iwaniuk, use Set rather than SetDelayed in your definitions for efficiency (e.g., each infinite sum is calculated just once). Use exact constants to facilitate simplifications (e.g., 7/10 vice .7)

ν = 10; c = 1/2; α = 7/10; β = 55/100;
vA[x_] = -c + x;
vB[x_] = c + x;
F[x_] = 1/2 + x/2;


Note that in your equation for σA you left out a space in α F[xaA]

σA[α_, xaA_, xbA_] = α F[xaA] + (1 - α) F[xbA];

σB[α_, xaB_,
xbB_] = α (1 - F[xaB]) + (1 - α) (1 - F[xbB]);
τA[β_, xaA_, xbA_] = (1 - β) F[xaA] + β F[xbA];
τB[β_, xaB_,
xbB_] = (1 - β) (1 - F[xaB]) + β (1 - F[xbB]);

pAA[ν_Integer, α_, xaA_, xbA_, xaB_, xbB_] =
Exp[-ν (σA[α, xaA, xbA] + σB[α, xaB,
xbB])] Sum[(ν σA[α, xaA, xbA])^
k (ν σB[α, xaB, xbB])^(k + 1)/(k! (k + 1)!), {k,
0, ∞}] //
Simplify[#, {-1 < xbA < 1, -1 < xaA < 1, -1 < xbB < 1, -1 < xaB < 1}] &;

pAB[ν_Integer, α_, xaA_, xbA_, xaB_, xbB_] =
Exp[-ν (σA[α, xaA, xbA] + σB[α, xaB,
xbB])] Sum[(ν σA[α, xaA, xbA])^
k (ν σB[α, xaB, xbB])^k/(k! k!), {k,
0, ∞}] //
Simplify[#, {-1 < xbA < 1, -1 < xaA < 1, -1 < xbB < 1, -1 < xaB < 1}] &;

pBA[ν_Integer, β_, xaA_, xbA_, xaB_, xbB_] =
Exp[-ν (τA[β, xaA, xbA] + τB[β, xaB,
xbB])] Sum[(ν τA[β, xaA, xbA])^
k (ν τB[β, xaB, xbB])^(k + 1)/(k! (k + 1)!), {k,
0, ∞}] //
Simplify[#, {-1 < xbA < 1, -1 < xaA < 1, -1 < xbB < 1, -1 < xaB < 1}] &;

pBB[ν_Integer, β_, xaA_, xbA_, xaB_, xbB_] =
Exp[-ν (τA[β, xaA, xbA] + τB[β, xaB,
xbB])] Sum[(ν τA[β, xaA, xbA])^
k (ν τB[β, xaB, xbB])^k/(k! k!), {k, 0, ∞}] //
Simplify[#, {-1 < xbA < 1, -1 < xaA < 1, -1 < xbB < 1, -1 < xaB < 1}] &;

eqns = {vA[xaA]/vB[xaA] α/(1 - β) ==
vA[xbA]/vB[xbA] (1 - α)/β && -vA[xbA]/vB[xbA] ==
pBA[ν, β, xaA, xbA, xaB, xbB]/
pAA[ν, α, xaA, xbA, xaB, xbB]*β/(1 - α) &&
vA[xaB]/vB[xaB] α/(1 - β) ==
vA[xbB]/vB[xbB] (1 - α)/β && -vA[xbB]/vB[xbB] ==
pBB[ν, β, xaA, xbA, xaB, xbB]/
pAB[ν, α, xaA, xbA, xaB, xbB]*β/(1 - α)} //
Simplify[#, {-1 < xbA < 1, -1 < xaA < 1, -1 < xbB < 1, -1 < xaB < 1}] &;

FindRoot[eqns, {{xaA, .1}, {xbA, .1}, {xaB, .1}, {xbB, .1}}]

(*  {xaA -> 0.0942508 + 0. I, xbA -> -0.160697 + 0. I, xaB -> 0.109664 + 0. I,
xbB -> -0.14613 + 0. I}  *)


As in this example, you may need to use Chop to remove imaginary artifacts due to machine precision.

W = % // Chop

(*  {xaA -> 0.0942508, xbA -> -0.160697, xaB -> 0.109664, xbB -> -0.14613}  *)

• Yours answer is better than mine :) – Mariusz Iwaniuk May 28 '16 at 13:54
• Many thanks, It works and is way faster now! – JKJ May 29 '16 at 16:17
• Regarding the "Set" instead of "Set Delayed", you did this for the sigma functions as well which then go into the infinite sums as arguments; don't they have to be re-evaluated when using FindRoot? – JKJ May 29 '16 at 16:22
• @JKJ - look at the stored definition for eqns. The equations no longer have any explicit dependence on σA or σB. The eqns contain only functions/expressions containing {xaA, xbA, xaB, xbB}. FindRoot will manipulate these four variables to locate the approximate root. – Bob Hanlon May 30 '16 at 0:31

With FindRoot works.

 w = FindRoot[{vA[xaA]/vB[xaA] α/(1 - β) ==
vA[xbA]/vB[xbA] (1 - α)/β, -vA[xbA]/vB[xbA] ==
pBA[ν, β, xaA, xbA, xaB, xbB]/
pAA[ν, α, xaA, xbA, xaB, xbB]*β/(1 - α),
vA[xaB]/vB[xaB] α/(1 - β) ==
vA[xbB]/vB[xbB] (1 - α)/β, -vA[xbB]/vB[xbB] ==
pBB[ν, β, xaA, xbA, xaB, xbB]/
pAB[ν, α, xaA, xbA, xaB,
xbB]*β/(1 - α)}, {{xaA, 1/2}, {xbA, 1/2}, {xaB, 1/4},{xbB, 1/4}}]


{xaA -> 0.0942508, xbA -> -0.160697, xaB -> 0.109664, xbB -> -0.14613}

If you change starting search point:

  w = FindRoot[{vA[xaA]/vB[xaA] α/(1 - β) ==
vA[xbA]/vB[xbA] (1 - α)/β, -vA[xbA]/vB[xbA] ==
pBA[ν, β, xaA, xbA, xaB, xbB]/
pAA[ν, α, xaA, xbA, xaB,
xbB]*β/(1 - α),
vA[xaB]/vB[xaB] α/(1 - β) ==
vA[xbB]/vB[xbB] (1 - α)/β, -vA[xbB]/vB[xbB] ==
pBB[ν, β, xaA, xbA, xaB, xbB]/
pAB[ν, α, xaA, xbA, xaB,
xbB]*β/(1 - α)}, {{xaA, #1}, {xbA, #2}, {xaB, #3}, {xbB, #4}}] &;

w[1/1, 1/2, 1/2, 1/2]


{xaA -> 0.401414, xbA -> 0.262252, xaB -> 0.709528, xbB -> 1.4764}

xbB does not meet the inequalities -1 < xbB < 1.

• I'm a little embarrassed by the time I have tried unsuccessfully. (At some point i tried FindRoot but must have done it wrong. Thanks a thousand times!! – JKJ May 28 '16 at 13:10
• @JKJ.NSOLVE deals primarily with linear and polynomial equations,therefore does not work with yours equation. – Mariusz Iwaniuk May 28 '16 at 13:27