Hi everyone,
I need to solve a Volterra integral equation as follows:
Here f(y) and K(x, y) are known functions (nonlinear)
As an example, take f(y) = 3x^3; K(x, y) = 1 - exp(x).
I saw this post about solving Volterra eq of the first kind:
Volterra of the first kind using ApproxFun.jl
But, I do not know what to do for Volterra of the second kind (I can not find documentation on this).
Thanks for checking!
The post you linked is also talking about second kind Volterra. That said, for problems like yours you could use what we have in the examples in GitHub - TSGut/SparseVolterraExamples.jl: A number of examples built on the method described in https://arxiv.org/abs/2005.06081 for solving nonlinear and integro-differential Volterra equations
It isn’t a full fledged package (it is companion software for papers) but it contains all the functionality to adapt it to your problem. Here is what I think is your example running in said package:
julia> using SparseVolterraExamples
julia> K(x, y) = 1 - exp(x)
K (generic function with 1 method)
julia> f(x) = 3x^3
f (generic function with 1 method)
julia> ffun = Fun(x->f(x), Jacobi(0,1,0..1)) # expand in appropriate Jacobi polynomial basis
Fun(Jacobi(0.0,1.0,0..1),[0.3, 0.6, 0.38571428571428573, 0.08571428571428573])
julia> N = 50 # ~approximation order
50
julia> u = triVolterraEQ2FullKernelSolver(ffun,K,N,true)
Fun(Jacobi(0.0,1.0,0..1),[0.25347222053158275, 0.4825526289261275, 0.26627555961222515, 0.020486497311035, -0.018897275106918167, -0.001981957028532612, 0.00027882041579321316, 7.603853364807388e-5, 2.9990342492582807e-6, -1.050172124132875e-6 … 7.394984261821294e-38, 5.200085135033818e-39, 6.765683559556091e-41, -2.820773009293137e-41, -3.60714367988023e-42, -2.119222434012716e-43, 4.548285999212194e-46, 1.4469052747132033e-45, 1.5733346327763277e-46, 7.881027094171445e-48])
plot(u,legend=false)
Does that do what you wanted?
Thank you very much!
I have noticed your package and the paper, going to try it.
Thanks again for the demo!
Hello, I have tried to use your package.
But I run into an error for reproducing your demo, I got:
using SparseVolterraExamples
using IntervalSets
K(x, y) = 1 - exp(x)
f(x) = 3x^3
ffun = Fun(x->f(x), Jacobi(0,1,0..1))
N = 50
u = triVolterraEQ2FullKernelSolver(ffun,K,N,true)
plot(u, legend=false)
Output exceeds the size limit. Open the full output data in a text editor
UndefVarError: _isless not defined
Stacktrace:
[1] getproperty
@ ./Base.jl:31 [inlined]
[2] isless
@ ~/.julia/packages/BlockArrays/8a1Kg/src/blockindices.jl:71 [inlined]
[3] <(x::Block{1, Int64}, y::Block{1, Int64})
@ Base ./operators.jl:356
[4] >(x::Block{1, Int64}, y::Block{1, Int64})
@ Base ./operators.jl:382
[5] BandedBlockBandedMatrix(V::ApproxFunBase.SubOperator{Float64, TimesOperator{Float64, Tuple{InfiniteArrays.Infinity, InfiniteArrays.Infinity}}, Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}, BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}}, Tuple{Int64, Int64}, Tuple{Int64, Int64}})
@ ApproxFunBase ~/.julia/packages/ApproxFunBase/92sMu/src/Operators/general/algebra.jl:487
[6] (AbstractMatrix)(V::ApproxFunBase.SubOperator{Float64, TimesOperator{Float64, Tuple{InfiniteArrays.Infinity, InfiniteArrays.Infinity}}, Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}, BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}}, Tuple{Int64, Int64}, Tuple{Int64, Int64}})
@ ApproxFunBase ~/.julia/packages/ApproxFunBase/92sMu/src/Operators/Operator.jl:781
[7] defaultgetindex(B::TimesOperator{Float64, Tuple{InfiniteArrays.Infinity, InfiniteArrays.Infinity}}, k::BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}, j::BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}})
@ ApproxFunBase ~/.julia/packages/ApproxFunBase/92sMu/src/Operators/Operator.jl:234
[8] getindex(B::TimesOperator{Float64, Tuple{InfiniteArrays.Infinity, InfiniteArrays.Infinity}}, k::BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}, j::BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}})
@ ApproxFunBase ~/.julia/packages/ApproxFunBase/92sMu/src/Operators/Operator.jl:208
[9] plan_evaluate(::Fun{JacobiTriangle, Float64, Vector{Float64}})
@ MultivariateOrthogonalPolynomials ~/.julia/packages/MultivariateOrthogonalPolynomials/kZS4T/src/Triangle/Triangle.jl:422
[10] evaluate(f::Vector{Float64}, K::JacobiTriangle, x::StaticArrays.SVector{2, Float64})
@ MultivariateOrthogonalPolynomials ~/.julia/packages/MultivariateOrthogonalPolynomials/kZS4T/src/Triangle/Triangle.jl:433
[11] evaluate
@ ~/.julia/packages/ApproxFunBase/92sMu/src/Fun.jl:216 [inlined]
...
@ SparseVolterraExamples ~/SparseVolterraExamples/src/triVolterraFullKernel.jl:86
[22] triVolterraEQ2FullKernelSolver(g::Fun{Jacobi{ClosedInterval{Int64}, Float64}, Float64, Vector{Float64}}, Kfun::Function, depth::Int64, flip::Bool)
@ SparseVolterraExamples ~/SparseVolterraExamples/src/triVolterraFullKernel.jl:36
[23] top-level scope
Try using ApproxFun before running it and let me know if that fixes your issue.
No, It does not work. Got the same errors.
It is strange, I have cloned your project and use
Pkg.instantiate() to use your settings.
Sounds like there might be a version or dependency introduced bug. What I ran was on Julia 1.5.3 by just instantiating as you say. Can you tell me your Julia version and installed package versions so I can try to reproduce your bug (feel free to open an issue on the repository for this)?
I am in Julia 1.8.5
I will have a look and see if I can reproduce your error in that version.
Thank you very much!
I have managed to rerun in Julia 1.5.3, it produced the same results as yours.
But what I want is to be able to run your package in the latest Julia version, since I am using the
obtained function results for further calculation.
Yes, I was able to reproduce your bug. I will see what I can do about patching it for 1.8.5.
Thank you very much!
I have tested a relatively complex function which is close to what I need:
using ApproxFun, SparseVolterraExamples, Plots, Distributions
# K(x, y) = 1 - exp(x)
# f(x) = 3x^3
A, B = 6, 2
α, β = 1, 1
m(x) = 1 + A - (A/B) * x
f(x) = pdf(Gamma(α*m(0), β), x)
K(x, y) = pdf(Gamma(α*m(x), β), y-x)
ffun = Fun(x->f(x), Jacobi(0,1,0..1))
N = 50
u = triVolterraEQ2FullKernelSolver(ffun,K,N,true)
plot(u, legend=false)
But I got error says
ERROR: BlockBoundsError: attempt to access 1-blocked 1-element BlockArrays.BlockedUnitRange{Array{Int64,1}} at block index [0]
Stacktrace:
[1] getindex at /home/zuoj/.julia/packages/BlockArrays/8a1Kg/src/blockaxis.jl:161 [inlined]
[2] unblock at /home/zuoj/.julia/packages/BlockArrays/8a1Kg/src/views.jl:10 [inlined]
[3] to_indices at /home/zuoj/.julia/packages/BlockArrays/8a1Kg/src/views.jl:22 [inlined]
[4] to_indices at /home/zuoj/.julia/packages/BlockArrays/8a1Kg/src/views.jl:44 [inlined]
[5] view(::BlockArrays.PseudoBlockArray{Float64,1,Array{Float64,1},Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}}}}, ::BlockArrays.Block{1,Int64}) at ./subarray.jl:157
[6] triVolterraFullKernelOpP01(::typeof(K), ::Int64, ::Bool) at /home/zuoj/SparseVolterraExamples/src/triVolterraFullKernel.jl:96
[7] triVolterraEQ2FullKernelSolver(::Fun{Jacobi{Interval{:closed,:closed,Int64},Float64},Float64,Array{Float64,1}}, ::Function, ::Int64, ::Bool) at /home/zuoj/SparseVolterraExamples/src/triVolterraFullKernel.jl:36
[8] top-level scope at REPL[55]:1
Okay I managed to get it working in 1.8.5, I will push an update in a few minutes and tag a new version for that. Then I will check what is going wrong in the example you posted.
The previous example now works like this
using SparseVolterraExamples
K(x, y) = 1 - exp(x)
f(x) = 3x^3
ffun = Fun(x->f(x), Jacobi(0,1,0..1))
N = 50
M = 7 # total degree 7 approximation in triangle polynomials
u = triVolterraEQ2FullKernelSolver(ffun,K,N,true,M)
plot(u, legend=false)
Note in order to get to a working option for you faster I opted for manual kernel approximation order rather than trying to make it adaptive - if you absolutely need it to be adaptive we can make that work but it would take longer since a lot changed in Julia since that code was written.
Now for the example you gave. It works but there is a strange bug (somewhere in MultivariateOrthogonalPolynomials or FastTransforms I think) which causes functions with globally defined variables to fail when being expanded on the triangle. I will look into this but for now, this already works for your example:
using Distributions, ApproxFun, SparseVolterraExamples
f(x) = pdf(Gamma(1*m(0), 1), x)
K(x, y) = pdf(Gamma(1*(1 + 6 - (6/2) * x), 1), y-x)
ffun = Fun(x->f(x), Jacobi(0,1,0..1))
N = 50
M = 7
u = triVolterraEQ2FullKernelSolver(ffun,K,N,true,M)
plot(u, legend=false)
Thank you very much for making these possible!
I am going to update the package in the latest Julia env.
For the kernel approximation change, I will check and learn more about this.
Thanks for providing the good demo, I will also try to see possible issues.
Glad to help! Let me know (either here, on an issue or also via email) if you need something more elaborate that goes beyond this implementation. I am always interested in Volterra integral applications.
Note that I have started an issue about the bug you had with the pdf example here. It may take some extra time but once that is fixed, it will work with the syntax you used as well.
Hello, I have a question about the u function results.
Here is the problem program:
using SparseVolterraExamples, Plots, ApproxFun
using Distributions
# function m(x)
# n(x) = A - (A/B) * x
# return 1 + max(0, n(x))
# end
# A = 6; B = 10
# 0<=x<4
# y>x
f(x) = pdf(Gamma(1*7, 1), x)
K(x, y) = pdf(Gamma(1*(1 + max(6 - (10/6) * x, 0.0)), 1), y-x)
ffun = Fun(x->f(x), Jacobi(0,1,0..1))
N = 50
M = 4
u = triVolterraEQ2FullKernelSolver(ffun,K,N,true,M)
plot(u, legend=false)
This is one parameter setting for the problem I need to solve. In the problem, the range for x is [0, 4[
I post the literature results as follows:
Yes, the range for x in the spectral method is always 0 to 1 since the resolving basis is Jacobi polynomials shifted to (0,1). But you should be able to simply re-scale the Volterra equation from any finite range to (0,1) prior to plugging it into the solver and then scale the result back afterwards.
