Updated 2012-05-26 04:36:17 by RLE ▲

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 MathematicsCategory Numerical Analysis