Symbolics and substitution?

I’m trying to use Symbolics.jl to fit a polynomial to some data points (\ldots,(x_j,y_j),\ldots) by Lagrange interpolation, then take the derivative of the polynomial, and substitute some relationships between the data points.

I get an error message when I do the substitution. To make some sense of the syntax, I use D to denote \Delta x, xim3d2 to denote x_{i-\frac{3}{2}}, etc. [Suggestions for better notation would be appreciated…]

So here is what happens:

using Symbolics
#
@variables x D
@variables xim5d2 xim3d2 xim1d2 xip1d2
@variables uim2 uim1 ui
#
xd = [xim5d2,xim3d2,xim1d2,xip1d2]
yd = [0,uim2*D, (uim2+uim1)*D, (uim2+uim1+ui)*D]
#
function lagrange_polynomial(xd,yd)
    n = length(xd)
    LL = 0
    for j in 1:n
        L = 1
        for k in setdiff(1:n,j)
            L = L*(x-xd[k])/(xd[j]-xd[k])
        end
        LL = LL + yd[j]*L
    end
    return LL            
end
#
L = lagrange_polynomial(xd,yd)
u = derivative(L,x)
#
xsubs = Dict([xim3d2=>xim5d2+D, xim1d2=>xim3d2+D, xip1d2=>xim1d2+D])
#
substitute(u,xsubs)

The last line gives an error message:

TypeError: non-boolean (Num) used in boolean context.....

Question: What is it that I do wrong in the substitute command?

It seems to work in Julia Version 1.6.0-beta1.0 on Win 10.
PS: had to specify: u = Symbolics.derivative(L,x)

julia> substitute(u,xsubs)
uim2*(x - xim5d2)*(x - (D + xim1d2))*((D + xim5d2 - (D + xim1d2))^-1)*((D + xim5d2 - (D + xim3d2))^-1) + uim2*(x - xim5d2)*(x - (D + xim3d2))*((D + xim5d2 - (D + xim1d2))^-1)*((D + xim5d2 - (D + xim3d2))^-1) + uim2*(x - (D + xim1d2))*(x - (D + xim3d2))*((D + xim5d2 - (D + xim1d2))^-1)*((D + xim5d2 - (D + xim3d2))^-1) + D*(uim1 + uim2)*(x - xim5d2)*(x - (D + xim1d2))*((D + xim3d2 - xim5d2)^-1)*((D + xim3d2 - (D + xim1d2))^-1)*((D + xim3d2 - (D + xim5d2))^-1) + D*(uim1 + uim2)*(x - xim5d2)*(x - (D + xim5d2))*((D + xim3d2 - xim5d2)^-1)*((D + xim3d2 - (D + xim1d2))^-1)*((D + xim3d2 - (D + xim5d2))^-1) + D*(uim1 + uim2)*(x - (D + xim1d2))*(x - (D + xim5d2))*((D + xim3d2 - xim5d2)^-1)*((D + xim3d2 - (D + xim1d2))^-1)*((D + xim3d2 - (D + xim5d2))^-1) + D*(x - xim5d2)*(x - (D + xim3d2))*(ui + uim1 + uim2)*((D + xim1d2 - xim5d2)^-1)*((D + xim1d2 - (D + xim3d2))^-1)*((D + xim1d2 - (D + xim5d2))^-1) + D*(x - xim5d2)*(x - (D + xim5d2))*(ui + uim1 + uim2)*((D + xim1d2 - xim5d2)^-1)*((D + xim1d2 - (D + xim3d2))^-1)*((D + 
xim1d2 - (D + xim5d2))^-1) + D*(x - (D + xim3d2))*(x - (D + xim5d2))*(ui + uim1 + uim2)*((D + xim1d2 - xim5d2)^-1)*((D + xim1d2 - (D + xim3d2))^-1)*((D + xim1d2 - (D + xim5d2))^-1)      
1 Like

Thanks – good to know. I’m on Julia v1.5.0 on my job desktop… I’m kind of waiting for v1.6.0 official release.

Still, I’m curious as to whether there is a known reason why it doesn’t work on v1.5.x.

I don’t get an error on Julia 1.5.3 on my Windows 10 PC.
Perhaps it is some kind of name conflict?

