On Sunday, November 6, 2016 at 3:00:30 PM UTC, Francis Banks wrote:
>
> I am solving a polynomial which arises from plotting titration cures in
> chemistry. The rule of signs suggests it has one positive root. Find_root
> seems to find it. Solve with poly_solve=true does not. Instead it gives 4
> complex roots, which don't appear to satisfy the equation. They do not even
> form two conjugate pairs. Here is some code which illustrates the problem:
>
> var('Kb','b0','K1','Ks','a0')
> fx(x)=-Kb*x^4 - (b0*Kb + K1*Kb + Ks)*x^3 - (b0*K1*Kb - a0*K1*Kb + K1*Ks -
> Kb*Ks)*x^2 + (a0*K1*Ks + K1*Kb*Ks + Ks^2)*x + K1*Ks^2
>
> fx(x)=fx.subs(K1=10^-2.15,Ks=10^-14,Kb=10^5.0,a0=0.1,b0=0)
> s=find_root(fx,0,0.025,1E-14)
> show(s)
> s1=solve(fx,x,to_poly_solve=true)
> show('roots from solve:
> ',N(s1[0].rhs(),3));show(N(s1[1].rhs(),3));show(N(s1[2].rhs(),3));show(N(s1[3].rhs(),3))
> show('Substituting first root from solve: ',N(fx(s1[0].rhs()),10))
> show('Substituting first root from find_root: ',N(fx(s),10))
> plot(fx,0.001,0.03)
>
>
> Can someone shed light on this?
>
You are working very close to the numerical precision allowed by the usual
float doubles (in particular with Ks^2 being 10^(-28)).
With degree 4 polynomials you might be better off using the exact formulae
for their roots.
--
You received this message because you are subscribed to the Google Groups
"sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to sage-support+unsubscr...@googlegroups.com.
To post to this group, send email to sage-support@googlegroups.com.
Visit this group at https://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.
