1
$\begingroup$

I am trying to solve numerically an equation and generate some results. I use the following code

u[c_] := (c^(1 - σ) - 1)/(1 - σ)
f[s_] := g s (1 - s/sbar1)
h[s_] := (2 hbar)/(1 + Exp[η (s/sbar - 1)])
co[a_] := ϕ (a^2)/2
ψ[k_] := wbar (ω + (1 - ω) Exp[-γ k])

The equation I try to solve is the following

adap[k_, s_] := (ρ + δ) u'[f[s] - priceadap δ k] co'[δ k] + ψ'[k] h[s]

I have the following constant parameter set

paramFinal2 = {σ -> 1.7, ρ -> 0.025, g -> 0.05, sbar -> 10, η -> 11, hbar -> 0.5,  priceadap -> 0.0006, γ -> 0.6, χ -> 1000, ϕ -> 0.05, ω -> 0.35, β -> 0.8, δ -> 0.065, sbar1 -> 10, wbar -> 1000};

So, for different values of $s$, I try to generate the corresponding values of $k$.

For this, I use the following code

tmax1 = 10;
solK[i_] := Solve[adap[k, i] == 0 /. paramFinal2, k];
Table[solK[i], {i, 1, tmax1}];

Unfortunately, this does not give any result. Mathematica is always on mode "Running...".

P.S I am using Mathematica 9.0

$\endgroup$
1
$\begingroup$

I always find it a good idea to plot the functions you are investigating. Thus

tmax1 = 10;
Column@Table[
  Plot[Evaluate[adap[k, i] /. paramFinal2], {k, 0, 20}], {i, 1 tmax1}]

Mathematica graphics

They all look fine except for the last one. The zero is at about k = 10 for most of the functions. The equation is not a polynomial and thus may not be solvable using Solve. I don't know what range of equations Solve can deal with but some functions are just too complicated. The solution for this is to use FindRoot.

rts = Table[
  FindRoot[adap[k, i] /. paramFinal2, {k, 10}], {i, 1, tmax1}]

which gives

{{k -> 10.7363}, {k -> 12.1674}, {k -> 12.8497}, {k -> 13.185}, {k -> 
   13.2841}, {k -> 13.169}, {k -> 12.797}, {k -> 12.0139}, {k -> 
   10.3238}, {k -> 3.43075*10^39 - 3.28305*10^36 I}}

The last root is clearly wrong and is complex which gives us a clue. If we do

adap[k, 10] /. paramFinal2 /. k -> 1

We get

5269.92 + 7400.72 I

So this looks like a complex function. Thus

Plot[Evaluate[ReIm[adap[k, 10] /. paramFinal2]], {k, 0, 20}]

Mathematica graphics

This suggests that you might not have a real root for this case.

Hope that helps.

$\endgroup$
1
$\begingroup$

Neither Solve nor NSolve will handle your equation. As Hugh has shown, you can use FindRoot. You could also use FindInstance as follows.

solK[i_] := FindInstance[adap[k, i] == 0 /. paramFinal2, k, Reals]
With[{tmax1 = 10}, Flatten @ Table[solK[i], {i, tmax1}]][[All, 2]]
{10.7363, 12.1674, 12.8497, 13.185, 13.2841, 13.169, 12.797, 12.0139, 10.3238}
$\endgroup$
0
$\begingroup$

one approach here is to use ContourPlot

p = ContourPlot[(adap[k, s] /. paramFinal2) == 0, {s, 0, 10}, 
  {k, 0, 30}, PlotPoints -> 100]

enter image description here

interpolate the curve and extract specific values:

g = Interpolation[First@Cases[Normal@p, Line[x_] :> x, Infinity]];
Table[{s, g[s]}, {s, 9}]

{{1, 10.7315}, {2, 12.1675}, {3, 12.8475}, {4, 13.1713}, {5, 13.2854}, {6, 13.1591}, {7, 12.7972}, {8, 12.0124}, {9, 10.3233}}

note ContourPlot only finds points to adequate precison for plotting purposes. If you need better precision use as initial guesses for FindRoot:

Table[{s, k /. FindRoot[adap[k, s] /. paramFinal2, {k, g[s]}]},
 {s,9}]

{{1, 10.7363}, {2, 12.1674}, {3, 12.8497}, {4, 13.185}, {5, 13.2841}, {6, 13.169}, {7, 12.797}, {8, 12.0139}, {9, 10.3238}}

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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