I think a bug in Symbolics package, but not sure

General Usage
1

I try the code

using Symbolics
D(f,x)=expand_derivatives(Differential(x)(f))
@variables H,S,r1
K1=10000.0; K2=6000.0;
a1=0.5; a2=0.5; r2=0.25;
dHdt(H,S)=r1*H*(1-(a1*S+H)/K1)
dSdt(H,S)=r2*S*(1-(S+a2*H)/K2)
f(H,S)=dHdt(H,S)+dSdt(H,S)

dfdH=expand(D(f(H,S),H))
dfdS=expand(D(f(H,S),S))
Asym=-[D(dfdH,H) D(dfdH,S);D(dfdS,H) D(dfdS,S)]
bsym=[substitute(dfdH,[H=>0,S=>0]),
      substitute(dfdS,[H=>0,S=>0])]
Xsym=Asym\bsym
substitute(Xsym[1],r1=>0.10)
simplify(Xsym[1])

and the last two evaluations return

julia> substitute(Xsym[1],r1=>0.10)
1876.3029881862406

julia> simplify(Xsym[1])
0

This looks like a bug to me because the substitute clearly indicates Xsym[1] is nonzero while simplify simplifies it to zero. The answer 1876.3… is correct. Is it possible I’m just using the simplify function incorrectly. What is simplify supposed to do?

2 Likes
2

I’m pretty sure something has gone wrong with simplify. I simplified by hand to the following expression assigned y, and not only is simplification to 0 more obviously unlikely, simplify affects it differently.

julia> y = ((8.333333333333333e-5 - 0.25 * 5e-5)*r1 - 0.25 * 2.0833e-5)/(8.333333333333333e-5 * 0.0002 * r1 -(2.0833333333333333e-5 + 5e-5 * r1)^2)
(-5.20825e-6 + 7.083333333333332e-5r1) / (1.6666666666666667e-8r1 - ((2.0833333333333333e-5 + 5.0e-5r1)^2))

julia> simplify(y)
1 / (0.00020328723875432528 - 3.529411764705883e-5r1)

To justify y being “equal” to Xsym[1], I substituted some values for r1 for comparison. With the exception of the Num(0)/Num(0) === Num(1) issue of Xsym[1], the values of y are close enough that the difference is probably caused by floating point precision and my manual entry of floating point literals for y. simplify(y) however is really off.

julia> substitute.([Xsym[1] y simplify(y)], [r1=>x for x in 0.0:0.1:1.0])
11×3 Matrix{Num}:
       1  11999.8  4919.15
  1876.3  1876.39  5006.06
 3759.84  3759.87   5096.1
 4316.95  4316.97  5189.44
 4625.64  4625.66  5286.26
  4846.8  4846.81  5386.76
 5028.56  5028.57  5491.16
 5190.48  5190.49  5599.69
 5342.08  5342.09  5712.59
 5488.64  5488.65  5830.14
 5633.38  5633.39  5952.62