Analysis_LinearVolterraIntegralEquation.jl

include("LinearVolterraIntegralEquation.jl")

using PyPlot 

problem = "1"
Lt,K,g,y_exc = Problems_LVIE(problem)

N_NI = 3

Nm = 3
CollocationParamType = "Gauss"
# CollocationParamType = "Radau"
# CollocationParamType = "Lobatto"

Nt = 100
Lt = 1.0 

t,y = LVIE(Nm,Nt,N_NI,Lt,K,g,CollocationParamType)

Y_exc = y_exc.(t)

figure()
plot(t,Y_exc,"k",lw=3,label="exc")
plot(t,y,"sC3",ms=3,label="num")

getCollocationParam.jl

using FastGaussQuadrature

function getCollocationParam(M,CollocationParamType)::Array{Float64,1}
    if CollocationParamType == "Gauss"
        c,_=gausslegendre(M)
    elseif CollocationParamType == "Lobatto"
        c,_=gausslobatto(M)
    elseif CollocationParamType == "Radau"
        c,_=gaussradau(M)
    end
    # [-1,1] → [0,1]
    c = (c .+ 1)./2
    return c 
end 

LinearVolterraIntegralEquation.jl

include("getCollocationParam.jl")
include("Problems_LVIE.jl")

using FastGaussQuadrature:gausslegendre

function GL1d(f,N,xlm,args)::Float64
    xi, wi = gausslegendre(N)
    xa,xb = xlm

    detJ = (xb-xa)/2.0

    S = 0.0
    for i = 1:N
        x  = (xb-xa)/2.0 * xi[i] + (xb+xa)/2.0
        S += wi[i] * f(x,args) * abs(detJ)
    end 
    return S 
end

function f1(x,args)
    tni,tn,h,j,c,K = args
    return K(tni,tn+x*h)*lj(x,j,c)
end 

function f2(x,args)
    tni,tl,h,j,c,K = args
    return K(tni,tl+x*h)*lj(x,j,c)
end 

function lj(x,j,c)
    m = length(c)

    y = 1.0
    for i = 1:m
        if i != j
            y *= (x - c[i]) / (c[j] - c[i])
        end
    end
    return y
end

function LVIE(Nm,Nt,N_NI,Lt,K,g,CollocationParamType)
    h = Lt/Nt
    t = collect(0:Nt) * h

    c = getCollocationParam(Nm,CollocationParamType)

    Y = zeros(Nt+1,Nm)
    A = zeros(Nt+1,Nm,Nm)

    AA = zeros(Nm,Nm)

    B = zeros(Nt+1,Nt+1,Nm,Nm);

    Im = Array{Float64}(I, Nm, Nm)

    xlm1 = [0.0,1.0]
    G    = zeros(Nm)


    for n = 1:Nt+1
        ########################################################
        for i = 1:Nm
            tni = t[n] + c[i] * h
            for j = 1:Nm
                xlm = [0.0,c[i]]
                args = (tni,t[n],h,j,c,K)
                A[n,i,j] = GL1d(f1,N_NI,xlm,args)

                AA[i,j] = Im[i,j] - h * A[n,i,j]
            end 
        end
        ########################################################
        S1   = zeros(Nm)
        for l = 1:n-1
            for i = 1:Nm 
                tni = t[n] + c[i] * h
                for j = 1:Nm
                    args = (tni,t[l],h,j,c,K)
                    B[n,l,i,j] = GL1d(f2,N_NI,xlm1,args)

                    S1[i] += B[n,l,i,j] * Y[l,j]
                end
            end 
        end
        for i = 1:Nm 
            tni = t[n] + c[i] * h
            G[i] = g(tni)
        end 
        BB = G .+ h * S1
        ########################################################

        Y[n,:] = AA \ BB
    end
    #############################################################
    y = zeros(Nt+1)
    for j = 1:Nm
        y[1] += lj(0.0,j,c) * Y[1,j]
    end

    for n = 1:Nt
        for j = 1:Nm
            y[n+1] += lj(1.0,j,c) * Y[n,j]
        end 
    end

    return t,y
end

