diff --git a/src/arrays.jl b/src/arrays.jl index 3c8b7b84e..5aecc990b 100644 --- a/src/arrays.jl +++ b/src/arrays.jl @@ -251,6 +251,7 @@ end # turn `f(x...)` into `term(f, x...)` # function call2term(expr, arrs=[]) + (expr isa QuoteNode) && return expr !(expr isa Expr) && return :($unwrap($expr)) if expr.head == :call if expr.args[1] == :(:) diff --git a/src/diff.jl b/src/diff.jl index 11c523dab..3a5eb07bd 100644 --- a/src/diff.jl +++ b/src/diff.jl @@ -33,17 +33,17 @@ struct Differential <: Operator x Differential(x) = new(value(x)) end -(D::Differential)(x) = Term{symtype(x)}(D, [x]) -(D::Differential)(x::Num) = Num(D(value(x))) -(D::Differential)(x::Complex{Num}) = wrap(ComplexTerm{Real}(D(unwrap(real(x))), D(unwrap(imag(x))))) +(D::Operator)(x) = Term{symtype(x)}(D, [x]) +(D::Operator)(x::Num) = Num(D(value(x))) +(D::Operator)(x::Complex{Num}) = wrap(ComplexTerm{Real}(D(unwrap(real(x))), D(unwrap(imag(x))))) SymbolicUtils.promote_symtype(::Differential, x) = x is_derivative(x) = istree(x) ? operation(x) isa Differential : false -Base.:*(D1, D2::Differential) = D1 ∘ D2 -Base.:*(D1::Differential, D2) = D1 ∘ D2 -Base.:*(D1::Differential, D2::Differential) = D1 ∘ D2 -Base.:^(D::Differential, n::Integer) = _repeat_apply(D, n) +Base.:*(D1, D2::Operator) = D1 ∘ D2 +Base.:*(D1::Operator, D2) = D1 ∘ D2 +Base.:*(D1::Operator, D2::Operator) = D1 ∘ D2 +Base.:^(D::Operator, n::Integer) = _repeat_apply(D, n) Base.show(io::IO, D::Differential) = print(io, "Differential(", D.x, ")") @@ -785,3 +785,95 @@ end function SymbolicUtils.substitute(op::Differential, dict; kwargs...) @set! op.x = substitute(op.x, dict; kwargs...) end + + +####################################################################################################################### +# Vector Calculus +####################################################################################################################### + +struct ArrayDifferentialOperator <: Operator + """The variables to differentiate with respect to.""" + vars + """The differentials, can be other functions if composite""" + differentials + """name""" + name + ArrayDifferentialOperator(vars, differentials, name) = new(vars, differentials, name) +end +Nabla(vars) = ArrayDifferentialOperator(value.(vars), map(Differential, scalarize(value.(vars))), "∇") +const Grad = Nabla +Div(vars) = (x) -> Nabla(vars) ⋅ x +Curl(vars) = (x) -> Nabla(vars) × x +Laplacian(vars) = Nabla(vars) ⋅ Nabla(vars) + +#? How to get transpose and Jac working? + +function (D::ArrayDifferentialOperator)(x::SymVec) + @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)." + @arrayop (i,) (D.differentials)[i](x[i]) term=D(x) +end +(D::ArrayDifferentialOperator)(x::Arr) = Arr(D(value(x))) + +function (D1::ArrayDifferentialOperator)(D2::ArrayDifferentialOperator) + @assert all(x -> any(isequal.((x,), D2.vars)), D1.vars) + + ArrayDifferentialOperator(D1.vars, scalarize(D1.differentials .∘ D2.differentials), "("*D1.name*"∘"*D2.name*")") +end + +function LinearAlgebra.dot(D::ArrayDifferentialOperator, x::SymVec) + @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)." + @show D(x), scalarize(D(x)) + sum(scalarize(D(x))) +end +LinearAlgebra.dot(D::ArrayDifferentialOperator, x::Arr) = Num(D ⋅ value(x)) + +function LinearAlgebra.dot(x::SymVec, D::ArrayDifferentialOperator) + @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)." + (y) -> sum(@arrayop (i,) x[i]*D.differentials[i](y) term = (x⋅D)(y)) +end +LinearAlgebra.dot(x::Arr, D::ArrayDifferentialOperator) = value(x) ⋅ D + +function LinearAlgebra.dot(D1::ArrayDifferentialOperator, D2::ArrayDifferentialOperator) + @assert all(scalarize(isequal.(D1.vars, D2.vars))) "Operators have different variables and cannot be composed." + lap = x -> sum((D1.differentials[i] ∘ D2.differentials[i])(x) for i in 1:length(D1.vars)) + (x) -> @arrayop (i,) lap(x[i]) term=(D1⋅D2)(x) reduce=+ +end + +function crosscompose(a, b) + v1 = x -> (a[2] ∘ b[3])(x) - (a[3] ∘ b[2])(x) + v2 = x -> (a[3] ∘ b[1])(x) - (a[1] ∘ b[3])(x) + v3 = x -> (a[1] ∘ b[2])(x) - (a[2] ∘ b[1])(x) + return [v1, v2, v3] +end + +function crosscall(a, b) + v1 = a[2](b[3]) - a[3](b[2]) + v2 = a[3](b[1]) - a[1](b[3]) + v3 = a[1](b[2]) - a[2](b[1]) + return [v1, v2, v3] +end +function LinearAlgebra.cross(D::ArrayDifferentialOperator, x::SymVec) + @assert length(D.vars) == length(x) == 3 "Cross product is only defined in 3 dimensions." + curl = crosscall(D.differentials, x) + @arrayop (i,) curl[i] term=D×x +end +LinearAlgebra.cross(D::ArrayDifferentialOperator, x::Arr) = Arr(D × value(x)) + +function LinearAlgebra.cross(D1::ArrayDifferentialOperator, D2::ArrayDifferentialOperator) + @assert length(D1.vars) == length(D2.vars) == 3 "Cross product is only defined in 3 dimensions." + @assert all(scalarize(isequal.(D1.vars, D2.vars))) "Operators have different variables and cannot be composed." + + ArrayDifferentialOperator(D1.vars, crosscompose(D1.differentials, D2.differentials), "("*D1.name*"×"*D2.name*")") +end + +SymbolicUtils.promote_symtype(::ArrayDifferentialOperator, x) = x + +Base.show(io::IO, D::ArrayDifferentialOperator) = print(io, D.name) +Base.nameof(D::ArrayDifferentialOperator) = Symbol(D.name) + +function Base.:(==)(D1::ArrayDifferentialOperator, D2::ArrayDifferentialOperator) + @variables x[1:length(D1.vars)] + all(scalarize(isequal.(D1.vars, D2.vars))) && all(scalarize(isequal.(D1(x), D2(x)))) +end + +# TODO: Add simplification rules for dot and cross products to remove 0 terms and simplify \ No newline at end of file diff --git a/test/diff.jl b/test/diff.jl index 7486a78e7..532f4cd5d 100644 --- a/test/diff.jl +++ b/test/diff.jl @@ -1,7 +1,7 @@ using Symbolics using Test using IfElse -using Symbolics: value +using Symbolics: value, scalarize, Nabla, Div, Curl, Laplacian # Derivatives @variables t σ ρ β @@ -37,32 +37,32 @@ dcsch = D(csch(t)) # Chain rule dsinsin = D(sin(sin(t))) -test_equal(expand_derivatives(dsinsin), cos(sin(t))*cos(t)) +test_equal(expand_derivatives(dsinsin), cos(sin(t)) * cos(t)) -d1 = D(sin(t)*t) -d2 = D(sin(t)*cos(t)) -@test isequal(expand_derivatives(d1), simplify(t*cos(t)+sin(t))) -@test isequal(expand_derivatives(d2), cos(t)^2-sin(t)^2) +d1 = D(sin(t) * t) +d2 = D(sin(t) * cos(t)) +@test isequal(expand_derivatives(d1), simplify(t * cos(t) + sin(t))) +@test isequal(expand_derivatives(d2), cos(t)^2 - sin(t)^2) -eqs = [σ*(y-x), - x*(ρ-z)-y, - x*y - β*z] +eqs = [σ * (y - x), + x * (ρ - z) - y, + x * y - β * z] jac = Symbolics.jacobian(eqs, [x, y, z]) -test_equal(jac[1,1], -1σ) -test_equal(jac[1,2], σ) -test_equal(jac[1,3], 0) -test_equal(jac[2,1], ρ - z) -test_equal(jac[2,2], -1) -test_equal(jac[2,3], -1x) -test_equal(jac[3,1], y) -test_equal(jac[3,2], x) -test_equal(jac[3,3], -1β) +test_equal(jac[1, 1], -1σ) +test_equal(jac[1, 2], σ) +test_equal(jac[1, 3], 0) +test_equal(jac[2, 1], ρ - z) +test_equal(jac[2, 2], -1) +test_equal(jac[2, 3], -1x) +test_equal(jac[3, 1], y) +test_equal(jac[3, 2], x) +test_equal(jac[3, 3], -1β) # issue #545 z = t + t^2 #test_equal(expand_derivatives(D(z)), 1 + t * 2) -z = t-2t +z = t - 2t #test_equal(expand_derivatives(D(z)), -1) # Variable dependence checking in differentiation @@ -73,8 +73,8 @@ z = t-2t @variables x(t) y(t) z(t) -@test isequal(expand_derivatives(D(x * y)), simplify(y*D(x) + x*D(y))) -@test isequal(expand_derivatives(D(x * y)), simplify(D(x)*y + x*D(y))) +@test isequal(expand_derivatives(D(x * y)), simplify(y * D(x) + x * D(y))) +@test isequal(expand_derivatives(D(x * y)), simplify(D(x) * y + x * D(y))) @test isequal(expand_derivatives(D(2t)), 2) @test isequal(expand_derivatives(D(2x)), 2D(x)) @@ -88,16 +88,16 @@ z = t-2t @test all(iszero, Symbolics.gradient(42, [t, x, y, z])) @test all(iszero, Symbolics.hessian(42, [t, x, y, z])) @test isequal(Symbolics.jacobian([t, x, 42], [t, x]), - Num[1 0 - Differential(t)(x) 1 - 0 0]) + Num[1 0 + Differential(t)(x) 1 + 0 0]) # issue 252 @variables beta, alpha, delta @variables x1, x2, x3 # expression -tmp = beta * (alpha * exp(x1) * x2 ^ (alpha - 1) + 1 - delta) / x3 +tmp = beta * (alpha * exp(x1) * x2^(alpha - 1) + 1 - delta) / x3 # derivative w.r.t. x1 and x2 t1 = Symbolics.gradient(tmp, [x1, x2]) @test_nowarn Symbolics.gradient(t1[1], [beta]) @@ -111,7 +111,7 @@ using Symbolics @variables t x(t) ∂ₜ = Differential(t) ∂ₓ = Differential(x) -L = .5 * ∂ₜ(x)^2 - .5 * x^2 +L = 0.5 * ∂ₜ(x)^2 - 0.5 * x^2 @test isequal(expand_derivatives(∂ₓ(L)), -1 * x) @test isequal(expand_derivatives(Differential(x)(L) - ∂ₜ(Differential(∂ₜ(x))(L))), -1 * (∂ₜ(∂ₜ(x)) + x)) @test isequal(expand_derivatives(Differential(x)(L) - ∂ₜ(Differential(∂ₜ(x))(L))), (-1 * x) - ∂ₜ(∂ₜ(x))) @@ -123,9 +123,9 @@ L = .5 * ∂ₜ(x)^2 - .5 * x^2 @variables u(..) Dy = Differential(y) Dx = Differential(x) -dxyu = Dx(Dy(u(x,y))) +dxyu = Dx(Dy(u(x, y))) @test isequal(expand_derivatives(dxyu), dxyu) -dxxu = Dx(Dx(u(x,y))) +dxxu = Dx(Dx(u(x, y))) @test isequal(expand_derivatives(dxxu), dxxu) using Symbolics, LinearAlgebra, SparseArrays @@ -137,112 +137,114 @@ canonequal(a, b) = isequal(simplify(a), simplify(b)) @variables t σ ρ β @variables x y z @test isequal( - (Differential(z) * Differential(y) * Differential(x))(t), - Differential(z)(Differential(y)(Differential(x)(t))) + (Differential(z) * Differential(y) * Differential(x))(t), + Differential(z)(Differential(y)(Differential(x)(t))), ) @test canonequal( - Symbolics.derivative(sin(cos(x)), x), - -sin(x) * cos(cos(x)) - ) + Symbolics.derivative(sin(cos(x)), x), + -sin(x) * cos(cos(x)), +) Symbolics.@register_symbolic no_der(x) @test canonequal( - Symbolics.derivative([sin(cos(x)), hypot(x, no_der(x))], x), - [ - -sin(x) * cos(cos(x)), - x/hypot(x, no_der(x)) + - no_der(x)*Differential(x)(no_der(x))/hypot(x, no_der(x)) - ] - ) + Symbolics.derivative([sin(cos(x)), hypot(x, no_der(x))], x), + [ + -sin(x) * cos(cos(x)), + x / hypot(x, no_der(x)) + + no_der(x) * Differential(x)(no_der(x)) / hypot(x, no_der(x)), + ], +) Symbolics.@register_symbolic intfun(x)::Int @test Symbolics.symtype(intfun(x).val) === Int -eqs = [σ*(y-x), - x*(ρ-z)-y, - x*y - β*z] +eqs = [σ * (y - x), + x * (ρ - z) - y, + x * y - β * z] -∂ = Symbolics.jacobian(eqs,[x,y,z]) +∂ = Symbolics.jacobian(eqs, [x, y, z]) for i in 1:3 - ∇ = Symbolics.gradient(eqs[i],[x,y,z]) - @test canonequal(∂[i,:],∇) + ∇ = Symbolics.gradient(eqs[i], [x, y, z]) + @test canonequal(∂[i, :], ∇) end -@test all(canonequal.(Symbolics.gradient(eqs[1],[x,y,z]),[σ * -1,σ,0])) -@test all(canonequal.(Symbolics.hessian(eqs[1],[x,y,z]),0)) +@test all(canonequal.(Symbolics.gradient(eqs[1], [x, y, z]), [σ * -1, σ, 0])) +@test all(canonequal.(Symbolics.hessian(eqs[1], [x, y, z]), 0)) -du = [x^2, y^3, x^4, sin(y), x+y, x+z^2, z+x, x+y^2+sin(z)] -reference_jac = sparse(Symbolics.jacobian(du, [x,y,z])) +du = [x^2, y^3, x^4, sin(y), x + y, x + z^2, z + x, x + y^2 + sin(z)] +reference_jac = sparse(Symbolics.jacobian(du, [x, y, z])) -@test findnz(Symbolics.jacobian_sparsity(du, [x,y,z]))[[1,2]] == findnz(reference_jac)[[1,2]] +@test findnz(Symbolics.jacobian_sparsity(du, [x, y, z]))[[1, 2]] == findnz(reference_jac)[[1, 2]] let - @variables t x(t) y(t) z(t) - @test Symbolics.exprs_occur_in([x,y,z], x^2*y) == [true, true, false] + @variables t x(t) y(t) z(t) + @test Symbolics.exprs_occur_in([x, y, z], x^2 * y) == [true, true, false] end -@test isequal(Symbolics.sparsejacobian(du, [x,y,z]), reference_jac) +@test isequal(Symbolics.sparsejacobian(du, [x, y, z]), reference_jac) -function f!(res,u) - (x,y,z)=u - res.=[x^2, y^3, x^4, sin(y), x+y, x+z^2, z+x, x+y^2+sin(z)] +function f!(res, u) + (x, y, z) = u + res .= [x^2, y^3, x^4, sin(y), x + y, x + z^2, z + x, x + y^2 + sin(z)] end -function f1!(res,u,a,b;c) - (x,y,z)=u - res.=[a*x^2, y^3, b*x^4, sin(y), c*x+y, x+z^2, a*z+x, x+y^2+sin(z)] +function f1!(res, u, a, b; c) + (x, y, z) = u + res .= [a * x^2, y^3, b * x^4, sin(y), c * x + y, x + z^2, a * z + x, x + y^2 + sin(z)] end -input=rand(3) -output=rand(8) +input = rand(3) +output = rand(8) -findnz(Symbolics.jacobian_sparsity(f!, output, input))[[1,2]] == findnz(reference_jac)[[1,2]] -findnz(Symbolics.jacobian_sparsity(f1!, output, input,1,2,c=3))[[1,2]] == findnz(reference_jac)[[1,2]] +findnz(Symbolics.jacobian_sparsity(f!, output, input))[[1, 2]] == findnz(reference_jac)[[1, 2]] +findnz(Symbolics.jacobian_sparsity(f1!, output, input, 1, 2, c = 3))[[1, 2]] == findnz(reference_jac)[[1, 2]] -input = rand(2,2) -function f2!(res,u,a,b,c) - (x,y,z)=u[1,1],u[2,1],u[3,1] - res.=[a*x^2, y^3, b*x^4, sin(y), c*x+y, x+z^2, a*z+x, x+y^2+sin(z)] +input = rand(2, 2) +function f2!(res, u, a, b, c) + (x, y, z) = u[1, 1], u[2, 1], u[3, 1] + res .= [a * x^2, y^3, b * x^4, sin(y), c * x + y, x + z^2, a * z + x, x + y^2 + sin(z)] end -findnz(Symbolics.jacobian_sparsity(f!, output, input))[[1,2]] == findnz(reference_jac)[[1,2]] +findnz(Symbolics.jacobian_sparsity(f!, output, input))[[1, 2]] == findnz(reference_jac)[[1, 2]] # Check for failures due to du[4] undefined -function f_undef(du,u) - du[1] = u[1] - du[2] = u[2] - du[3] = u[3] + u[4] +function f_undef(du, u) + du[1] = u[1] + du[2] = u[2] + du[3] = u[3] + u[4] end u0 = rand(4) du0 = similar(u0) -sparsity_pattern = Symbolics.jacobian_sparsity(f_undef,du0,u0) +sparsity_pattern = Symbolics.jacobian_sparsity(f_undef, du0, u0) udef_ref = sparse([1 0 0 0 - 0 1 0 0 - 0 0 1 1 - 0 0 0 0]) -findnz(sparsity_pattern)[[1,2]] == findnz(udef_ref)[[1,2]] + 0 1 0 0 + 0 0 1 1 + 0 0 0 0]) +findnz(sparsity_pattern)[[1, 2]] == findnz(udef_ref)[[1, 2]] using Symbolics -rosenbrock(X) = sum(1:length(X)-1) do i - 100 * (X[i+1] - X[i]^2)^2 + (1 - X[i])^2 -end -rosenbrock2(X,a;b) = sum(1:length(X)-1) do i - a * (X[i+1] - X[i]^2)^2 + (b - X[i])^2 -end - -@variables a,b -X = [a,b] +rosenbrock(X) = + sum(1:length(X)-1) do i + 100 * (X[i+1] - X[i]^2)^2 + (1 - X[i])^2 + end +rosenbrock2(X, a; b) = + sum(1:length(X)-1) do i + a * (X[i+1] - X[i]^2)^2 + (b - X[i])^2 + end + +@variables a, b +X = [a, b] input = rand(2) -spoly(x) = simplify(x, expand=true) +spoly(x) = simplify(x, expand = true) rr = rosenbrock(X) reference_hes = Symbolics.hessian(rr, X) @test findnz(sparse(reference_hes))[1:2] == findnz(Symbolics.hessian_sparsity(rr, X))[1:2] @test findnz(sparse(reference_hes))[1:2] == findnz(Symbolics.hessian_sparsity(rosenbrock, input))[1:2] -@test findnz(sparse(reference_hes))[1:2] == findnz(Symbolics.hessian_sparsity(rosenbrock2, input,100,b=1))[1:2] +@test findnz(sparse(reference_hes))[1:2] == findnz(Symbolics.hessian_sparsity(rosenbrock2, input, 100, b = 1))[1:2] sp_hess = Symbolics.sparsehessian(rr, X) @test findnz(sparse(reference_hes))[1:2] == findnz(sp_hess)[1:2] @@ -252,24 +254,24 @@ sp_hess = Symbolics.sparsehessian(rr, X) @variables t x(t)[1:4] ẋ(t)[1:4] expression = sin(x[1] + x[2] + x[3] + x[4]) |> Differential(t) |> expand_derivatives expression2 = substitute(expression, Dict(collect(Differential(t).(x) .=> ẋ))) -@test isequal(expression2, (ẋ[1] + ẋ[2] + ẋ[3] + ẋ[4])*cos(x[1] + x[2] + x[3] + x[4])) +@test isequal(expression2, (ẋ[1] + ẋ[2] + ẋ[3] + ẋ[4]) * cos(x[1] + x[2] + x[3] + x[4])) @test isequal( - Symbolics.derivative(IfElse.ifelse(signbit(b), b^2, sqrt(b)), b), - IfElse.ifelse(signbit(b), 2b,(SymbolicUtils.unstable_pow(2Symbolics.unwrap(sqrt(b)), -1))) + Symbolics.derivative(IfElse.ifelse(signbit(b), b^2, sqrt(b)), b), + IfElse.ifelse(signbit(b), 2b, (SymbolicUtils.unstable_pow(2Symbolics.unwrap(sqrt(b)), -1))), ) # Chain rule # let - @variables t x(..) y(..) z(..) - Dt = Differential(t) + @variables t x(..) y(..) z(..) + Dt = Differential(t) - @test isequal(expand_derivatives(Dt(x(y(t)))), - Dt(y(t))*Differential(y(t))(x(y(t)))) + @test isequal(expand_derivatives(Dt(x(y(t)))), + Dt(y(t)) * Differential(y(t))(x(y(t)))) - @test isequal(expand_derivatives(Dt(x(y(t), z(t)))), - Dt(y(t))*Differential(y(t))(x(y(t), z(t))) + Dt(z(t))*Differential(z(t))(x(y(t), z(t)))) + @test isequal(expand_derivatives(Dt(x(y(t), z(t)))), + Dt(y(t)) * Differential(y(t))(x(y(t), z(t))) + Dt(z(t)) * Differential(z(t))(x(y(t), z(t)))) end @variables x y @@ -282,69 +284,120 @@ D = Differential(t) eqs = [D(x) ~ x, D(y) ~ y + x] -sub_eqs = substitute(eqs, Dict([D(x)=>D(x), x=>1])) +sub_eqs = substitute(eqs, Dict([D(x) => D(x), x => 1])) @test sub_eqs == [D(x) ~ 1, D(y) ~ 1 + y] @variables x y -@test substitute([x + y; x - y], Dict(x=>1, y=>2)) == [3, -1] +@test substitute([x + y; x - y], Dict(x => 1, y => 2)) == [3, -1] # 530#discussion_r825125589 let - using Symbolics - @variables u[1:2] y[1:1] t - u = collect(u) - y = collect(y) - @test isequal(Symbolics.jacobian([u;u[1]^2; y], u), Num[1 0 - 0 1 - 2u[1] 0 - 0 0]) + using Symbolics + @variables u[1:2] y[1:1] t + u = collect(u) + y = collect(y) + @test isequal( + Symbolics.jacobian([u; u[1]^2; y], u), + Num[1 0 + 0 1 + 2u[1] 0 + 0 0], + ) end # make sure derivative(x[1](t), y) does not fail let - @variables t a(t) - vars = collect(@variables(x(t)[1:1])[1]) - ps = collect(@variables(ps[1:1])[1]) - @test Symbolics.derivative(ps[1], vars[1]) == 0 - @test Symbolics.derivative(ps[1], a) == 0 - @test Symbolics.derivative(x[1], a) == 0 + @variables t a(t) + vars = collect(@variables(x(t)[1:1])[1]) + ps = collect(@variables(ps[1:1])[1]) + @test Symbolics.derivative(ps[1], vars[1]) == 0 + @test Symbolics.derivative(ps[1], a) == 0 + @test Symbolics.derivative(x[1], a) == 0 end # 580 let - @variables x[1:3] - y = [sin.(x); cos.(x)] - dj = Symbolics.jacobian(y, x) - @test !iszero(dj) - @test isequal(dj, Symbolics.jacobian(Symbolics.scalarize.(y), x)) - sj = Symbolics.sparsejacobian(y, x) - @test !iszero(sj) - @test isequal(sj, Symbolics.jacobian(Symbolics.scalarize.(y), x)) + @variables x[1:3] + y = [sin.(x); cos.(x)] + dj = Symbolics.jacobian(y, x) + @test !iszero(dj) + @test isequal(dj, Symbolics.jacobian(Symbolics.scalarize.(y), x)) + sj = Symbolics.sparsejacobian(y, x) + @test !iszero(sj) + @test isequal(sj, Symbolics.jacobian(Symbolics.scalarize.(y), x)) end # substituting iv of differentials @variables t t2 x(t) D = Differential(t) ex = D(x) -ex2 = substitute(ex, [t=>t2]) +ex2 = substitute(ex, [t => t2]) @test isequal(operation(Symbolics.unwrap(ex2)).x, t2) -ex3 = substitute(D(x) * 2 + x / t, [t=>t2]) +ex3 = substitute(D(x) * 2 + x / t, [t => t2]) xt2 = substitute(x, [t => t2]) @test isequal(ex3, xt2 / t2 + 2Differential(t2)(xt2)) # 581 # let - @variables x(t)[1:3] - @test iszero(Symbolics.derivative(x[1], x[2])) + @variables x(t)[1:3] + @test iszero(Symbolics.derivative(x[1], x[2])) end #908 # let - using Symbolics - @variables t - @test isequal(expand_derivatives(Differential(t)(im*t)), im) - @test isequal(expand_derivatives(Differential(t)(t^2 + im*t)), 2t + im) + using Symbolics + @variables t + @test isequal(expand_derivatives(Differential(t)(im * t)), im) + @test isequal(expand_derivatives(Differential(t)(t^2 + im * t)), 2t + im) end + +## Vector Calculus + + +@variables x[1:3] y[1:3] + +∇ = Nabla(x) + +out = ∇(y) +@test all(isequal.(scalarize(out), [Differential(x[1])(y[1]), Differential(x[2])(y[2]), Differential(x[3])(y[3])])) + + +out = ∇ ⋅ y +@test isequal(out, Div(x)(y)) +@test isequal(out, Differential(x[1])(y[1]) + Differential(x[2])(y[2]) + Differential(x[3])(y[3])) + +out = ∇ × y +@test isequal(out, Curl(x)(y)) +@test all(isequal.(scalarize(out), [Differential(x[2])(y[3]) - Differential(x[3])(y[2]), + Differential(x[3])(y[1]) - Differential(x[1])(y[3]), + Differential(x[1])(y[2]) - Differential(x[2])(y[1])])) + + + +out = (∇ ⋅ ∇)(y) +@test all(isequal.(scalarize(out), scalarize(Laplacian(x)(y)))) +@test all( + isequal.( + scalarize(out), + [Differential(x[1])(Differential(x[1])(y[1])) + Differential(x[2])(Differential(x[2])(y[1])) + Differential(x[3])(Differential(x[3])(y[1])), + Differential(x[1])(Differential(x[1])(y[2])) + Differential(x[2])(Differential(x[2])(y[2])) + Differential(x[3])(Differential(x[3])(y[2])), + Differential(x[1])(Differential(x[1])(y[3])) + Differential(x[2])(Differential(x[2])(y[3])) + Differential(x[3])(Differential(x[3])(y[3]))], + ), +) + + +out = (∇ × ∇)(y) + +@test_broken isequal(out, Curl(x)(Curl(x)(y))) +@test all( + isequal.( + scalarize(out), + [Differential(x[2])(Differential(x[3])(y[1])) - Differential(x[3])(Differential(x[2])(y[1])), + Differential(x[3])(Differential(x[1])(y[2])) - Differential(x[1])(Differential(x[3])(y[2])), + Differential(x[1])(Differential(x[2])(y[3])) - Differential(x[2])(Differential(x[1])(y[3]))], + ), +) +