julia> versioninfo()
Julia Version 1.5.3
Commit 788b2c77c1 (2020-11-09 13:37 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: AMD Ryzen 9 3900X 12-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-9.0.1 (ORCJIT, znver2)

julia> substitute(u,xsubs)
uim2*(x - xim5d2)*(x - (D + xim1d2))*((D + xim5d2 - (D + xim1d2))^-1)*((D + xim5d2 - (D + xim3d2))^-1) + uim2*(x - xim5d2)*(x - (D + xim3d2))*((D + xim5d2 - (D + xim1d2))^-1)*((D + xim5d2 - (D + xim3d2))^-1) + uim2*(x - (D + xim1d2))*(x - (D + xim3d2))*((D + xim5d2 - (D + xim1d2))^-1)*((D + xim5d2 - (D + xim3d2))^-1) + D*(uim1 + uim2)*(x - xim5d2)*(x - (D + xim1d2))*((D + xim3d2 - xim5d2)^-1)*((D + xim3d2 - (D + xim1d2))^-1)*((D + xim3d2 - (D + xim5d2))^-1) + D*(uim1 + uim2)*(x - xim5d2)*(x - (D + xim5d2))*((D + xim3d2 - xim5d2)^-1)*((D + xim3d2 - (D + xim1d2))^-1)*((D + xim3d2 - (D + xim5d2))^-1) + D*(uim1 + uim2)*(x - (D + xim1d2))*(x - (D + xim5d2))*((D + xim3d2 - xim5d2)^-1)*((D + xim3d2 - (D + xim1d2))^-1)*((D + xim3d2 - (D + xim5d2))^-1) + D*(x - xim5d2)*(x - (D + xim3d2))*(ui + uim1 + uim2)*((D + xim1d2 - xim5d2)^-1)*((D + xim1d2 - (D + xim3d2))^-1)*((D + xim1d2 - (D + xim5d2))^-1) + D*(x - xim5d2)*(x - (D + xim5d2))*(ui + uim1 + uim2)*((D + xim1d2 - xim5d2)^-1)*((D + xim1d2 - (D + xim3d2))^-1)*((D + xim1d2 - (D + xim5d2))^-1) + D*(x - (D + xim3d2))*(x - (D + xim5d2))*(ui + uim1 + uim2)*((D + xim1d2 - xim5d2)^-1)*((D + xim1d2 - (D + xim3d2))^-1)*((D + xim1d2 - (D + xim5d2))^-1)

(@v1.5) pkg> st
Status `C:\Users\oheil\.julia\environments\v1.5\Project.toml`
...
  [0c5d862f] Symbolics v0.1.8
...

Hm. So it could be that Julia 1.5.0 downgrades the package, perhaps. The code I included above is the entire code I have in my Jupyter Notebook.

[I can try to run it in Visual Code to see if there is a problem with Jupyter Lab.]

I ran the code in the REPL, and I get the same error message, so the problem must be with my Julia set-up…perhaps.

The Polynomials package does Lagrange interpolation. If you’re just trying to fit and manipulate polynomials that might be a better alternative. Of course you might be exploring Symbolics, learning Lagrange interpolation, or something like that.

julia> using Polynomials

julia> x = [-1.0; 0.0; 1.0];

julia> y = [2.0; 3.0; 6.0]; # numerial values of x^2 + 2x + 3

julia> f = fit(x,y)
Polynomial(3.0 + 2.0*x + 1.0*x^2)

julia> f(2)
11.0

julia> dfdx = derivative(f)
Polynomial(2.0 + 2.0*x)

julia> dfdx(0)
2.0
4 Likes

I’m looking into symbolic expressions for reconstruction of interpolation values in PDE discretization.

Ah. I’ve used Polynomials for a similar purpose, interpolating numerical integrations of ODEs to compute intersections of trajectories with Poincare sections. For PDEs you’d probably need MultivariatePolynomials.jl. I haven’t used that but it might work for your purpose.

1 Like

I’m looking into symbolic schemes, so I need symbolic capabilities. I already did it by hand for a simple case, and I got the “text book” correct solution. I just want to automate it using Symbolics, so that I can find more schemes and for higher order without doing the “boring” manual computation.

In essence, I’m looking for interpolation schemes for the interface between nodes, and based on solutions at spatial nodes.

I’ve now tested the code on my laptop, too (Julia 1.5.0), with the same problem. So it seems that I need to update my version to Julia 1.5.3, or wait for the release of Julia 1.6.0.

1 Like

My set-up is on Symbolics v.0.1.2.

What happens if you do
] up

1 Like

I’m trying to do it right now… it takes some time on my 3 year old laptop…

1 Like

OK – Symbolics was updated to v 0.1.9 … I didn’t even think of this: I installed Symbolics 1-2 weeks ago. I guess I should have known that the developers work quickly… Anyway, now substitute works.

Next: my construct xsubs

xsubs = Dict([xim3d2=>xim5d2+D, xim1d2=>xim3d2+D, xip1d2=>xim1d2+D])

