# 3D plot using Table & FindRoot

I have two expressions involving terms $$S_1$$ and $$S_2$$, call the expressions $$F1, F2$$.

I cannot solve $$F_i=1$$ for $$S_i$$ so instead I numerically approximate using FindRoot and giving specific parameter values for $$\alpha_{ij},\beta_i,\delta_i,\gamma$$.

I am interested in making a 3D plot, ListPlot3D where the x-axis is $$\delta_2$$ the y-axis is $$\beta_2$$ and the z-axis will be my $$S_i$$ values, but I haven't gotten this far yet.

For this I believe I need a list of $$S_i$$ values where $$\beta_2,\delta_2$$ are varying, so I attempt this with Table but I am not getting any results. It just runs forever. Am I doing something wrong? Is there a way to speed this up?

Here is the code I am using edited for stack exchange:

F1[S1_, S2_, γ_, β1_, β2_, δ1_, δ2_,α11_, α12_] := (γ*β1)/δ1*Integrate[Exp[-γ*a]*(1 - Exp[-δ1*a])*Exp[-S1*(α11*β1)/δ1 (1 - Exp[-δ1*a]) - S2*(α12*β2)/δ2 (1 - Exp[-δ2*a])], {a,0,[Infinity]}]

F2[S1_, S2_, γ_, β1_, β2_, δ1_, δ2_, α22_, α21_] := (γ*β2)/δ2*Integrate[Exp[-γ*a]*(1 - Exp[-δ2*a])*Exp[-S2*(α22*β2)/δ2 (1 - Exp[-δ2*a]) - S1*(α21*β1)/δ1 (1 - Exp[-δ1*a])], {a,0, \[Infinity]}]

V= With[{β1 = 1.1, γ = .65, δ1 = 1, α11 = 1, α12 = 1, α22 = 1, α21 = 1}, Table[FindRoot[{F1[S1, S2, γ, β1, β2, δ1, δ2, α11, α12] - 1, F2[S1, S2, γ, β1, β2, δ1, δ2, α22, α21] - 1}, {{S1, 0.1}, {S2, 0.1}}], {δ2, 0.1, 1.5, 0.1}, {β2, 0.1, 1.5, 0.1}]]


Here is my actual code:

F1[S1_, S2_, \[Gamma]_, \[Beta]1_, \[Beta]2_, \[Delta]1_, \[Delta]2_, \[Alpha]11_, \[Alpha]12_] := (\[Gamma]*\[Beta]1)/\[Delta]1*Integrate[Exp[-\[Gamma]*a]*(1 - Exp[-\[Delta]1*a])*Exp[-S1*(\[Alpha]11*\[Beta]1)/\[Delta]1 (1 - Exp[-\[Delta]1*a]) - S2*(\[Alpha]12*\[Beta]2)/\[Delta]2 (1 - Exp[-\[Delta]2*a])], {a,0, \[Infinity]}]

F2[S1_, S2_, \[Gamma]_, \[Beta]1_, \[Beta]2_, \[Delta]1_, \[Delta]2_, \[Alpha]22_, \[Alpha]21_] := (\[Gamma]*\[Beta]2)/\[Delta]2*Integrate[Exp[-\[Gamma]*a]*(1 - Exp[-\[Delta]2*a])*Exp[-S2*(\[Alpha]22*\[Beta]2)/\[Delta]2 (1 - Exp[-\[Delta]2*a]) - S1*(\[Alpha]21*\[Beta]1)/\[Delta]1 (1 - Exp[-\[Delta]1*a])], {a,0, \[Infinity]}]

V=With[{\[Beta]1 = 1.1, \[Gamma] = .65, \[Delta]1 = 1, \[Alpha]11 = 1, \[Alpha]12 = 1, \[Alpha]22 = 1, \[Alpha]21 = 1}, Table[FindRoot[{F1[S1, S2, \[Gamma], \[Beta]1, \[Beta]2, \[Delta]1, \[Delta]2, \[Alpha]11, \[Alpha]12] - 1, F2[S1, S2, \[Gamma], \[Beta]1, \[Beta]2, \[Delta]1, \[Delta]2, \[Alpha]22, \[Alpha]21] - 1}, {{S1, 0.1}, {S2, 0.1}}], {\[Delta]2, 0.1, 1.5, 0.1}, {\[Beta]2, 0.1, 1.5, 0.1}]]


Doing this for specific values using With gives me back $$S_i$$ vals but the Table does not.