LinearVolterraIntegralEquation.jl.55355.mem

        - include("getCollocationParam.jl")
        - include("Problems_LVIE.jl")
        - 
        - using FastGaussQuadrature:gausslegendre
        - 
        - function GL1d(f,N,xlm,args)::Float64
        0     xi, wi = gausslegendre(N)
        0     xa,xb = xlm
        - 
        0     detJ = (xb-xa)/2.0
        - 
        -     S = 0.0
        0     for i = 1:N
        0         x  = (xb-xa)/2.0 * xi[i] + (xb+xa)/2.0
        0         S += wi[i] * f(x,args) * abs(detJ)
        -     end 
        0     return S 
        - end
        - 
        - function f1(x,args)
        0     tni,tn,h,j,c,K = args
        0     return K(tni,tn+x*h)*lj(x,j,c)
        - end 
        - 
        - function f2(x,args)
        0     tni,tl,h,j,c,K = args
        0     return K(tni,tl+x*h)*lj(x,j,c)
        - end 
        - 
        - function lj(x,j,c)
        -     m = length(c)
        - 
        -     y = 1.0
        -     for i = 1:m
        -         if i != j
        -             y *= (x - c[i]) / (c[j] - c[i])
        -         end
        -     end
        -     return y
        - end
        - 
        - function LVIE(Nm,Nt,N_NI,Lt,K,g,CollocationParamType)
     2176     h = Lt/Nt
     1792     t = collect(0:Nt) * h
        - 
        0     c = getCollocationParam(Nm,CollocationParamType)
        - 
     2560     Y = zeros(Nt+1,Nm)
     7424     A = zeros(Nt+1,Nm,Nm)
        - 
      160     AA = zeros(Nm,Nm)
        - 
   734656     B = zeros(Nt+1,Nt+1,Nm,Nm);
        -     # 単位行列
        0     Im = Array{Float64}(I, Nm, Nm)
        -     #####################
       96     xlm1 = [0.0,1.0]
      112     G    = zeros(Nm)
        -     #####################
        - 
        0     for n = 1:Nt+1
        -         ########################################################
        0         for i = 1:Nm
        0             tni = t[n] + c[i] * h
        0             for j = 1:Nm
    87264                 xlm = [0.0,c[i]]
    87264                 args = (tni,t[n],h,j,c,K)
    14544                 A[n,i,j] = GL1d(f1,N_NI,xlm,args)
        - 
        0                 AA[i,j] = Im[i,j] - h * A[n,i,j]
        -             end 
        -         end
        -         ########################################################
    11312         S1   = zeros(Nm)
        0         for l = 1:n-1
        0             for i = 1:Nm 
        0                 tni = t[n] + c[i] * h
        0                 for j = 1:Nm
  4363200                     args = (tni,t[l],h,j,c,K)
   727200                     B[n,l,i,j] = GL1d(f2,N_NI,xlm1,args)
        - 
        0                     S1[i] += B[n,l,i,j] * Y[l,j]
        -                 end
        -             end 
        -         end
        0         for i = 1:Nm 
        0             tni = t[n] + c[i] * h
        0             G[i] = g(tni)
        -         end 
    22624         BB = G .+ h * S1
        -         ########################################################
        - 
        0         Y[n,:] = AA \ BB
        -     end
        -     #############################################################
      896     y = zeros(Nt+1)
        0     for j = 1:Nm
        0         y[1] += lj(0.0,j,c) * Y[1,j]
        -     end
        - 
        0     for n = 1:Nt
        0         for j = 1:Nm
        0             y[n+1] += lj(1.0,j,c) * Y[n,j]
        -         end 
        -     end
        - 
       32     return t,y
        - end
        - 

MemoryProfiling_LVIE1.jl

include("LinearVolterraIntegralEquation.jl")

problem="1"
Lt,K,g,y_exc = Problems_LVIE(problem)
N_NI = 3

Nm = 3
CollocationParamType = "Gauss"
# CollocationParamType = "Radau"
# CollocationParamType = "Lobatto"

Nt = 100

LVIE(Nm,Nt,N_NI,Lt,K,g,CollocationParamType)

using Profile

Profile.clear_malloc_data()
LVIE(Nm,Nt,N_NI,Lt,K,g,CollocationParamType)

Problems_LVIE.jl

function Problems_LVIE(problem)
    ##############################################################
    if problem == "1"
        Lt    = 1.0
        K     = (t,s) -> 2.0 * cos( t - s )
        g     = (t) -> exp(t)
        y_exc = (t) -> exp(t) * ( 1.0 + t )^2 
    ##############################################################
    elseif problem == "2"
        Lt    = 3.0
        K     = (t,s) -> (t-s)
        g     = (t) -> 1.0 - t - t^2/2.0
        y_exc = (t) -> 1.0 - sinh(t)
    end
    return Lt,K,g,y_exc
end

Problems_LVIE.jl.55355.mem

        - function Problems_LVIE(problem)
        -     ##############################################################
        0     if problem == "1"
        -         Lt    = 1.0
        0         K     = (t,s) -> 2.0 * cos( t - s )
        0         g     = (t) -> exp(t)
        0         y_exc = (t) -> exp(t) * ( 1.0 + t )^2 
        -     ##############################################################
        0     elseif problem == "2"
        -         Lt    = 3.0
        0         K     = (t,s) -> (t-s)
        0         g     = (t) -> 1.0 - t - t^2/2.0
        0         y_exc = (t) -> 1.0 - sinh(t)
        -     end
        0     return Lt,K,g,y_exc
        - end