I would like to solve convolution equation numerically by the so-called convolution quadrature (CQ) (see the picture) based on Runge-Kutta CQ method. However, I didn’t arrive to get the expected solution as in the attached picture (TextBook; Integral Equation Methods for Evolutionary PDE, page 139) with matlab code. You can see my attempt with julia where all the steps are true,
using TaylorSeries, LinearAlgebra, Plots
function weights( Ks::Function , N , T )
A = [(88 - 7*√6)/360 (296-169* √6)/1800 (-2+3* √6) / 225
(296+169* √6) /1800 (88+7*√6) /360 (-2-3* √6) / 225
(16-√6 )/36 (16+√6)/36 1/9 ] ;
b = [(16-√6)/36 , (16+√6)/36 , 1/9];
c = [(4-√6)/10 , (4 + √6)/10 , 1] ;
Δt = T/N;
m = length(c) ;
ω=Matrix{Vector{Float64}}(undef,m,m)
z = Taylor1(Float64, N+1)
Δ(z) = inv( A + ( z / (1-z) ) * ones(m) * b' );
y = Ks.(Δ(z) / Δt);
for i in 1:3
for j in 1:3
ω[i,j] = y[i,j].coeffs
end
end
return ω, A,b,c;
end
function evalRKCQ( g::Function , Ks::Function , N , T)
Δt = T/N ; ts = (0:N) * Δt ;
ω , A , b , c = weights( Ks , N , T )
m = length(c)
gs = zeros(m , N+1) ;
for j in 1:m
gs[j,:] = g.(ts .+ Δt * c[j]);
end
U=zeros(m,N+1)
for i in 1:m
for j in 1:m
for k in 1:N+1
U[i,k] = sum( ω[i,j][k+1-n] * gs[j,n] for n in 1:k)
end
end
end
return U
end
N = 100; T = 4; ts = (0:N)*T/N;
Ks(s) = 1 /(1 - exp(-s));
g(t) = exp(-100*(t-0.5) ^2);
exact= @. g(ts)+g(ts-1)+g(ts-2)+g(ts-3)
U = evalRKCQ( g , Ks , N , T)
plot(ts[2:N],exact[2:N])
scatter!(ts[2:N], U[3,1:N-1])
This is my plot’s attempt
