These are chat archives for symengine/symengine

25th
Jan 2017
Isuru Fernando
@isuruf
Jan 25 2017 03:05
@mkarakoc, not yet
Christopher Rackauckas
@ChrisRackauckas
Jan 25 2017 03:09
Isuru Fernando
@isuruf
Jan 25 2017 03:17
@ChrisRackauckas, does Mathematica succeed for complex examples in calculating expm?
Christopher Rackauckas
@ChrisRackauckas
Jan 25 2017 03:17
For complex examples?
Haven't tried anything too complex
I am looking to use it for enhancing DifferentialEquations.jl's stiff solvers (for certain PDEs) though
so I haven't used Mathematica for that too much
Isuru Fernando
@isuruf
Jan 25 2017 03:20
Examples you sent were all polynomials with integer coefficients or polynomials with order <=4. This will not work for higher order polynomials with symbolic coefficients
Even sympy can only handle polynomials with integer/real coefficients or polynomials with order <=4
Christopher Rackauckas
@ChrisRackauckas
Jan 25 2017 03:27
For expm?
Mathematica worked with whatever I've thrown at it
the polynomials from DiffEqs tend to be very low order
like no more than 2
I'm looking to expm the Jacobian of an ODE system
Christopher Rackauckas
@ChrisRackauckas
Jan 25 2017 03:35
If you're curious, it to symbolically pre-calculate expressions for $e^{hL}$ where $L$ is the system Jacobian (for part of the system) for this type of method: https://en.wikipedia.org/wiki/Exponential_integrator
Isuru Fernando
@isuruf
Jan 25 2017 03:39
I meant the characteristic polynomial of the matrix
Christopher Rackauckas
@ChrisRackauckas
Jan 25 2017 03:40
oh hmm
Isuru Fernando
@isuruf
Jan 25 2017 03:40
But if the polynomials from DiffEqs are real/integer coefficient ones, that'll not be too hard to solve
Christopher Rackauckas
@ChrisRackauckas
Jan 25 2017 03:42
Here's a pretty standard equation
More complicated forms of things like that
but yeah, polynomials with real coefficients
that's super standard and is what I'd like to do some pre-calculations on
I know I won't be able to handle all cases, but the most standard cases will be pretty vital
Isuru Fernando
@isuruf
Jan 25 2017 03:44
In that example, all the coefficients are integers. Do you also get polynomials with float coefficients
?
Christopher Rackauckas
@ChrisRackauckas
Jan 25 2017 03:47
yes
in that example one of the coefficients is a=1.5
Isuru Fernando
@isuruf
Jan 25 2017 03:48
Okay, so a, b, c, d are not symbolic
Christopher Rackauckas
@ChrisRackauckas
Jan 25 2017 03:48
yes
they can be treated as constants
real constants
and the symbols are x and y
the right hand side is are the equations
derivatives in x and y makes the Jacobian
(each row is for each function)
and then I'd like to expm that
Isuru Fernando
@isuruf
Jan 25 2017 03:54
we could use Nemo to solve the polynomials
but Nemo is a pain to install
Christopher Rackauckas
@ChrisRackauckas
Jan 25 2017 03:55
yeah
I don't want a Nemo dependency
I can wait quite awhile: this is just for speedup a method, and that method is niche enough I don't plan to implement it for awhile
Isuru Fernando
@isuruf
Jan 25 2017 04:29
julia> using SymEngine

julia> import Base.expm

julia> function roots(A::Array{SymEngine.Basic,2})
b = A[1,1] + A[2,2]
c = A[1,1]*A[2,2]-A[1,2]*A[2,1]
l1 = (b + sqrt(b^2-4*c))/2
l2 = (b - sqrt(b^2-4*c))/2
return l1, l2
end
roots (generic function with 1 method)

julia> function expm(A::Array{SymEngine.Basic,2})
l1, l2 = roots(A)
P = SymEngine.Basic[A[1,2] A[1,2]; l1-A[1,1] l2-A[1,1]]
D = SymEngine.Basic[exp(l1) 0; 0 exp(l2)]
P*D*inv(P)
end
expm (generic function with 7 methods)

julia> A = SymEngine.Basic[a - b*y -b*x; y -c + x]
ERROR: UndefVarError: a not defined

julia> @vars a b c y x
(a,b,c,y,x)

julia> A = SymEngine.Basic[a - b*y -b*x; y -c + x]
2×2 Array{SymEngine.Basic,2}:
a - b*y    -b*x
y  -c + x

