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]