Information Technology Reference
In-Depth Information
% main program for solution of an ODE:
tstop = 2.0;
dt = 0.004;
n = tstop/dt;
u0 = 1.0;
u = heun(@f1, u0, dt, n);
m = length(u);
fprintf('end value=%g error=%g\n', u(m), u(m)-exp(-dt
*
n));
fid = fopen('sol.dat', 'w');
dump(fid, u, 0.0, dt);
fclose(fid);
Matlab has a wide collection of visualization tools. In particular, it is straightfor-
ward to plot Matlab arrays without going through a data file. The commands listed
in Sect.
6.4.3
can be issued right after
u
is computed:
% visualization:
% plot t and u array, i.e., (t(i),u(i)) data points:
t = [0:dt:tstop];
plot(t, u);
title('Heun method for du/dt = -u');
legend('u(t)');
% attach curve label
6.4.7
Python
The Python code for solving the present ODE problem is quite similar to the Matlab
code. In general, these two languages usually lead to the same nature of compact,
easy-to-understand programs.
Algorithm
6.12
can be implemented directly as a Python function called
heun
:
def heun (f, u0, dt, n):
"""numerical integration of ODEs by Heun's method"""
u = zeros(n+1) # solution array
u[0] = u0 # initial condition
# advance n steps:
for i in range(0,n,1):
v = f(u[i])
u[i+1] = u[i] + 0.5
*
dt
*
(v + f(u[i] + dt
*
v));
return u
In Python, we create arrays when we need them. Arrays are efficiently transferred
to and from functions using references, and Python automatically deletes arrays that
are no longer in use (these two features are shared with Matlab and Java). Python
arrays have zero as the base index, as in C
CC
. The programmer can, however,
create her own arrays with arbitrary base indices. We emphasize that the arrays we
use in Python for numerical computing are created by the
Numeric
module; the
built-in Python lists (also sometimes called arrays) are not well suited for scientific
computations.
File writing in Python follows this setup:
file = open("sol.dat", "w") # open file for writing
file.write("some text, u[%d]=%g\n" % (i,u[i])
file.close()
