Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/auxiliary.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function FastConv(x, y)
function fast_conv(x, y)
Lx = length(x)
#Ly = size(y, 2)
problem_size = size(y, 1)
Expand Down
23 changes: 13 additions & 10 deletions src/fodesystem/PIPECE.jl
Original file line number Diff line number Diff line change
Expand Up @@ -135,7 +135,7 @@ function solve(prob::FODESystem, h, ::PECE)
t = t0 .+ collect(0:N)*h
y[:, 1] = u0[:, 1]
fy[:, 1] = f_temp
(y, fy) = ABMTriangolo(1, r-1, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, p, f) ;
(y, fy) = ABM_triangolo(1, r-1, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, p, f) ;

# Main process of computation by means of the FFT algorithm
ff = zeros(1, 2^(Qr+2)); ff[1:2] = [0; 2] ; card_ff = 2
Expand Down Expand Up @@ -177,9 +177,9 @@ function DisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn_pred, zn_corr, N, M

stop = (nxi+r-1 == nx0+L*r-1) || (nxi+r-1>=Nr-1)

(zn_pred, zn_corr) = ABMQuadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length)
(zn_pred, zn_corr) = ABM_quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length)

(y, fy) = ABMTriangolo(nxi, nxi+r-1, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, p, f)
(y, fy) = ABM_triangolo(nxi, nxi+r-1, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, p, f)
i_triangolo = i_triangolo + 1

if stop == false
Expand All @@ -200,10 +200,11 @@ function DisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn_pred, zn_corr, N, M
return y, fy
end

function ABMQuadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length)
function ABM_quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length)
coef_end::Int = nxf-nyi+1
i_fft::Int = log2(coef_end/METH.r)
funz_beg::Int = nyi+1; funz_end::Int = nyf+1
funz_beg::Int = nyi+1
funz_end::Int = nyf+1
Nnxf::Int = min(N, nxf)

# Evaluation convolution segment for the predictor
Expand Down Expand Up @@ -240,7 +241,7 @@ end



function ABMTriangolo(nxi, nxf, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, p, f)
function ABM_triangolo(nxi, nxf, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, p, f)

for n = nxi:min(N, nxf)

