Skip to content

Commit

Permalink
Merge pull request #75 from ErikQQY/name
Browse files Browse the repository at this point in the history 
Rename PIIMTrap and PIIMRect
  • Loading branch information
@ErikQQY
ErikQQY authored Apr 29, 2023
2 parents 9d5ddf0 + 6717074 commit 7eee4aa
Show file tree
Hide file tree
Showing 8 changed files with 64 additions and 69 deletions.
4 changes: 2 additions & 2 deletions docs/src/algorithms.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,8 @@ FractionalDiffEq.FODEMatrixDiscrete
FractionalDiffEq.ClosedForm
FractionalDiffEq.ClosedFormHankelM
FractionalDiffEq.PIPECE
FractionalDiffEq.PIIMRect
FractionalDiffEq.PIIMTrap
FractionalDiffEq.PIRect
FractionalDiffEq.PITrap
```

### System of FODE
Expand Down
2 changes: 1 addition & 1 deletion src/FractionalDiffEq.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@ export FODESolution, FDifferenceSolution, DODESolution, FFMODESolution
export FODESystemSolution, FDDESystemSolution, FFMODESystem

# FODE solvers
export PIEX, PIPECE, PIIMRect, PIIMTrap
export PIEX, PIPECE, PIRect, PITrap
export PECE, FODEMatrixDiscrete, ClosedForm, ClosedFormHankelM, GL
export AtanganaSeda, AtanganaSedaAB
export Euler
Expand Down
21 changes: 10 additions & 11 deletions src/fodesystem/bdf.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,11 +17,10 @@ Use [Backward Differentiation Formula](https://en.wikipedia.org/wiki/Backward_di
"""
struct FLMMBDF <: FODESystemAlgorithm end

function solve(prob::FODESystem, h, ::FLMMBDF; reltol=1e-6, abstol=1e-6)
function solve(prob::FODESystem, h, ::FLMMBDF; reltol=1e-6, abstol=1e-6, maxiters = 100)
@unpack f, α, u0, tspan, p = prob
t0 = tspan[1]; tfinal = tspan[2]
alpha = α[1]
itmax = 100

# issue [#64](https://github.com/SciFracX/FractionalDiffEq.jl/issues/64)
max_order = findmax(α)[1]
Expand Down Expand Up @@ -61,16 +60,16 @@ function solve(prob::FODESystem, h, ::FLMMBDF; reltol=1e-6, abstol=1e-6)
f(temp, u0[:, 1], p, t0)
#fy[:, 1] = BDFf_vectorfield(t0, y0[:, 1], fdefun)
fy[:, 1] = temp
(y, fy) = BDF_first_approximations(t, y, fy, abstol, itmax, s, halpha, omega, w, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = BDF_triangolo(s+1, r-1, 0, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = BDF_first_approximations(t, y, fy, abstol, maxiters, s, halpha, omega, w, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = BDF_triangolo(s+1, r-1, 0, t, y, fy, zn, N, abstol, maxiters, s, w, omega, halpha, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)

# Main process of computation by means of the FFT algorithm
nx0::Int = 0; ny0::Int = 0
ff = zeros(1, 2^(Q+2), 1)
ff[1:2] = [0 2]
for q = 0:Q
L::Int = 2^q
(y, fy) = BDF_disegna_blocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = BDF_disegna_blocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, zn, N, abstol, maxiters, s, w, omega, halpha, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
ff[1:4*L] = [ff[1:2*L]; ff[1:2*L-1]; 4*L]
end
# Evaluation solution in TFINAL when TFINAL is not in the mesh
Expand All @@ -84,7 +83,7 @@ function solve(prob::FODESystem, h, ::FLMMBDF; reltol=1e-6, abstol=1e-6)
end


function BDF_disegna_blocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N , abstol, itmax, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
function BDF_disegna_blocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N , abstol, maxiters, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
nxi::Int = copy(nx0); nxf::Int = copy(nx0 + L*r - 1)
nyi::Int = copy(ny0); nyf::Int = copy(ny0 + L*r - 1)
is::Int = 1
Expand All @@ -97,7 +96,7 @@ function BDF_disegna_blocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N , abstol, i
while ~stop
stop = (nxi+r-1 == nx0+L*r-1) || (nxi+r-1>=Nr-1)
zn = BDF_quadrato(nxi, nxf, nyi, nyf, fy, zn, omega, problem_size)
(y, fy) = BDF_triangolo(nxi, nxi+r-1, nxi, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = BDF_triangolo(nxi, nxi+r-1, nxi, t, y, fy, zn, N, abstol, maxiters, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
i_triangolo = i_triangolo + 1

if ~stop
Expand Down Expand Up @@ -128,7 +127,7 @@ function BDF_quadrato(nxi, nxf, nyi, nyf, fy, zn, omega, problem_size)
return zn
end

function BDF_triangolo(nxi, nxf, j0, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
function BDF_triangolo(nxi, nxf, j0, t, y, fy, zn, N, abstol, maxiters, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
for n = nxi:min(N, nxf)
n1::Int = n+1
St = ABM_starting_term(t[n1], y0, m_alpha, t0, m_alpha_factorial)
Expand Down Expand Up @@ -158,7 +157,7 @@ function BDF_triangolo(nxi, nxf, j0, t, y, fy, zn, N, abstol, itmax, s, w, omega
it = it + 1

stop = (norm(yn1-yn0, Inf) < abstol) || (norm(Gn1, Inf)<abstol)
if it > itmax && ~stop
if it > maxiters && ~stop
@warn "Non Convergence"
stop = true
end
Expand All @@ -175,7 +174,7 @@ function BDF_triangolo(nxi, nxf, j0, t, y, fy, zn, N, abstol, itmax, s, w, omega
return y, fy
end

function BDF_first_approximations(t, y, fy, abstol, itmax, s, halpha, omega, w, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
function BDF_first_approximations(t, y, fy, abstol, maxiters, s, halpha, omega, w, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
m = problem_size
Im = zeros(m, m)+I; Ims = zeros(m*s, m*s)+I
Y0 = zeros(s*m, 1); F0 = copy(Y0); B0 = copy(Y0)
Expand Down Expand Up @@ -219,7 +218,7 @@ function BDF_first_approximations(t, y, fy, abstol, itmax, s, halpha, omega, w,
it = it + 1

stop = (norm(Y1-Y0, Inf) < abstol) || (norm(G1, Inf) < abstol)
if it > itmax && ~stop
if it > maxiters && ~stop
@warn "Non Convergence"
stop = 1
end
Expand Down
21 changes: 10 additions & 11 deletions src/fodesystem/newton_gregory.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,11 +17,10 @@ Use [Newton Gregory](https://www.geeksforgeeks.org/newton-forward-backward-inter
"""
struct FLMMNewtonGregory <: FODESystemAlgorithm end

function solve(prob::FODESystem, h, ::FLMMNewtonGregory; reltol=1e-6, abstol=1e-6)
function solve(prob::FODESystem, h, ::FLMMNewtonGregory; reltol=1e-6, abstol=1e-6, maxiters = 100)
@unpack f, α, u0, tspan, p = prob
t0 = tspan[1]; tfinal = tspan[2]
alpha = α[1]
itmax = 100

# issue [#64](https://github.com/SciFracX/FractionalDiffEq.jl/issues/64)
max_order = findmax(α)[1]
Expand Down Expand Up @@ -63,16 +62,16 @@ function solve(prob::FODESystem, h, ::FLMMNewtonGregory; reltol=1e-6, abstol=1e-
temp = zeros(length(u0[:, 1]))
f(temp, u0[:, 1], p, t0)
fy[:, 1] = temp#f_vectorfield(t0, y0[:, 1], fdefun)
(y, fy) = NG_first_approximations(t, y, fy, abstol, itmax, s, halpha, omega, w, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = NG_triangolo(s+1, r-1, 0, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = NG_first_approximations(t, y, fy, abstol, maxiters, s, halpha, omega, w, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = NG_triangolo(s+1, r-1, 0, t, y, fy, zn, N, abstol, maxiters, s, w, omega, halpha, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)

# Main process of computation by means of the FFT algorithm
nx0 = 0; ny0 = 0
ff = zeros(1, 2^(Q+2), 1)
ff[1:2] = [0 2]
for q = 0:Q
L = Int64(2^q)
(y, fy) = NG_disegna_blocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = NG_disegna_blocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, zn, N, abstol, maxiters, s, w, omega, halpha, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
ff[1:4*L] = [ff[1:2*L]; ff[1:2*L-1]; 4*L]
end
# Evaluation solution in TFINAL when TFINAL is not in the mesh
Expand All @@ -86,7 +85,7 @@ function solve(prob::FODESystem, h, ::FLMMNewtonGregory; reltol=1e-6, abstol=1e-
end


function NG_disegna_blocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N , abstol, itmax, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
function NG_disegna_blocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N , abstol, maxiters, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
nxi::Int = copy(nx0); nxf::Int = copy(nx0 + L*r - 1)
nyi::Int = copy(ny0); nyf::Int = copy(ny0 + L*r - 1)
is::Int = 1
Expand All @@ -99,7 +98,7 @@ function NG_disegna_blocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N , abstol, it
while ~stop
stop = (nxi+r-1 == nx0+L*r-1) || (nxi+r-1>=Nr-1)
zn = NG_quadrato(nxi, nxf, nyi, nyf, fy, zn, omega, problem_size)
(y, fy) = NG_triangolo(nxi, nxi+r-1, nxi, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = NG_triangolo(nxi, nxi+r-1, nxi, t, y, fy, zn, N, abstol, maxiters, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
i_triangolo = i_triangolo + 1

if ~stop
Expand Down Expand Up @@ -130,7 +129,7 @@ function NG_quadrato(nxi, nxf, nyi, nyf, fy, zn, omega, problem_size)
return zn
end

function NG_triangolo(nxi, nxf, j0, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
function NG_triangolo(nxi, nxf, j0, t, y, fy, zn, N, abstol, maxiters, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
for n = nxi:min(N, nxf)
n1::Int = Int64(n+1)
St = NG_starting_term(t[n1], y0, m_alpha, t0, m_alpha_factorial)
Expand Down Expand Up @@ -160,7 +159,7 @@ function NG_triangolo(nxi, nxf, j0, t, y, fy, zn, N, abstol, itmax, s, w, omega,
it = it + 1

stop = (norm(yn1-yn0, Inf) < abstol) || (norm(Gn1, Inf)<abstol)
if it > itmax && ~stop
if it > maxiters && ~stop
@warn "Non Convergence"
stop = true
end
Expand All @@ -177,7 +176,7 @@ function NG_triangolo(nxi, nxf, j0, t, y, fy, zn, N, abstol, itmax, s, w, omega,
return y, fy
end

function NG_first_approximations(t, y, fy, abstol, itmax, s, halpha, omega, w, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
function NG_first_approximations(t, y, fy, abstol, maxiters, s, halpha, omega, w, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
m = problem_size
Im = zeros(m, m)+I ; Ims = zeros(m*s, m*s)+I
Y0 = zeros(s*m, 1); F0 = copy(Y0); B0 = copy(Y0)
Expand Down Expand Up @@ -221,7 +220,7 @@ function NG_first_approximations(t, y, fy, abstol, itmax, s, halpha, omega, w, p
it = it + 1

stop = (norm(Y1-Y0, Inf) < abstol) || (norm(G1, Inf) < abstol)
if it > itmax && ~stop
if it > maxiters && ~stop
@warn "Non Convergence"
stop = 1
end
Expand Down
20 changes: 10 additions & 10 deletions src/fodesystem/trapezoid.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ function solve(prob::FODESystem, h, ::FLMMTrap; reltol=1e-6, abstol=1e-6)
t0 = tspan[1]; tfinal = tspan[2]
u0 = u0
alpha = α[1]
itmax = 100
maxiters = 100

# issue [#64](https://github.com/SciFracX/FractionalDiffEq.jl/issues/64)
max_order = findmax(α)[1]
Expand Down Expand Up @@ -63,16 +63,16 @@ function solve(prob::FODESystem, h, ::FLMMTrap; reltol=1e-6, abstol=1e-6)
temp = zeros(problem_size)
f(temp, u0[:, 1], p, t0)
fy[:, 1] = temp
(y, fy) = TrapFirstApproximations(t, y, fy, abstol, itmax, s, halpha, omega, w, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = TrapTriangolo(s+1, r-1, 0, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = TrapFirstApproximations(t, y, fy, abstol, maxiters, s, halpha, omega, w, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = TrapTriangolo(s+1, r-1, 0, t, y, fy, zn, N, abstol, maxiters, s, w, omega, halpha, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)

# Main process of computation by means of the FFT algorithm
nx0 = 0; ny0 = 0
ff = zeros(1, 2^(Q+2), 1)
ff[1:2] = [0 2]
for q = 0:Q
L::Int = 2^q
(y, fy) = TrapDisegnaBlocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = TrapDisegnaBlocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, zn, N, abstol, maxiters, s, w, omega, halpha, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
ff[1:4*L] = [ff[1:2*L]; ff[1:2*L-1]; 4*L]
end
# Evaluation solution in TFINAL when TFINAL is not in the mesh
Expand All @@ -86,7 +86,7 @@ function solve(prob::FODESystem, h, ::FLMMTrap; reltol=1e-6, abstol=1e-6)
end


function TrapDisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N::Int, abstol, itmax, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
function TrapDisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N::Int, abstol, maxiters, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
nxi::Int = copy(nx0); nxf::Int = copy(nx0 + L*r - 1)
nyi::Int = copy(ny0); nyf::Int = copy(ny0 + L*r - 1)
is::Int = 1
Expand All @@ -99,7 +99,7 @@ function TrapDisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N::Int, abstol
while ~stop
stop = (nxi+r-1 == nx0+L*r-1) || (nxi+r-1>=Nr-1)
zn = TrapQuadrato(nxi, nxf, nyi, nyf, fy, zn, omega, problem_size)
(y, fy) = TrapTriangolo(nxi, nxi+r-1, nxi, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = TrapTriangolo(nxi, nxi+r-1, nxi, t, y, fy, zn, N, abstol, maxiters, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
i_triangolo = i_triangolo + 1

if ~stop
Expand Down Expand Up @@ -130,7 +130,7 @@ function TrapQuadrato(nxi, nxf, nyi, nyf, fy, zn, omega, problem_size)
return zn
end

function TrapTriangolo(nxi, nxf, j0, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
function TrapTriangolo(nxi, nxf, j0, t, y, fy, zn, N, abstol, maxiters, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
for n = nxi:min(N, nxf)
n1::Int = n+1
St = TrapStartingTerm(t[n1], y0, m_alpha, t0, m_alpha_factorial)
Expand Down Expand Up @@ -160,7 +160,7 @@ function TrapTriangolo(nxi, nxf, j0, t, y, fy, zn, N, abstol, itmax, s, w, omega
it = it + 1

stop = (norm(yn1-yn0, Inf) < abstol) || (norm(Gn1, Inf)<abstol)
if it > itmax && ~stop
if it > maxiters && ~stop
@warn "Non Convergence"
stop = true
end
Expand All @@ -177,7 +177,7 @@ function TrapTriangolo(nxi, nxf, j0, t, y, fy, zn, N, abstol, itmax, s, w, omega
return y, fy
end

function TrapFirstApproximations(t, y, fy, abstol, itmax, s, halpha, omega, w, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
function TrapFirstApproximations(t, y, fy, abstol, maxiters, s, halpha, omega, w, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial, p)
m = problem_size
Im = zeros(m, m)+I ; Ims = zeros(m*s, m*s)+I
Y0 = zeros(s*m, 1); F0 = copy(Y0); B0 = copy(Y0)
Expand Down Expand Up @@ -221,7 +221,7 @@ function TrapFirstApproximations(t, y, fy, abstol, itmax, s, halpha, omega, w, p
it = it + 1

stop = (norm(Y1-Y0, Inf) < abstol) || (norm(G1, Inf) < abstol)
if it > itmax && ~stop
if it > maxiters && ~stop
@warn "Non Convergence"
stop = 1
end
Expand Down
Loading

0 comments on commit 7eee4aa

Please sign in to comment.