Skip to content

Commit

Permalink
symbolic division by h and latexify (#72)
Browse files Browse the repository at this point in the history
* test symbolic division by h and latexify

* symbolic division by h in tutorial

* set version to v0.1.30

* format
  • Loading branch information
ranocha authored Aug 16, 2022
1 parent 0335ad4 commit 15fe575
Show file tree
Hide file tree
Showing 3 changed files with 108 additions and 23 deletions.
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "BSeries"
uuid = "ebb8d67c-85b4-416c-b05f-5f409e808f32"
authors = ["Hendrik Ranocha <[email protected]> and contributors"]
version = "0.1.29"
version = "0.1.30"

[deps]
Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
Expand Down
42 changes: 38 additions & 4 deletions docs/src/tutorials/symbolic_computations.md
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,22 @@ using Latexify
latexify(series, reduce_order_by=1, dt=SymEngine.symbols("h"), cdot=false) |> println
```

The argument `dt=SymEngine.symbols("h")` is not required here.

```@example modified-equation-symengine
using Latexify
latexify(series, reduce_order_by=1, cdot=false) |> println
```

We can obtain basically the same result (up to other notations of the coefficients)
by dividing the B-series `series` by a symbolic variable for the time step size
as follows.

```@example modified-equation-symengine
using Latexify
h = SymEngine.symbols("h")
latexify(series / h, cdot=false) |> println
```

We can also use other packages for the symbolic computations, of course.
SymPy.jl often provides very clean expressions.
Expand All @@ -92,28 +108,46 @@ using Latexify
latexify(series, reduce_order_by=1, dt=SymPy.symbols("h"), cdot=false) |> println
```

We can also use the simplified versions.

```@example modified-equation-sympy
using Latexify
latexify(series, reduce_order_by=1, cdot=false) |> println
```

```@example modified-equation-sympy
using Latexify
h = SymPy.symbols("h", real=true)
latexify(series / h, cdot=false) |> println
```

Alternatively, we can also use Symbolics.jl.

```@example modified-equation-symbolics
using BSeries, Symbolics
Symbolics.@variables α h
Symbolics.@variables α
A = [0 0; 1/(2α) 0]
b = [1-α, α]
c = [0, 1/(2α)]
series = modified_equation(A, b, c, 3)
```

```@example modified-equation-symbolics
using Latexify
Symbolics.@variables h
latexify(series / h, cdot=false) |> println
```

## Working with rational coefficients

High-order Runge-Kutta methods are often given in terms of rational coefficients that have large denominators.
It's often best to use the rational coefficients in calculations, rather than converting to floating-point, since otherwise
It's often best to use the rational coefficients in calculations, rather than converting to floating-point, since otherwise
terms that ought to cancel may not cancel exactly. Rational coefficients can be entered by using `//` for the division operation.

By default, Julia uses 64-bit integers for the numerator and denominator of a rational on 64-bit systems. In practical
calculations with high-order RK methods, the denominators may become too large to be represented with 64 bits,
By default, Julia uses 64-bit integers for the numerator and denominator of a rational on 64-bit systems. In practical
calculations with high-order RK methods, the denominators may become too large to be represented with 64 bits,
leading to overflow. This can be remedied by specifying higher precision when entering the coefficients:

```@example int128-coefficients-symbolics
Expand Down
87 changes: 69 additions & 18 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -31,30 +31,81 @@ using Aqua: Aqua
@test_nowarn latexify(series_integrator, dt = SymEngine.symbols("h"))

@testset "SymEngine.jl" begin
α = SymEngine.symbols("α")
A = [0 0; 1/(2 * α) 0]
b = [1 - α, α]
c = [0, 1 / (2 * α)]
series_integrator = bseries(A, b, c, 3)
@test_nowarn latexify(series_integrator)
@testset "Symbolic coefficients" begin
α = SymEngine.symbols("α")
A = [0 0; 1/(2 * α) 0]
b = [1 - α, α]
c = [0, 1 / (2 * α)]
series_integrator = bseries(A, b, c, 3)
@test_nowarn latexify(series_integrator)
end

@testset "Divide by h" begin
A = [0 0; 1//2 0]
b = [0, 1]
c = [0, 1 // 2]
series_integrator = bseries(A, b, c, 1)
@test_nowarn latexify(series_integrator)

h = SymEngine.symbols("h")
coefficients = modified_equation(series_integrator)
val1 = @test_nowarn latexify(coefficients, reduce_order_by = 1,
cdot = false)
val2 = @test_nowarn latexify(coefficients / h, cdot = false)
@test val1 == val2
end
end

@testset "SymPy.jl" begin
α = SymPy.symbols("α", real = true)
A = [0 0; 1/(2 * α) 0]
b = [1 - α, α]
c = [0, 1 / (2 * α)]
series_integrator = bseries(A, b, c, 3)
@test_nowarn latexify(series_integrator)
@testset "Symbolic coefficients" begin
α = SymPy.symbols("α", real = true)
A = [0 0; 1/(2 * α) 0]
b = [1 - α, α]
c = [0, 1 / (2 * α)]
series_integrator = bseries(A, b, c, 3)
@test_nowarn latexify(series_integrator)
end

@testset "Divide by h" begin
A = [0 0; 1//2 0]
b = [0, 1]
c = [0, 1 // 2]
series_integrator = bseries(A, b, c, 1)
@test_nowarn latexify(series_integrator)

h = SymPy.symbols("h", real = true)
coefficients = modified_equation(series_integrator)
val1 = @test_nowarn latexify(coefficients, reduce_order_by = 1,
cdot = false)
val2 = @test_nowarn latexify(coefficients / h, cdot = false)
@test val1 == val2
end
end

@testset "Symbolics.jl" begin
Symbolics.@variables α
A = [0 0; 1/(2 * α) 0]
b = [1 - α, α]
c = [0, 1 / (2 * α)]
series_integrator = bseries(A, b, c, 3)
@test_nowarn latexify(series_integrator)
@testset "Symbolic coefficients" begin
Symbolics.@variables α
A = [0 0; 1/(2 * α) 0]
b = [1 - α, α]
c = [0, 1 / (2 * α)]
series_integrator = bseries(A, b, c, 3)
@test_nowarn latexify(series_integrator)
end

@testset "Divide by h" begin
A = [0 0; 1//2 0]
b = [0, 1]
c = [0, 1 // 2]
series_integrator = bseries(A, b, c, 1)
@test_nowarn latexify(series_integrator)

Symbolics.@variables h
coefficients = modified_equation(series_integrator)
val1 = @test_nowarn latexify(coefficients, reduce_order_by = 1,
cdot = false)
val2 = @test_nowarn latexify(coefficients / h, cdot = false)
@test val1 == val2
end
end
end # @testset "latexify"

Expand Down

2 comments on commit 15fe575

@ranocha
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/66358

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.1.30 -m "<description of version>" 15fe5759799f47e4b3c04f07f18952af607e79f1
git push origin v0.1.30

Please sign in to comment.