seems to be ill chosen. The substitute command doesn’t understand that it should do a recursive substitute. So I need to change this to:

xsubs = Dict([xim3d2=>xim5d2+D, xim1d2=>xim3d2+2D, xip1d2=>xim1d2+3D])

When I do this, I get the desired result:

julia> xsubs = Dict(xim3d2=>xim5d2+D, xim1d2=>xim5d2+2D, xip1d2=>xim5d2+3D);
julia> uu = simplify(substitute(u,xsubs);polynorm=true);
julia> simplify(substitute(uu,Dict(x=>xim5d2+3D));polynorm=true)
(11//6)∗𝑢𝑖+(1//3)∗𝑢𝑖𝑚2−((7//6)∗𝑢𝑖𝑚1)

which is consistent with

u_{i+\frac{1}{2}}^{(1)} = \frac {1}{3}\bar{u}_{i-2} - \frac{7}{6}\bar{u}_{i-1} + \frac{11}{6}\bar{u}_i

in Eq. 5, p. 106 of Y.-T. Zhang and C.-W. Shu: “ENO and WENO Schemes”, Chapter 5 in Handbook of Numerical Analysis, Vol. 5, Elsevier, 2016. http://dx.doi.org/10.1016/bs.hna.2016.09.009

It should be relatively simple to extend my code above to generate such reconstruction expressions for general stencils.

So I did a ] up now (after installing it a today) and got the 0.1.9 too, crazy speed :wink:

2 Likes

Well, Polynomials struggled a bit out of the box with this, but if that is the hammer you have, then well, you can get there through the Lagrange polynomial type, as follows:


using Polynomials, SpecialPolynomials

fit(Polynomial, xd, yd); # error

u = fit(Lagrange, xd, yd); # show error with `Num` coefficients...
uu = convert(Polynomial, u);
uu′ = Polynomials.derivative(uu);

xsubs = Dict(xim3d2=>xim5d2+D, xim1d2=>xim5d2+2D, xip1d2=>xim5d2+3D);
v = map(pⱼ -> simplify(substitute(pⱼ, xsubs), polynorm=true), uu′);
vv = v(xim5d2+3D)

simplify(vv, polynorm=true)

This results in the same answer as already shown:

(11//6)*ui + (1//3)*uim2 - ((7//6)*uim1)
3 Likes

A couple of questions:

  1. fit(Lagrange, xd, yd) doesn’t work for my set-up. I use Julia v. 1.5.0, with Symbolics v. 0.1.11, Polynomials v. 1.2.0, and SpecialPolynomials v. 0.1.10. What could be wrong?
  2. My own lagrange_polynomial() function works. Would it be possible to instead convert the output from that function to a polynomial? (convert(Polynomial,u)? If so, how do I specify the variable in the polynomial (x) in contrast to my symbolic parameters?
  3. Would it be possible to integrate the resulting polynomial using Polynomials?
1 Like

I have a slightly updated set of packages

(fit) pkg> st
Status `/private/tmp/fit/Project.toml`
  [f27b6e38] Polynomials v1.2.0
  [a25cea48] SpecialPolynomials v0.1.2
  [0c5d862f] Symbolics v0.1.11

And the fitting worked.

As for how to convert the expression to a polynomial, I had thought it would be easy. I tried something along the lines of this:

@variables a b c x
p = a*x^2 + b*x + c
z = variable(Polynomial)
q = substitute(p, Dict(x => z))

But that produced a Num expression, not a polynomial, as I had hoped. I would guess you would need someway to collect the expression (in SymPy’s use of that word) and then pick off each coefficient. This then could be passed along to Polynomials.

As to integrating, before this call vv = v(xim5d2+3D), v is a polynomial which can be integrated. But vv becomes a symbolic expression after the polynomial evaluation.

  • OK: I upgraded SpecialPolynomials to version v. 0.1.2, and now have the same version numbering as you report (except perhaps for the Julia v.). I still get an error message when I do fit(Lagrange,xd,yd).

  • I haven’t found a collect method in Symbolics. Since we are talking polynomials, I can always differentiate a number of times to find the coefficients of the polynomial, although this is somewhat clumsy.

There is only 0.1.0, 0.1.1, and 0.1.2 tagged. To test the above, I activated a new environment and just installed the three packages with add. (I too am waiting for 1.6 to sort out my mess of a default environment :slight_smile: As for collect, SymPy’s collect finds the common powers of a term in an expression. I’m not sure how that is implemented, perhaps through pattern matching, but differentiating might be one way to go about it if you know you have a polynomial with non-negative powers only.

Is the error in the fitting or the display? I have a semicolon to prevent calling show which needs updating to accommodate the Num type.