Computing a modified equation with BigInt

[ketch]

I'm posting this issue on behalf of Rob Corless, who wrote to me about it. I don't have time to look into it right now but perhaps @ranocha can see immediately what the problem is.

This works:

using BSeries
using Latexify
import SymPyPythonCall; sp = SymPyPythonCall;
α = sp.symbols("α", real = true)
s = sp.symbols("s", real=true )
As = [0 0; α 0]; bs = [s - s^2/(2*α), s^2/(2*α)]; cs = [0, α]
coeffss = bseries(As, bs, cs, 4)
latexify(coeffss, cdot=false)
huh = modified_equation(coeffss)
latexify(huh)

I can even get a modified equation for that, which tells me how the residual behaves asymptotically.

But the one that I want to do, which even has fewer symbols, fails. Since I am a Julia novice this isn't very surprising. This is the continuous extension of the ode45 code from Matlab. Probably Rational{Int128} is enough, but currently I'm trying BigInt.

bcext = [s - 183//64 * s ^ 2 + 37//12 * s ^ 3 - 145//128 * s ^ 4,0,1500//371 * s ^ 2 - 1000//159 * s ^ 3 + 1000//371 * s ^ 4,-125//32 * s ^ 2 + 125//12 * s ^ 3 - 375//64 * s ^ 4,9477//3392 * s ^ 2 - 729//106 * s ^ 3 + 25515//6784 * s ^ 4,-11//7 * s ^ 2 + 11//3 * s ^ 3 - 55//28 * s ^ 4,3//2 * s ^ 2 - 4 * s ^ 3 + 5//2 * s ^ 4]
A7 = Rational{BigInt}[0 0 0 0 0 0 0; 1//5 0 0 0 0 0 0; 3//40 9//40 0 0 0 0 0; 44//45 -56//15 32//9 0 0 0 0; 19372//6561 -25360//2187 64448//6561 -212//729 0 0 0; 9017//3168 -355//33 46732//5247 49//176 -5103//18656 0 0; 35//384 0 500//1113 125//192 -2187//6784 11//84 0];
c7 = Rational{BigInt}[0, 1//5, 3//10, 4//5, 8//9,1,1];
coeffs7 = bseries(A7, bcext, c7, 6);
latexify( coeffs7, cdot=false )

yields an error message:

MethodError: Cannot convert an object of type Rational{BigInt} to an object of type PythonCall.Py

[ranocha]

I cannot reproduce the error locally:

julia> using BSeries

julia> using Latexify

julia> import SymPyPythonCall; sp = SymPyPythonCall;

julia> α = sp.symbols("α", real = true)
α

julia> s = sp.symbols("s", real=true )
s

julia> As = [0 0; α 0]; bs = [s - s^2/(2*α), s^2/(2*α)]; cs = [0, α]
2-element Vector{SymPyCore.Sym{PythonCall.Py}}:
 0
 α

julia> coeffss = bseries(As, bs, cs, 4)
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, SymPyCore.Sym{PythonCall.Py}} with 9 entries:
  RootedTree{Int64}: Int64[]      => 1
  RootedTree{Int64}: [1]          => s
  RootedTree{Int64}: [1, 2]       => s^2/2
  RootedTree{Int64}: [1, 2, 3]    => 0
  RootedTree{Int64}: [1, 2, 2]    => s^2*α/2
  RootedTree{Int64}: [1, 2, 3, 4] => 0
  RootedTree{Int64}: [1, 2, 3, 3] => 0
  RootedTree{Int64}: [1, 2, 3, 2] => 0
  RootedTree{Int64}: [1, 2, 2, 2] => s^2*α^2/2

julia> latexify(coeffss, cdot=false)
L"$F_{f}\mathopen{}\left( \varnothing \right)\mathclose{} + h s F_{f}\mathopen{}\left( \rootedtree[.] \right)\mathclose{} + \frac{h^{2} s^{2}}{2} F_{f}\mathopen{}\left( \rootedtree[.[.]] \right)\mathclose{} + \frac{h^{3} s^{2} \alpha}{4} F_{f}\mathopen{}\left( \rootedtree[.[.][.]] \right)\mathclose{} + \frac{h^{4} s^{2} \alpha^{2}}{12} F_{f}\mathopen{}\left( \rootedtree[.[.][.][.]] \right)\mathclose{}$"

julia> huh = modified_equation(coeffss)
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, SymPyCore.Sym{PythonCall.Py}} with 9 entries:
  RootedTree{Int64}: Int64[]      => 0
  RootedTree{Int64}: [1]          => s
  RootedTree{Int64}: [1, 2]       => 0
  RootedTree{Int64}: [1, 2, 3]    => -s^3/6
  RootedTree{Int64}: [1, 2, 2]    => -s^3/3 + s^2*α/2
  RootedTree{Int64}: [1, 2, 3, 4] => s^4/8
  RootedTree{Int64}: [1, 2, 3, 3] => s^4/12 - s*(-s^3/6 + s^2*α/4)
  RootedTree{Int64}: [1, 2, 3, 2] => -s^4/24 - s*(-s^3/3 + s^2*α/2)/2
  RootedTree{Int64}: [1, 2, 2, 2] => -s^4/4 + s^2*α^2/2 - 3*s*(-s^3/3 + s^2*α/2)/2

julia> latexify(huh)
"\$h \\cdot s \\cdot F_{f}\\mathopen{}\\left( \\rootedtree[.] \\right)\\mathclose{} + \\frac{ - h^{3} \\cdot s^{3}}{6} \\cdot F_{f}\\mathopen{}\\left( \\rootedtree[.[.[.]]] \\right)\\mathclose{} + h^{3} \\cdot \\left( \\frac{ - s^{3}}{6} + \\frac{s^{2} \\cdot \\alpha}{4} \\right) \\cdot F_{f}\\mathopen{}\\left( \\rootedtree[.[.][.]] \\right)\\mathclose{} + \\frac{h^{4} " ⋯ 403 bytes ⋯ "{2} \\right)}{2} \\right) \\cdot F_{f}\\mathopen{}\\left( \\rootedtree[.[.[.]][.]] \\right)\\mathclose{} + h^{4} \\cdot \\left( \\frac{ - s^{4}}{24} + \\frac{s^{2} \\cdot \\alpha^{2}}{12} - \\frac{s \\cdot \\left( \\frac{ - s^{3}}{3} + \\frac{s^{2} \\cdot \\alpha}{2} \\right)}{4} \\right) \\cdot F_{f}\\mathopen{}\\left( \\rootedtree[.[.][.][.]] \\right)\\mathclose{}\$"

julia> # next the code that was reported to fail

julia> bcext = [s - 183//64 * s ^ 2 + 37//12 * s ^ 3 - 145//128 * s ^ 4,0,1500//371 * s ^ 2 - 1000//159 * s ^ 3 + 1000//371 * s ^ 4,-125//32 * s ^ 2 + 125//12 * s ^ 3 - 375//64 * s ^ 4,9477//3392 * s ^ 2 - 729//106 * s ^ 3 + 25515//6784 * s ^ 4,-11//7 * s ^ 2 + 11//3 * s ^ 3 - 55//28 * s ^ 4,3//2 * s ^ 2 - 4 * s ^ 3 + 5//2 * s ^ 4]
7-element Vector{SymPyCore.Sym{PythonCall.Py}}:
    -145*s^4/128 + 37*s^3/12 - 183*s^2/64 + s
                                            0
   1000*s^4/371 - 1000*s^3/159 + 1500*s^2/371
        -375*s^4/64 + 125*s^3/12 - 125*s^2/32
 25515*s^4/6784 - 729*s^3/106 + 9477*s^2/3392
             -55*s^4/28 + 11*s^3/3 - 11*s^2/7
                    5*s^4/2 - 4*s^3 + 3*s^2/2

julia> A7 = Rational{BigInt}[0 0 0 0 0 0 0; 1//5 0 0 0 0 0 0; 3//40 9//40 0 0 0 0 0; 44//45 -56//15 32//9 0 0 0 0; 19372//6561 -25360//2187 64448//6561 -212//729 0 0 0; 9017//3168 -355//33 46732//5247 49//176 -5103//18656 0 0; 35//384 0 500//1113 125//192 -2187//6784 11//84 0];

julia> c7 = Rational{BigInt}[0, 1//5, 3//10, 4//5, 8//9,1,1];

julia> coeffs7 = bseries(A7, bcext, c7, 6);

julia> latexify( coeffs7, cdot=false )
"\$F_{f}\\mathopen{}\\left( \\varnothing \\right)\\mathclose{} + h s F_{f}\\mathopen{}\\left( \\rootedtree[.] \\right)\\mathclose{} + \\frac{h^{2} s^{2}}{2} F_{f}\\mathopen{}\\left( \\rootedtree[.[.]] \\right)\\mathclose{} + \\frac{h^{3} s^{3}}{6} F_{f}\\mathopen{}\\left( \\rootedtree[.[.[.]]] \\right)\\mathclose{} + \\frac{h^{3} s^{3}}{6} F_{f}\\mathopen{}\\left( " ⋯ 4788 bytes ⋯ "\\frac{2873 s^{4}}{48600} - \\frac{6089 s^{3}}{97200} + \\frac{3383 s^{2}}{194400} \\right) F_{f}\\mathopen{}\\left( \\rootedtree[.[.[.]][.][.][.]] \\right)\\mathclose{} + h^{6} \\left( \\frac{2873 s^{4}}{486000} - \\frac{6089 s^{3}}{972000} + \\frac{3383 s^{2}}{1944000} \\right) F_{f}\\mathopen{}\\left( \\rootedtree[.[.][.][.][.][.]] \\right)\\mathclose{}\$"

I am using Julia v1.10 and the following package versions:

julia> import Pkg; Pkg.status()
Status `~/.julia/dev/BSeries/run/Project.toml`
  [ebb8d67c] BSeries v0.1.67 `..`
  [23fbe1c1] Latexify v0.16.10
  [295af30f] Revise v3.9.0
  [bc8888f7] SymPyPythonCall v0.5.1

[rcorless]

Thank you very much---I am updating my packages (using Julia v1.9.3, so that may be an issue as well). Had a dirty registry so I had to delete registries/General first, so it's taking it's time in rebuilding All The Things.

I'm mostly happy that the code that I thought should work actually does.

And as a byproduct I have now learned about Issues and reporting on GitHub, which I ought to have known about before.

Nice to meet you, Hendrik. And thank you!

[rcorless]

Updating the packages fixed the problem. For the record, the modified equation indeed has only fifth order corrections (and higher). I don't know if anywhere there exists a publication with this information, so even though it's just an example I will put the results in. Therefore you both will be acknowledged in the paper (I am writing it now with Yalçın Kaya of UniSA in Adelaide.

Thank you!

[ketch]

@rcorless I'm glad that worked and that it has led to something interesting! Will be interested to see the paper.

[ranocha]

Thanks, @rcorless. Please let us know if there are any other issues.