# 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 ]