Merge pull request #54 from mikelehu/master

complex equations
SciML · Oct 14, 2023 · 45dfc0d · 45dfc0d
2 parents b10c810 + 3864f50
commit 45dfc0d
Show file tree

Hide file tree

Showing 9 changed files with 2,205 additions and 1,181 deletions.
diff --git a/Examples/.ipynb_checkpoints/Burrau example using Complex Numbers (Github)-checkpoint.ipynb b/Examples/.ipynb_checkpoints/Burrau example using Complex Numbers (Github)-checkpoint.ipynb
diff --git a/Examples/Burrau example using Complex Numbers (Github).ipynb b/Examples/Burrau example using Complex Numbers (Github).ipynb
diff --git a/src/IRKCoefficients.jl b/src/IRKCoefficients.jl
@@ -6,8 +6,8 @@
 #   EstimateCoeffs!
 
 function PolInterp(X::AbstractVector{ctype},
-                   Y::AbstractMatrix{ctype},
-                   Z::AbstractVector{ctype}) where {ctype}
+    Y::AbstractMatrix{ctype},
+    Z::AbstractVector{ctype}) where {ctype}
     N = length(X)
     M = length(Z)
     K = size(Y, 1)
@@ -16,26 +16,28 @@ function PolInterp(X::AbstractVector{ctype},
     end
     pz = zeros(ctype, K, M)
 
-    @inbounds begin for i in 1:N
-        lag = 1.0
-        for j in 1:N
-            if (j != i)
-                lag *= X[i] - X[j]
-            end
-        end
-        lag = 1 / lag
-        for m in 1:M
-            liz = lag
+    @inbounds begin
+        for i in 1:N
+            lag = 1.0
             for j in 1:N
                 if (j != i)
-                    liz *= Z[m] - X[j]
+                    lag *= X[i] - X[j]
                 end
             end
-            for k in 1:K
-                pz[k, m] += Y[k, i] * liz
+            lag = 1 / lag
+            for m in 1:M
+                liz = lag
+                for j in 1:N
+                    if (j != i)
+                        liz *= Z[m] - X[j]
+                    end
+                end
+                for k in 1:K
+                    pz[k, m] += Y[k, i] * liz
+                end
             end
         end
-    end end
+    end
     return pz
 end
 
@@ -85,11 +87,13 @@ function MuCoefficients!(mu, T::Type{<:CompiledFloats})
     mu[8, 8] = convert(T, 0.5)
 
     s = 8
-    @inbounds begin for i in 1:s
-        for j in (i + 1):s
-            mu[i, j] = 1 - mu[j, i]
+    @inbounds begin
+        for i in 1:s
+            for j in (i + 1):s
+                mu[i, j] = 1 - mu[j, i]
+            end
         end
-    end end
+    end
 end
 
 function MuCoefficients!(mu, T)
@@ -141,37 +145,39 @@ function MuCoefficients!(mu, T)
     mu[8, 8] = parse(T, "0.5")
 
     s = 8
-    @inbounds begin for i in 1:s
-        for j in (i + 1):s
-            mu[i, j] = 1 - mu[j, i]
+    @inbounds begin
+        for i in 1:s
+            for j in (i + 1):s
+                mu[i, j] = 1 - mu[j, i]
+            end
         end
-    end end
+    end
 end
 
 function HCoefficients!(mu, hc, hb, nu, h, hprev, T::Type{<:CompiledFloats})
     hb .= convert.(T,
-                   [
-                       +5.0614268145188129576265677154981094e-02,
-                       +1.1119051722668723527217799721312045e-01,
-                       +1.5685332293894364366898110099330067e-01,
-                       +1.8134189168918099148257522463859781e-01,
-                       +1.8134189168918099148257522463859781e-01,
-                       +1.5685332293894364366898110099330067e-01,
-                       +1.1119051722668723527217799721312045e-01,
-                       +5.0614268145188129576265677154981094e-02,
-                   ])
+        [
+            +5.0614268145188129576265677154981094e-02,
+            +1.1119051722668723527217799721312045e-01,
+            +1.5685332293894364366898110099330067e-01,
+            +1.8134189168918099148257522463859781e-01,
+            +1.8134189168918099148257522463859781e-01,
+            +1.5685332293894364366898110099330067e-01,
+            +1.1119051722668723527217799721312045e-01,
+            +5.0614268145188129576265677154981094e-02,
+        ])
 
     hc .= convert.(T,
-                   [
-                       +1.9855071751231884158219565715263505e-02,
-                       +1.0166676129318663020422303176208480e-01,
-                       +2.3723379504183550709113047540537686e-01,
-                       +4.0828267875217509753026192881990801e-01,
-                       +5.9171732124782490246973807118009203e-01,
-                       +7.6276620495816449290886952459462321e-01,
-                       +8.9833323870681336979577696823791522e-01,
-                       +9.8014492824876811584178043428473653e-01,
-                   ])
+        [
+            +1.9855071751231884158219565715263505e-02,
+            +1.0166676129318663020422303176208480e-01,
+            +2.3723379504183550709113047540537686e-01,
+            +4.0828267875217509753026192881990801e-01,
+            +5.9171732124782490246973807118009203e-01,
+            +7.6276620495816449290886952459462321e-01,
+            +8.9833323870681336979577696823791522e-01,
+            +9.8014492824876811584178043428473653e-01,
+        ])
 
     s = length(hb)
 
