Arjen Markus (9 october 2017) This little program is meant to illustrate the use of the Tclode package, and in particular its use for solving a boundary value problem. Tclode uses the LSODE solver for initial value problems, but can be used for this type of problems as well.
Consider a metal bar through which an electric current is running. The current causes the metal to warm up. However, due to heat exchange with the environment (both conduction and thermal radiation), an equilibrium develops. A simple equation to model this situation is:
d2T
--- - T + 1 = 0
dx2where x is the distance along the bar and T is the temperature.
In this equation we have done away with complication like the actual law for radiation losses (for small temperature differences this can be linearised) and all the coefficients involved.
For our boundary value problem we assume that the temperature at the endpoints of the bar at x = 0 and x = 1 is kept constant at zero. In order to apply Tclode to this problem, we need to:
The complication is that we know the temperature but not the gradient. We have to estimate that in such a way that the temperature at the other end becomes zero (or very near to zero). The estimation method takes most of the code below.
(An alternative is to discretise the differential equation as a set of linear equations, but achieving a high-order approximation is rather involved.)
The way Tclode works is:
The solver command takes the initial values and calculates the new values at the end of the given period. The independent variable t is used for this. At the end the variables associated with the solver are automatically set to the end values. This method of interfacing is used in the code below.
# Use a point-and-shoot method to solve the following
# problem:
# y" - y + 1 = 0
# at t = 0 and t = 1: y = 0
#
# The second-order equation can be written as two first-order equations:
# y' = x
# x' = y - 1
#
# The problem has a straightforward solution, which makes it easy to
# check the numerical solution:
# e 1
# y = 1 - --- exp(-x) - --- exp(x)
# e+1 e+1
#
# The method below varies the initial value of x (=dy/dt) over a
# certain range, searching for a value such that y(1) = 0 or close to 0.
# The exact value for x(0) is (e-1)/(e+1) ~ 0.046212.
#
package require Tclode
#
# Create the solver:
# y is the temperature, x is the first derivative of y
# t is (implicitly) the independent variable
#
set solver \
[tclode::odesolv -atol {1e-6 1e-6} -rtol 1e-4 x {$y - 10.0} y {$x}]
#
# Try to enclose the endpoint we seek: fundamental assumption is that
# the value of y(t=1) is a monotone function of x(t=0).
#
#
set xstart 0.0
set dx 0.25
puts "Start estimate of dx: $dx"
while { $dx < 1000.0 } {
#
# Lower bound for x(t=0)?
#
set t 0.0
set y 0.0
set x [expr {$xstart - $dx}]
$solver run 1.0
set ydown $y
#
# Upper bound for x(t=0)?
#
set t 0.0
set y 0.0
set x [expr {$xstart + $dx}]
$solver run 1.0
set yup $y
#
# Have we got an interval for x(t=0) that includes y = 0
#
puts " Step: $dx -- $ydown -- $yup"
if { $ydown <= 0.0 && $yup >= 0.0 } {
break
}
if { $ydown >= 0.0 && $yup <= 0.0 } {
break
}
#
# Not yet ... increase the interval
#
set dx [expr {2.0 * $dx}]
}
#
# We have a large interval that encloses the solution, so now
# look for a narrower interval
#
puts "Start interval: [expr {$xstart - $dx}] - [expr {$xstart + $dx}]"
set xstart 0.0
while {$dx > 0.001} {
#
# Determine the end value y(t=0) for x(t=0) = xstart
# - we already have the y values for the endpoints
#
set t 0.0
set y 0.0
set x $xstart
$solver run 1.0
set ycentre $y
#
# Now decide which interval to use
#
if { $ycentre < 0.0 && $yup > 0.0 } {
set xstart [expr {$xstart + 0.5 * $dx}]
set ydown $ycentre
}
if { $ycentre > 0.0 && $yup < 0.0 } {
set xstart [expr {$xstart + 0.5 * $dx}]
set ydown $ycentre
}
if { $ycentre < 0.0 && $ydown > 0.0 } {
set xstart [expr {$xstart - 0.5 * $dx}]
set yup $ycentre
}
if { $ycentre > 0.0 && $ydown < 0.0 } {
set xstart [expr {$xstart - 0.5 * $dx}]
set yup $ycentre
}
set dx [expr {0.5 * $dx}]
}
puts "Result: $xstart -- $ycentre -- $yup -- $ydown"
puts ""
#
# Tabulate the solution
#
puts "Table of solution (t, y, dy/dt)"
set t 0.0
set x $xstart
set y 0.0
set tend 0.0
$solver run $tend
while {$tend <= 1.0+0.01} {
$solver continue $tend
puts "$t -- $y -- $x"
set tend [expr {$tend + 0.025}]
}