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?
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
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.
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.
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.
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. Redirecting
It should be relatively simple to extend my code above to generate such reconstruction expressions for general stencils.
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)
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?
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?
Would it be possible to integrate the resulting polynomial using Polynomials?
(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 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.
As for collect, SymPy’s