ここで、 y0 はxでの値、y_1 は x+Δxでの値で、ypは途中 で使うだけである。この方法のプログラムを作って、同じ問題を解いてみよう。 誤差はどうなるだろうか。
時間の関係があるので偏微分方程式は扱わない。
program ode1(input, output);
var x, xend, dx, y, k : real;
begin
write('Enter k, y0, xend, dx:');
readln(k, y, xend, dx);
x := 0;
while x <= xend do begin
y := y + dx*k*y;
x := x + dx;
writeln('x, y = ', x:20:10,' ', y:20:10);
end;
end.
上の例は、一階線形常微分方程式の初期値問題
program ode_system(input, output);
const nvar = 2;
type vector = array[1..nvar] of real;
var i : integer;
y : vector;
x : real;
xinit, xend, dx : real;
procedure calculate_dy(x: real; y : vector; var dy : vector);
begin
dy[1] := y[2];
dy[2] := -y[1];
end;
procedure integrate(var y : vector; var x : real; dx : real);
var yp,dy : vector;
i : integer;
begin
calculate_dy(x,y,dy);
for i:=1 to nvar do yp[i] := y[i]+0.5*dx*dy[i];
calculate_dy(x,yp,dy);
for i:=1 to nvar do y[i] := y[i]+dx*dy[i];
x := x + dx;
end;
begin
write('Enter xend, dx: ');
readln(xinit, xend, dx);
y[1]:=0;
y[2]:=1;
x := 0;
while x < xend do begin
integrate(y, x, dx);
write(x:16:7);
for i:= 1 to nvar do write(' ',y[i]:15);
write(y[1]-sin(x), y[2]-cos(x));
writeln;
end;
end.
このプログラムは、（手続き calculate_dy の中身さえ書き換えれば）
どのような連立常微分方程式でも解くことができる（はずである）。
今までに使ったことのない Pascal の言語機能としては
const nvar = 2; type vector = array[1..nvar];の
プログラムでやっていることは、手続き、配列を使っていることの他は上でやっ た1変数常微分方程式とあまり違わない。
procedure integrate(var y : vector;
var x : real;
dx : real);
var y1,dy0,dy1 : vector;
i,k : integer;
begin
calculate_dy(x,y,dy0);
for i:=1 to nvar do y1[i] := y[i]+dx*dy0[i];
for k:=1 to 3 do begin
calculate_dy(x,y1,dy1);
for i:=1 to nvar do y1[i] := y[i]+0.5*dx*(dy0[i] + dy1[i]);
end;
for i:=1 to nvar do y[i] := y1[i];
x := x + dx;
end;
台形公式は、ルンゲ・クッタ公式と同様に高い精度が得られる。この方式の特徴は
y1 = y + 0.5*dx*(dy0 + dy1)としているところで、 dy1 が実は y1 の関数になっているところである。こ のために、ここでは前回やった「直接代入法」で、 y1 についての方程式を解 いている。ただし、 答が正しく求まったかどうかの判定はしないで、最初の 値（オイラー法でもとめた）から数回（この例では3回）繰り返して打ちきっ ている。
台形公式は、計算量が増えるという欠点はあるが、重要な特長がいくつかあり、 この方法を発展させた様々な方法が広く用いられている。
begin
write('Enter xend, dx: ');
readln(xend, dx);
y[1]:=0;
y[2]:=1;
x := 0;
initgraph(gdriver, gmode, '');
moveto(400+round(200*y[1]),400-round(200*y[2]));
while x < xend do begin
integrate(y, x, dx);
lineto(400+round(200*y[1]),400-round(200*y[2]));
end;
readln;
end.
メインプログラムを上のように書き換えると、数値積分した答をグラフィック
表示することができる。（
（省略)
a, b, c,omega : real;
procedure calculate_dy(x : real;y : vector; var dy : vector);
begin
dy[1] := y[2];
dy[2] := -2*b*y[2] -c*y[1] + a* sin(omega*x);
end;
（省略)
begin
write('Enter xinit, xend, dx: ');
readln(xinit, xend, dx);
write('Enter a, b, c, omega ');
readln(a, b,c, omega);
y[1]:=1;
y[2]:=0;
x := xinit;
initgraph(gdriver, gmode, '');
moveto(400+round(200*y[1]),400-round(200*y[2]));
while x < xend do begin
integrate(y, x, dx);
lineto(400+round(200*y[1]),400-round(200*y[2]));
end;
readln;
end.
