## Integral equations

Arjen Markus (7 october 2004) Inspired by the recent pages on differential equations, I decided to set up a page on their counterpart: integral equations. They pop up in all manner of places, but are much less known than differential equations.

The script below is an amateur's attempt to solve a linear Volterra equation of the second kind (and I mean amateur in two meanings: as someone fond of mathematics and as someone not endowed with thorough knowledge of this particular subject). Surprisingly it gives very reasonable results.

(Volterra equations of the first kind require a different strategy - I will try with the midpoint rule instead)

Note:

Volterra equations are typified by the variable bounds on the integral. Fixed bounds typify equations of the Fredholm type.

``` # volterra.tcl --
#    Solve Volterra type integral equations
#
# The procedure volterra solves the following equation for f:
#             x
#            /
#     f(x) + | K(t,x) f(t) dt = G(t)
#            /
#            0
#
# where K(t,x) and G(t) are given functions. Both must be continuous
# and K(t,x) must be not be negative (otherwise the solution is
# unstable).
#

# namespace integralequations
#    Create a convenient namespace for the integral equations procedures
#
namespace eval ::math::integralequations {
#
# Export the various functions
#
namespace export volterra
}

# volterra --
#    Solve a Volterra integral equation
#
# Arguments:
#    K          Procedure to evaluate the kernel function
#    G          Procedure to evaluate the forcing function
#    xend       End of the interval
#    nsteps     Number of steps for sampling the interval (default 50)
#
# Result:
#    A list consisting of x and f(x) where x = 0, step, 2*step, ... xend
#
# Note:
#    The underlying method is the trapezoid rule for integrals.
#
proc ::math::integralequations::volterra {K G xend {nsteps 50} } {

set xstep [expr {\$xend/\$nsteps}]

set g0     [\$G 0.0]
set k0     [\$K 0.0 0.0]
set result [list]
for { set i 1 } { \$i <= \$nsteps } { incr i } {
set xe       [expr {\$i*\$xstep}]
set integral [expr {0.5*\$g0*\$k0}]

foreach {t1 f} \$result {
set k1       [\$K \$t1 \$xe]
set integral [expr {\$integral+\$k1*\$f}]
}
set integral [expr {\$integral*\$xstep}]
set g1       [\$G \$xe]
set k1       [\$K \$xe \$xe]
set f        [expr {(\$g1-\$integral)/(1.0+0.5*\$k1*\$xstep)}]

lappend result \$xe \$f
}

return [concat 0.0 \$g0 \$result]
}

# test --
#    Some tests:
#    K = constant, G = constant
#
proc K1 {t x} {return 1}
proc G1 {x} {return 1}
puts "Exponential solution expected:"
puts [::math::integralequations::volterra K1 G1 1.0 10]

puts "Linear solution expected:"
proc K2 {t x} {return 0}
proc G2 {x} {expr {\$x}}
puts [::math::integralequations::volterra K2 G2 1.0 10]

puts "Third equation (convolution; exact solution: f(x)=sin(x)):"
proc K3 {t x} {expr {\$x-\$t}}
proc G3 {x} {expr {\$x}}
puts [::math::integralequations::volterra K3 G3 1.0 10]
puts [::math::integralequations::volterra K3 G3 10.0]

puts "Fourth equation (convolution; exact solution: f(x)=cos(x)):"
puts [::math::integralequations::volterra K3 G1 10.0]```

 Category Mathematics Category Numerical Analysis