julia> expm(A)
@ChrisRackauckas, that method is for 2 x 2 matrices
result is pretty large
It'll get huge when the expressionss in A becomes larger
Christopher Rackauckas
@ChrisRackauckas
Jan 25 2017 04:31
what's it look like?
Isuru Fernando
@isuruf
Jan 25 2017 04:31
julia> expm(A)
2×2 Array{SymEngine.Basic,2}:
E^((1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))*(1 + (-(a - b*y) + (1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))/(-(a - b*y) - (-(a - b*y) + (1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2))) + (1/2)*(a - c + x - b*y - sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))) - E^((1/2)*(a - c + x - b*y - sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))*(-(a - b*y) + (1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))/(-(a - b*y) - (-(a - b*y) + (1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2))) + (1/2)*(a - c + x - b*y - sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))  …                                                                                                                                                                                                   -b*x*E^((1/2)*(a - c + x - b*y - sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))/(-(a - b*y) - (-(a - b*y) + (1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2))) + (1/2)*(a - c + x - b*y - sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2))) + b*x*E^((1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))/(-(a - b*y) - (-(a - b*y) + (1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2))) + (1/2)*(a - c + x - b*y - sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))
-E^((1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))*(-(a - b*y) + (1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))*(1 + (-(a - b*y) + (1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))/(-(a - b*y) - (-(a - b*y) + (1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2))) + (1/2)*(a - c + x - b*y - sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2))))/(b*x) + E^((1/2)*(a - c + x - b*y - sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))*(-(a - b*y) + (1/2)*(a - c + x - b*y - sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))*(-(a - b*y) + (1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))/(b*x*(-(a - b*y) - (-(a - b*y) + (1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2))) + (1/2)*(a - c + x - b*y - sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2))))     E^((1/2)*(a - c + x - b*y - sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))*(-(a - b*y) + (1/2)*(a - c + x - b*y - sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))/(-(a - b*y) - (-(a - b*y) + (1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2))) + (1/2)*(a - c + x - b*y - sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2))) - E^((1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))*(-(a - b*y) + (1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))/(-(a - b*y) - (-(a - b*y) + (1/2)*(a - c + x - b*y + sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2))) + (1/2)*(a - c + x - b*y - sqrt(-4*((-c + x)*(a - b*y) + b*x*y) + (a - c + x - b*y)^2)))
Christopher Rackauckas
@ChrisRackauckas
Jan 25 2017 04:31
I see
:worried:
Good to know
looks like it's impractical then
that's the end of the road there then :smile:
Looks like I just need a good expmv implementation then
Isuru Fernando
@isuruf
Jan 25 2017 05:02
It's still faster than expm though
julia> @time expm(B)
0.000192 seconds (80 allocations: 5.828 KB)
2×2 Array{Float64,2}:
2.69227   -0.067827
0.339135   0.047011

julia> @time f(0.1, 0.5)
0.000003 seconds (9 allocations: 336 bytes)
2×2 Array{Float64,2}:
2.69227   -0.067827
0.339135   0.047011
f is a function of x and y
Isuru Fernando
@isuruf
Jan 25 2017 05:11
there's numerical error though (0-2 decimals) assuming expm result is correct
Christopher Rackauckas
@ChrisRackauckas
Jan 25 2017 05:20
most systems will be much larger though
so if this case is already getting big, it sounds intractible
Isuru Fernando
@isuruf
Jan 25 2017 05:25
Yes, I agree
Isuru Fernando
@isuruf
Jan 25 2017 07:01
function h(a::Float64, b::Float64, c::Float64, d::Float64)
x0=sqrt(a^2 - 2*a*d + 4*b*c + d^2)
x1=a - d + x0
x2=a/2
x3=d/2
x4=x2 - x3
x5=x0/2
x6=-x5
x7=1/(x4 + x6)
x8=1/x0
x9=-a + d + x0
x10=x2 + x3
x11=exp(x10 + x5)
x12=x11*x8*x9/4
x13=b/(x4 + x5)
x14=exp(x10 + x6)
x15=1/b
x16=x1*x15
x17=x7*x8*x9
x18=x14*(-x1^2*x15*x17/8 - x16/2)
x19=x11*x8*x9/2
x20=x1*x14*x17/4
[-x1*x12*x7-x13*x18 -b*x19*x7+x13*x20; x12*x16+x18 x19-x20]
end

expm2(A::Array{Float64,2}) = h(A[1,1], A[1,2], A[2,1], A[2,2])
@ChrisRackauckas, julia's expm is like 10 times slower for 2x2 matrices
expm2 given above has no precision loss, and performs 10 times faster
Christopher Rackauckas
@ChrisRackauckas
Jan 25 2017 07:03
yeah, but I'd likely not use expm directly
instead an implementation of expmv