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]