With[{γ = 0.5,β1 = 2, β2 = 2,δ1 = 1, δ2 = 1,α11 = 1, α12 = 1, α22 = 1, α21 = 1},
FindRoot[{F1[S1, S2, γ, β1, β2, δ1, δ2,α11, α12] - 1, F2[S1, S2, γ, β1, β2, δ1, δ2, α22, α21] - 1}, {{S1, 0.1}, {S2, 0.1}}]]


gives me back: {S1 -> 0.0908753, S2 -> 0.0908753}

Update: I am now trying:

Table[FindRoot[{F1[S1, S2, .65, 1.1, β2, 1, δ2, 1, 1] - 1, F2[S1, S2, .65, 1.1, β2, 1, δ2, 1, 1] - 1}, {{S1, 0.1}, {S2, 0.1}}],{δ2, 0.1, 1.5, 0.1}, {β2, 0.1, 1.5, 0.1}]


10 mins and counting...

Try inspecting the plots of {s1,s2,f1[...]-1} and {s1,s2,f2[...]-1} where ... is your parameter set. (I've replaced your functions with f1 and f2 using NIntegrate to speed things up):

f1[S1_, S2_, γ_, β1_, β2_, δ1_, δ2_,α11_, α12_] := (γ*β1)/δ1*NIntegrate[Exp[-γ*a]*(1 - Exp[-δ1*a])*Exp[-S1*(α11*β1)/δ1 (1 - Exp[-δ1*a]) - S2*(α12*β2)/δ2 (1 - Exp[-δ2*a])], {a,0,[Infinity]}]

f2[S1_, S2_, γ_, β1_, β2_, δ1_, δ2_, α22_, α21_] := (γ*β2)/δ2*NIntegrate[Exp[-γ*a]*(1 - Exp[-δ2*a])*Exp[-S2*(α22*β2)/δ2 (1 - Exp[-δ2*a]) - S1*(α21*β1)/δ1 (1 - Exp[-δ1*a])], {a,0, ∞}]

t = Flatten[#, 1] &@With[
{β1 = 1.1, γ = .65, δ1 = 1, α11 = 1, α12 = 1, α22 = 1, α21 = 1, β2 = 0.1, δ2 = 0.1},
Table[
{s1, s2,
f1[s1,
s2, γ, β1, β2, δ1, δ2, α11, α12] - 1,
f2[s1,
s2, γ, β1, β2, δ1, δ2, α22, α21] - 1
}
, {s1, 0, 1, 0.1}, {s2, -1, 1, 0.1}]
];


Plotting these values you see that for the starting parameters in your Table there is no intersection of the two planes. I added a plane at 0 to make it obvious that there is likely no simultaneous solution of f1-1==0 and f2-1==0, even if they both might reach 0 at some point {s1,s2}.

ListPlot3D[{
t[[;; , {1, 2, 3}]],
t[[;; , {1, 2, 4}]], {1, 1, 0}*# & /@ t[[;; , {1, 2, 3}]]
}]


In the special case you gave outside of Table where you were able to find a solution all the parameters are equal, so the two planes lie on top of one another and any point on the line f1-1==0 works.

t = Flatten[#, 1] &@With[
{γ = 0.5, β1 = 2, β2 = 2, δ1 = 1, δ2 = 1, α11 = 1, α12 = 1, α22 = 1, α21 = 1},
Table[
{s1, s2,
f1[s1,
s2, γ, β1, β2, δ1, δ2, α11, α12] - 1,
f2[s1,
s2, γ, β1, β2, δ1, δ2, α22, α21] - 1
}
, {s1, 0, 1, 0.1}, {s2, 0, 1, 1}]
];

ListPlot3D[{
t[[;; , {1, 2, 3}]],
t[[;; , {1, 2, 4}]], {1, 1, 0}*# & /@ t[[;; , {1, 2, 3}]]
}]


It's not obvious, but the orange and blue planes are the same.

• Is β2 and δ2 varying here? it seems you have them fixed at β2=2 and δ2=1. I want to look at S_i values when we vary β2 and δ2 but keep the other parameters fixed. Apr 6 at 15:07
• Yes, the first plot has one set of parameters, and the second a different set (although some other parameters change as well). The point here is that for some values of your parameters there is no {s1,s2} where the solutions are both 0, so I wouldn't expect FindRoot to find an answer. This is why your table is taking so long - in the first case I used the first point the table would test, which has no solns. In the case you tested outside of the table there are an infinite number of solutions. Syntactically your code looks fine, I think the problem is mathematical. Apr 6 at 17:18