*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:

- write the equation as two first-order equations
- supply initial values at one of the sides of the bar for both the temperature itself and its gradient.

- define the equations that are to solved
- set the numerical options, such as the fault tolerance
- create a new command that solves the equations
- apply it to the initial values over a given "period"

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