@@ -253,16 +259,16 @@ end
 
 function EstimateCoeffs!(alpha, T::Type{<:CompiledFloats})
     c = convert.(T,
-                 [
-                     +1.9855071751231884158219565715263505e-02,
-                     +1.0166676129318663020422303176208480e-01,
-                     +2.3723379504183550709113047540537686e-01,
-                     +4.0828267875217509753026192881990801e-01,
-                     +5.9171732124782490246973807118009203e-01,
-                     +7.6276620495816449290886952459462321e-01,
-                     +8.9833323870681336979577696823791522e-01,
-                     +9.8014492824876811584178043428473653e-01,
-                 ])
+        [
+            +1.9855071751231884158219565715263505e-02,
+            +1.0166676129318663020422303176208480e-01,
+            +2.3723379504183550709113047540537686e-01,
+            +4.0828267875217509753026192881990801e-01,
+            +5.9171732124782490246973807118009203e-01,
+            +7.6276620495816449290886952459462321e-01,
+            +8.9833323870681336979577696823791522e-01,
+            +9.8014492824876811584178043428473653e-01,
+        ])
 
     s = length(c)
 

diff --git a/src/IRKGL16AuxFunctions.jl b/src/IRKGL16AuxFunctions.jl
@@ -6,41 +6,57 @@
 
 function ErrorEst(U, F, dt, beta, abstol, reltol)
     uiType = eltype(U[1])
+    realuiType = real(uiType)
 
     (s,) = size(F)
     D = length(U[1])
 
-    est = zero(uiType)
+    est = zero(realuiType)
 
-    @inbounds begin for k in eachindex(U[1])
-        sum = zero(uiType)
-        maxU = zero(uiType)
-        for is in 1:s
-            sum += beta[is] * F[is][k]
-            maxU = max(maxU, abs(U[is][k]))
-        end
+    @inbounds begin
+        for k in eachindex(U[1])
+            sum = zero(realuiType)
+            maxU = zero(realuiType)
+            for is in 1:s
+                sum += beta[is] * F[is][k]
+                maxU = max(maxU, abs(U[is][k]))
+            end
 
-        est += (abs(dt * sum))^2 / (abstol + maxU^2 * reltol)
-    end end
+            est += (abs(dt * sum))^2 / (abstol + maxU^2 * reltol)
+        end
+    end
 
     return (est / D)
 end
 
 function MyNorm(u, abstol, reltol)
     uiType = eltype(u)
-    norm = zero(uiType)
 
-    @inbounds begin for k in eachindex(u)
-        aux = u[k] / (abstol + abs(u[k]) * reltol)
-        norm += aux * aux
-    end end
+    realuiType = real(uiType)
+    norm = zero(realuiType)
+
+    if uiType <: Complex
+        @inbounds begin
+            for k in eachindex(u)
+                aux = u[k] / (abstol + abs(u[k]) * reltol)
+                norm += aux * conj(aux)
+            end
+        end
+    else
+        @inbounds begin
+            for k in eachindex(u)
+                aux = u[k] / (abstol + abs(u[k]) * reltol)
+                norm += aux * aux
+            end
+        end
+    end
 
     norm = sqrt(norm / (length(u)))
 
-    return (norm)
+    return (real(norm))
 end
 
-function Rdigits(x::Real, r::Integer)
+function Rdigits(x, r::Integer)
     mx = r * x
     mxx = mx + x
     return (mxx - mx)