Expand All @@ -254,7 +255,7 @@ function ABMTriangolo(nxi, nxf, t, y, fy, zn_pred, zn_corr, N, METH, problem_siz
for j = j_beg:n-1
Phi = Phi + METH.bn[1:alpha_length,n-j].*fy[:, j+1]
end
St = StartingTerm(t[n+1], u0, m_alpha, t0, m_alpha_factorial)
St = starting_term(t[n+1], u0, m_alpha, t0, m_alpha_factorial)
y_pred = St + METH.halpha1.*(zn_pred[:, n+1] + Phi)
f_pred = zeros(length(y_pred))
#f_pred = sysf_vectorfield(t[n+1], y_pred, f)
Expand All @@ -271,8 +272,10 @@ function ABMTriangolo(nxi, nxf, t, y, fy, zn_pred, zn_corr, N, METH, problem_siz
Phi = Phi + METH.an[1:alpha_length, n-j+1].*fy[:, j+1]
end
Phi_n = St + METH.halpha2.*(METH.a0[1:alpha_length, n+1].*fy[:, 1] + zn_corr[:, n+1] + Phi)
yn0 = y_pred; fn0 = f_pred
stop = false; mu_it = 0
yn0 = y_pred
fn0 = f_pred
stop = false
mu_it = 0
while stop == false
global yn1 = Phi_n + METH.halpha2.*fn0
mu_it = mu_it + 1
Expand All @@ -299,7 +302,7 @@ end

sysf_vectorfield(t, y, f_fun) = f_fun(t, y)

function StartingTerm(t, u0, m_alpha, t0, m_alpha_factorial)
function starting_term(t, u0, m_alpha, t0, m_alpha_factorial)
ys = zeros(size(u0, 1), 1)
for k = 1 : maximum(m_alpha)
if length(m_alpha) == 1
Expand Down
36 changes: 18 additions & 18 deletions src/fodesystem/bdf.jl
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ function solve(prob::FODESystem, h, ::FLMMBDF; reltol=1e-6, abstol=1e-6)

problem_size = size(u0, 1)

# Number of points in which to evaluate the solution or the BDFWeights
# Number of points in which to evaluate the solution or the BDF_weights
r::Int = 16
N::Int = ceil(Int, (tfinal-t0)/h)
Nr::Int = ceil(Int, (N+1)/r)*r
Expand All @@ -50,8 +50,8 @@ function solve(prob::FODESystem, h, ::FLMMBDF; reltol=1e-6, abstol=1e-6)
fy = zeros(problem_size, N+1)
zn = zeros(problem_size, NNr+1)

# Evaluation of convolution and starting BDFWeights of the FLMM
(omega, w, s) = BDFWeights(alpha, NNr+1)
# Evaluation of convolution and starting BDF_weights of the FLMM
(omega, w, s) = BDF_weights(alpha, NNr+1)
halpha = h^alpha

# Initializing solution and proces of computation
Expand All @@ -61,16 +61,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) = BDFFirstApproximations(t, y, fy, abstol, itmax, s, halpha, omega, w, problem_size, f, Jfdefun, u0, m_alpha, t0, m_alpha_factorial, p)
(y, fy) = BDFTriangolo(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, 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)

# 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) = BDFDisegnaBlocchi(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, itmax, 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 +84,7 @@ function solve(prob::FODESystem, h, ::FLMMBDF; reltol=1e-6, abstol=1e-6)
end


function BDFDisegnaBlocchi(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, itmax, 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 @@ -96,8 +96,8 @@ function BDFDisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N , abstol, itm
i_triangolo = 0; stop = false
while ~stop
stop = (nxi+r-1 == nx0+L*r-1) || (nxi+r-1>=Nr-1)
zn = BDFQuadrato(nxi, nxf, nyi, nyf, fy, zn, omega, problem_size)
(y, fy) = BDFTriangolo(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)
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)
i_triangolo = i_triangolo + 1

if ~stop
Expand All @@ -118,20 +118,20 @@ function BDFDisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N , abstol, itm
return y, fy
end

function BDFQuadrato(nxi, nxf, nyi, nyf, fy, zn, omega, problem_size)
function BDF_quadrato(nxi, nxf, nyi, nyf, fy, zn, omega, problem_size)
coef_beg::Int = nxi-nyf; coef_end::Int = nxf-nyi+1
funz_beg::Int = nyi+1; funz_end::Int = nyf+1
vett_coef = omega[coef_beg+1:coef_end+1]
vett_funz = [fy[:, funz_beg:funz_end] zeros(problem_size, funz_end-funz_beg+1)]
zzn = real(FastConv(vett_coef, vett_funz))
zzn = real(fast_conv(vett_coef, vett_funz))
zn[:, nxi+1:nxf+1] = zn[:, nxi+1:nxf+1] + zzn[:, nxf-nyf:end-1]
return zn
end

function BDFTriangolo(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, itmax, 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 = ABMStartingTerm(t[n1], y0, m_alpha, t0, m_alpha_factorial)
St = ABM_starting_term(t[n1], y0, m_alpha, t0, m_alpha_factorial)

Phi = zeros(problem_size, 1)
for j = 0:s
Expand Down Expand Up @@ -175,7 +175,7 @@ function BDFTriangolo(nxi, nxf, j0, t, y, fy, zn, N, abstol, itmax, s, w, omega,
return y, fy
end

function BDFFirstApproximations(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, itmax, 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 All @@ -184,7 +184,7 @@ function BDFFirstApproximations(t, y, fy, abstol, itmax, s, halpha, omega, w, pr
temp = zeros(length(y[:, 1]))
fdefun(temp, y[:, 1], p, t[j+1])
F0[(j-1)*m+1:j*m, 1] = temp#BDFf_vectorfield(t[j+1], y[:, 1], fdefun)
St = ABMStartingTerm(t[j+1], y0, m_alpha, t0, m_alpha_factorial)
St = ABM_starting_term(t[j+1], y0, m_alpha, t0, m_alpha_factorial)
B0[(j-1)*m+1:j*m, 1] = St + halpha*(omega[j+1]+w[1, j+1])*fy[:, 1]
end
W = zeros(s, s)
Expand Down Expand Up @@ -239,7 +239,7 @@ function BDFFirstApproximations(t, y, fy, abstol, itmax, s, halpha, omega, w, pr
return y, fy
end

function BDFWeights(alpha, N)
function BDF_weights(alpha, N)
# BDF-2 with generating function (2/3/(1-4x/3+x^2/3))^alpha
omega = zeros(1, N+1)
onethird = 1/3; fourthird = 4/3
Expand Down Expand Up @@ -276,7 +276,7 @@ function BDFWeights(alpha, N)
end
end

temp = FastConv([omega zeros(size(omega))], [jj_nu zeros(size(jj_nu))])
temp = fast_conv([omega zeros(size(omega))], [jj_nu zeros(size(jj_nu))])
temp = real.(temp)
b = nn_nu_alpha - temp[:, 1:N+1]
# Solution of the linear system with multiple right-hand side
Expand All @@ -286,7 +286,7 @@ end

Jf_vectorfield(t, y, Jfdefun) = Jfdefun(t, y)

function ABMStartingTerm(t,y0, m_alpha, t0, m_alpha_factorial)
function ABM_starting_term(t,y0, m_alpha, t0, m_alpha_factorial)
ys = zeros(size(y0, 1), 1)
for k = 1:Int64(m_alpha)
ys = ys + (t-t0)^(k-1)/m_alpha_factorial[k]*y0[:, k]
Expand Down
20 changes: 10 additions & 10 deletions src/fodesystem/explicit_pi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ function solve(prob::FODESystem, h, ::PIEX)

# Evaluation of coefficients of the PECE method
nvett = 0:NNr+1 |> collect
bn = zeros(alpha_length, NNr+1); an = copy(bn); a0 = copy(bn)
bn = zeros(alpha_length, NNr+1)#; an = copy(bn); a0 = copy(bn)
for i_alpha = 1:alpha_length
find_alpha = Float64[]
if α[i_alpha] == α[1:i_alpha-1]
Expand Down Expand Up @@ -113,14 +113,14 @@ function solve(prob::FODESystem, h, ::PIEX)
t = t0 .+ collect(0:N)*h
y[:, 1] = u0[:, 1]
fy[:, 1] = f_temp
(y, fy) = PIEXTriangolo(1, r-1, t, y, fy, zn, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f, α, p)
(y, fy) = PIEX_system_triangolo(1, r-1, t, y, fy, zn, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f, α, p)

# Main process of computation by means of the FFT algorithm
ff = zeros(1, 2^(Qr+2)); ff[1:2] = [0; 2] ; card_ff = 2
nx0::Int = 0; ny0::Int = 0
for qr = 0 : Qr
L = 2^qr
(y, fy) = PIDisegnaBlocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, zn, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f, α, p)
(y, fy) = PIEX_system_disegna_blocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, zn, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f, α, p)
ff[1:2*card_ff] = [ff[1:card_ff]; ff[1:card_ff]]
card_ff = 2*card_ff
ff[card_ff] = 4*L
Expand All @@ -138,7 +138,7 @@ function solve(prob::FODESystem, h, ::PIEX)
end


function PIDisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f, alpha, p)
function PIEX_system_disegna_blocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f, alpha, p)

nxi::Int = nx0; nxf::Int = nx0 + L*r - 1
nyi::Int = ny0; nyf::Int = ny0 + L*r - 1
Expand All @@ -154,9 +154,9 @@ function PIDisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N, METH, problem

stop = (nxi+r-1 == nx0+L*r-1) || (nxi+r-1>=Nr-1)

zn = PIEXQuadrato(nxi, nxf, nyi, nyf, fy, zn, N, METH, problem_size, alpha_length, alpha)
zn = PIEX_system_quadrato(nxi, nxf, nyi, nyf, fy, zn, N, METH, problem_size, alpha_length, alpha)

(y, fy) = PIEXTriangolo(nxi, nxi+r-1, t, y, fy, zn, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f, alpha, p)
(y, fy) = PIEX_system_triangolo(nxi, nxi+r-1, t, y, fy, zn, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f, alpha, p)
i_triangolo = i_triangolo + 1

if stop == false
Expand All @@ -177,7 +177,7 @@ function PIDisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N, METH, problem
return y, fy
end

function PIEXQuadrato(nxi, nxf, nyi, nyf, fy, zn, N, METH, problem_size, alpha_length, alpha)
function PIEX_system_quadrato(nxi, nxf, nyi, nyf, fy, zn, N, METH, problem_size, alpha_length, alpha)
coef_end::Int = nxf-nyi+1
i_fft::Int = log2(coef_end/METH.r)
funz_beg::Int = nyi+1; funz_end::Int = nyf+1
Expand All @@ -201,10 +201,10 @@ end



function PIEXTriangolo(nxi, nxf, t, y, fy, zn, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f, alpha, p)
function PIEX_system_triangolo(nxi, nxf, t, y, fy, zn, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f, alpha, p)

for n = nxi:min(N, nxf)
St = StartingTerm(t[n+1], u0, m_alpha, t0, m_alpha_factorial)
St = PIEX_system_starting_term(t[n+1], u0, m_alpha, t0, m_alpha_factorial)
# Evaluation of the predictor
Phi = zeros(problem_size, 1)
if nxi == 1 # Case of the first triangle
Expand Down Expand Up @@ -232,7 +232,7 @@ end

sysf_vectorfield(t, y, f_fun) = f_fun(t, y)

function StartingTerm(t, u0, m_alpha, t0, m_alpha_factorial)
function PIEX_system_starting_term(t, u0, m_alpha, t0, m_alpha_factorial)
ys = zeros(size(u0, 1), 1)
for k = 1 : maximum(m_alpha)
if length(m_alpha) == 1
Expand Down
Loading