Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to solve a equation with Newton's method via FindRoot, and the codes are:

Define the functions:

τs[Ωs_,t_]:=2*Ωs*Coth[Ωs/(2*t)]
τb[Ωb_,t_]:=2*Ωb*Coth[Ωb/(2*t)]
Xs[Js_,Jb_,Ωs_,t_,δ_]:=(τs[Ωs,t]-4*Js-4*Js/(Sqrt[2])^δ)/(Jb+4*Js/(Sqrt[2])^δ+4*Js/(Sqrt[3])^δ)
Xb[Jb_,Ωb_,t_,δ_]:=(τb[Ωb,t]-4*Jb-4*Jb/(Sqrt[2])^δ)/(Jb+4*Jb/(Sqrt[2])^δ+4*Jb/(Sqrt[3])^δ)
ϕ[Jb_,Ωb_,t_,δ_]:=ArcCosh[Xb[Jb,Ωb,t,δ]/2]
BN[Jb_,Ωb_,t_,δ_,n_]:=Sinh[(n+1)*ϕ[Jb,Ωb,t,δ]]/Sinh[ϕ[Jb,Ωb,t,δ]]
fumfa[Js_,Jb_,Ωs_,Ωb_,t_,δ_,n_]:=(Xs[Js,Jb,Ωs,t,δ])^2*BN[Jb,Ωb,t,δ,n-2]-2*Xs[Js,Jb,Ωs,t,δ]*BN[Jb,Ωb,t,δ,n-3]+BN[Jb,Ωb,t,δ,n-4]

Solve a series of solutions:

Table[FindRoot[fumfa[2.0,1.0,Ωs,2.0,3.5,t,4]==0.0,{Ωs,4.0}],{t,0.1,4.0,0.2}]

However, I got the following warning message:

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.

And some of the solutions have imaginary parts though the coefficients are actually zeros:

{{Ωs->11.6453 +0. I},{Ωs->10.5823 +0. I},{Ωs->9.58167 +0. I},{Ωs->8.63181 +0. I},{Ωs->7.72096 +0. I},{Ωs->6.83656 +0. I},{Ωs->5.96386 +0. I},{Ωs->5.08324 +0. I},{Ωs->4.16381 +0. I},{Ωs->3.14345 +0. I},{Ωs->1.82232},{Ωs->5.70687*10^-6},{Ωs->4.74748},{Ωs->3.29661},{Ωs->1.28711},{Ωs->-1.26412*10^-6},{Ωs->-3.05724*10^-6},{Ωs->-0.0000392875},{Ωs->0.0000453593},{Ωs->0.0000543535}}

My questions are:

  1. How to remove the warning message?
  2. why there are imaginary parts in the solutions and how to assure Reals only computation?
share|improve this question
    
Do you want a "Reals only computation" or is removing the imaginary numbers afterwards acceptable? –  Öskå Aug 1 at 12:03
    
Re[]/.solutions can remove imaginary parts easily; I am wondering how to compute via reals only method. –  LCFactorization Aug 1 at 12:05
2  
You can also use Chop to remove very small imaginaries caused by numerical methods applied to inexact numbers, which is your situation. –  m_goldberg Aug 1 at 12:11
1  
BTW for the first question, just use Quiet@FindRoot[]. –  Öskå Aug 1 at 12:15
1  
Often the warning message indicates that FindRoot got stuck at a non-root. You should check your solutions to see if they are accurate enough. –  Michael E2 Aug 1 at 14:57

1 Answer 1

up vote 2 down vote accepted

First, you have t in the slot for δ -- that may be a mistake. Be that as it may, the the question about the warning FindRoot::lstol has an explanation.

Second, you're getting complex solutions because the function evaluates to complex numbers:

fumfa[2.0, 1.0, Ωs, 2.0, 3.5, t, 4] /. {t -> 0.1, Ωs -> 4.0}
(* 1.02171 + 0. I *)

Finally, the explanation. The warning FindRoot::lstol means that Mathematica was having trouble making the function(s) close to zero. This can happen at a local extremum, for instance.

We can use Check to catch the values of t that produce warnings. The value of t and the corresponding "solution" sol0 will be stored in err. All solutions together with their corresponding values of t will be stored in sol.

Module[{sol0},
  err = {};
  sols = Table[
    Quiet@Check[
      sol0 = Flatten @
        {HoldPattern[t] -> t, FindRoot[fumfa[2.0, 1.0, Ωs, 2.0, 3.5, t, 4] == 0.0, {Ωs, 4.0}]},
      AppendTo[err, sol0]; sol0,
      {FindRoot::lstol}],
    {t, 0.1, 4.0, 0.2}]
  ];

Here are the values of t that produce warnings and the corresponding values of the function:

t /. err
fumfa[2.0, 1.0, Ωs, 2.0, 3.5, t, 4] /. err
(*
  {2.3, 3.1, 3.3, 3.5, 3.7, 3.9}
  {-0.0643143, 0.298996, 0.960447, 2.04808, 3.70264, 6.09244}
*)

You can have such large errors when the derivative of the objective function is very large, since $dy = f'(x)\;dx$. That is not the case here:

D[fumfa[2.0, 1.0, Ωs, 2.0, 3.5, t, 4], Ωs] /. err
(* {-2.75827*10^-7, -1.67904*10^-7, -7.33596*10^-7, -0.0000152322, 0.0000266041, 0.0000460977} *)

In fact, we see that FindMinimum has gotten stuck at relative extrema:

GraphicsGrid @ Partition[
  Hold[Plot[
      fumfa[2.0, 1.0, x, 2.0, 3.5, t, 4], {x, Ωs - 0.1, Ωs + 0.1}, 
      PlotLabel -> HoldForm["t" == t], 
      AxesLabel -> {"Ωs", "fumfa"}]] /. err // 
   ReleaseHold,
  3]

Mathematica graphics

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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