Programs in Pascal for solving various kinds of problems using Numerical Methods
Note: I don't guarantee anything about any of these programs. They are provided as-is. Also, most of them are hacks rather than the most elegant solutions. Sorry about the formatting, I hadn't read Linus' coding style back in 1998. Dated: pre-1999.
- For finding roots of equations of the form y=f(x)=0
- For solving simultaneous equations
- Interpolation
- Numerical Differentiation
- Numerical Integration
Bisection Method
program bisection_method;
var
x1,x2,x3:real;
acc:integer;
function f(t:real):real;
begin
f:=t*t*t-t-1;
{the required function goes here. in this case it's
f(x)= x^3 - x - 1 }
end;
begin
writeln('Bisection method');
writeln('enter the 1st start value: ');
readln(x1);
writeln('enter the 2nd start value: ');
readln(x2);
if (f(x1)*f(x2)>0) or (f(x1)=f(x2)) then
writeln('Invalid starting points')
else
begin
writeln('to how many decimal places? ');
readln(acc);
if f(x1)=0 then begin x2:=x1;x3:=x1;end;
if f(x2)=0 then begin x1:=x2;x3:=x2;end;
while abs(f(x1)-f(x2))>1/exp(acc*(ln(10))) do
begin
x3:=(x1+x2)/2;
if f(x3)*f(x2)>=0 then x2:=x3;
if f(x3)*f(x1)>=0 then x1:=x3;
end;
writeln('The answer is: ',x3:5:acc+1);
end;
readln;
end.
Regula Falsi (False Postion) Method
program false_position;
var
x1,x2,x3:real;
acc:integer;
function f(t:real):real;
begin
f:=t*t*t-2*t-5;
{the required function goes here. in this case it's
f(x)= x^3 - 2x - 5 }
end;
begin
writeln('False Position method');
writeln('enter the 1st start value: ');
readln(x1);
writeln('enter the 2nd start value: ');
readln(x2);
if (f(x1)*f(x2)>0) or (f(x1)=f(x2)) then
writeln('Invalid starting points')
else
begin
writeln('to how many decimal places? ');
readln(acc);
if f(x1)=0 then begin x2:=x1;x3:=x1;end;
if f(x2)=0 then begin x1:=x2;x3:=x2;end;
while (abs(f(x1))>1/exp(acc*(ln(10)))) and (abs(f(x2))>1/exp(acc*(ln(10)))) do
begin
x3:=x1-f(x1)*(x2-x1)/(f(x2)-f(x1));
if f(x3)*f(x2)>=0 then x2:=x3;
if f(x3)*f(x1)>=0 then x1:=x3;
end;
writeln('The answer is: ',x3:5:acc+1);
end;
readln;
end.
Newton-Raphson Method
program newton_rhapson;
var
x1:real;
acc:integer;
function f(t:real):real;
begin
f:=t*sin(t)+cos(t);
end;
function fd(t:real):real;
begin
fd:=t*cos(t);
end;
begin
writeln('Newton-Rhapson method');
writeln('enter the start value: ');
readln(x1);
if fd(x1)=0 then
writeln('Invalid starting points')
else
begin
writeln('to how many decimal places? ');
readln(acc);
while (abs(f(x1))>1/exp(acc*(ln(10)))) do
begin
x1:=x1-f(x1)/fd(x1);
end;
writeln('The answer is: ',x1:5:acc+1);
end;
readln;
end.
Secant Method
program secant;
var
x1,x2,x3:real;
acc:integer;
function f(t:real):real;
begin
f:=t*t*t-2*t*t+t+10;
end;
begin
writeln('Secant method');
writeln('enter the 1st start value: ');
readln(x1);
writeln('enter the 2nd start value: ');
readln(x2);
if f(x1)=f(x2) then
writeln('Invalid starting points')
else
begin
writeln('to how many decimal places? ');
readln(acc);
while (abs(f(x2))>1/exp(acc*(ln(10)))) do
begin
x3:=(x1*f(x2)-x2*f(x1)) / (f(x2)-f(x1));
x1:=x2;x2:=x3;
end;
writeln('The answer is: ',x3:5:acc+1);
end;
readln;
end.
Gauss Elimination Method
program gauss_elimination;
uses crt;
const MAXX=4;
const MAXEQNS=MAXX; {MAXEQNS should be = MAXX}
const SIGPLACES=4;
const DECPLACES=4;
var a:array[1..MAXEQNS,1..MAXX+1] of real;
x:array[1..MAXX] of real;
i,j,k:integer;
temp:real;
begin
clrscr;
{write the equations in general form}
for i:=1 to MAXEQNS do
begin
for j:=1 to MAXX do
begin
write('a',i,j,'.x',j);
if j < MAXX then write(' + ');
end;
writeln(' = b',i);
end;
{get the values of a[i,j]}
for i:=1 to MAXEQNS do
begin
for j:=1 to MAXX do
begin
write('enter a',i,j,' ');
readln(a[i,j]);
end;
write('enter b',i,' ');
readln(a[i,MAXX+1]);
end;
readln;
clrscr;
{write the given eqns}
for i:=1 to MAXEQNS do
begin
for j:=1 to MAXX do
begin
write(a[i,j]:SIGPLACES:DECPLACES,' x',j);
if j < MAXX then write(' + ');
end;
writeln(' = ',a[i,MAXX+1]:SIGPLACES:DECPLACES);
end;
{calculate the output matrix}
for i:=1 to MAXEQNS-1 do
for j:=i+1 to MAXEQNS do
begin
temp:=a[j,i];
for k:=1 to MAXX+1 do
a[j,k]:=a[j,k] - temp*a[i,k]/a[i,i];
end;
{calc the values of x1 x2 etc}
x[MAXEQNS] := a[MAXEQNS,MAXX+1] / a[MAXEQNS,MAXX];
for i:=MAXEQNS downto 1 do
begin
x[i]:=a[i,MAXX+1];
for j:=i+1 to MAXX do
begin
x[i]:=x[i] - a[i,j]*x[j];
end;
x[i]:=x[i]/a[i,i];
end;
{answers}
writeln;
for i:=1 to MAXX do
write('x',i,' = ',x[i]:SIGPLACES:DECPLACES,' ');
readln;
end.
Gauss-Seidel Method
program gauss_seidel;
uses crt;
const MAXX=4;
const MAXEQNS=MAXX; {MAXEQNS should be = MAXX}
const SIGPLACES=4;
var a:array[1..MAXEQNS,1..MAXX+1] of real;
x:array[1..MAXX] of real;
fl:array[1..MAXX] of integer;
max:real;
i,j,acc,flag:integer;
procedure iterate;
var k,m:integer;
xtemp:real;
begin
flag:=0;
for k:=1 to MAXEQNS do
begin
xtemp:=x[k];
x[k]:=a[k,MAXX+1] / a[k,k];
for m:=1 to MAXX do
if k<>m then x[k]:=x[k] - a[k,m]/a[k,k]*x[m];
write(x[k]:SIGPLACES:acc,' ');
if abs(xtemp-x[k]) < 1/(exp(acc*ln(10))) then fl[k]:=1;
end;
writeln;
for k:=1 to MAXX do if fl[k]=1 then flag:=flag+1
end; {proc iterate}
begin {main()}
clrscr;
{write the equations in general form}
for i:=1 to MAXEQNS do
begin
for j:=1 to MAXX do
begin
write('a',i,j,'.x',j);
if j < MAXX then write(' + ');
end;
writeln(' = b',i);
end;
{get the values of a[i,j]}
for i:=1 to MAXEQNS do
begin
max:=-99999999;
for j:=1 to MAXX do
begin
write('enter a',i,j,' ');
readln(a[i,j]);
if a[i,j]>max then max:=a[i,j];
end;
if a[i,i]<>max then
begin
writeln;
writeln('please restart and re-enter the equations ');
writeln('such that the max coefficient in each row ');
writeln('is the a(i,i) value');
readln;
halt;
end;
write('enter b',i,' ');
readln(a[i,MAXX+1]);
end;
readln;
clrscr;
write('enter the number of decimal places to which you want the answer: ');
readln(acc);
writeln;
writeln('the given equations are:');
writeln;
{write the given eqns}
for i:=1 to MAXEQNS do
begin
for j:=1 to MAXX do
begin
write(a[i,j]:SIGPLACES:acc,' x',j);
if j < MAXX then write(' + ');
end;
writeln(' = ',a[i,MAXX+1]:SIGPLACES:acc);
end;
readln;
clrscr;
{seed}
flag:=0;
for i:=1 to MAXX do
begin
x[i]:=0;
fl[i]:=0;
end;
{iterate}
while flag < MAXX do iterate;
readln;
clrscr;
for i:=1 to MAXX do
writeln('x',i,' = ',x[i]);
readln;
end.
Jacobi Method
program jacobi;
uses crt;
const MAXX=4;
const MAXEQNS=MAXX; {MAXEQNS should be = MAXX}
const SIGPLACES=4;
var a:array[1..MAXEQNS,1..MAXX+1] of real;
x:array[1..MAXX] of real;
i,j,acc,flag:integer;
procedure iterate;
var xtemp:array[1..MAXX] of real;
k,m:integer;
begin
for k:=1 to MAXEQNS do
begin
xtemp[k]:=a[k,MAXX+1] / a[k,k];
for m:=1 to MAXX do
if k<>m then xtemp[k]:=xtemp[k] - a[k,m]/a[k,k]*x[m];
write(xtemp[k]:SIGPLACES:acc,' ');
if abs(xtemp[k]-x[k]) < 1/(exp(acc*ln(10))) then flag:=flag+1;
end;
writeln;
for k:=1 to MAXX do x[k]:=xtemp[k];
end; {proc iterate}
begin
clrscr;
{write the equations in general form}
for i:=1 to MAXEQNS do
begin
for j:=1 to MAXX do
begin
write('a',i,j,'.x',j);
if j < MAXX then write(' + ');
end;
writeln(' = b',i);
end;
{get the values of a[i,j]}
for i:=1 to MAXEQNS do
begin
for j:=1 to MAXX do
begin
write('enter a',i,j,' ');
readln(a[i,j]);
end;
write('enter b',i,' ');
readln(a[i,MAXX+1]);
end;
readln;
clrscr;
write('enter the number of decimal places to which you want the answer: ');
readln(acc);
writeln;
writeln('the given equations are:');
writeln;
{write the given eqns}
for i:=1 to MAXEQNS do
begin
for j:=1 to MAXX do
begin
write(a[i,j]:SIGPLACES:acc,' x',j);
if j < MAXX then write(' + ');
end;
writeln(' = ',a[i,MAXX+1]:SIGPLACES:acc);
end;
readln;
clrscr;
{seed}
x[1]:=a[1,MAXX+1] / a[1,1];
for i:=2 to MAXX do x[i]:=0;
{iterate}
flag:=0;
while flag < MAXX do iterate;
readln;
clrscr;
for i:=1 to MAXX do
writeln('x',i,' = ',x[i]:SIGPLACES:acc);
readln;
end.
Newton's Forward Difference Interpolation Method
program forward_diff;
const SIGPLACES=4;
DECPLACES=4;
MAXDATA=100;
MAXDELS=10;
var
n,i,j:integer;
x:array[1..MAXDATA] of real;
y:array[1..MAXDATA,1..MAXDELS] of real;
r,xn,yn,h:real;
function fact(num:longint):longint;
begin
if num>0 then fact:=num*fact(num-1) else fact:=1;
end;
function mult(num:real;subtract:integer):real;
begin
if subtract>0 then mult:=(num-subtract)*mult(num,subtract-1) else mult:=num;
end;
begin
write('how many data points? ');
readln(n);
write('enter starting value of x: ');
readln(x[1]);
write('enter the x interval: ');
readln(h);
for i:=1 to n do
begin
x[i]:=x[1]+(i-1)*h;
write('enter y',i,': ');
readln(y[i,1]);
end;
write('what value of x do you want y for? ');
readln(xn);
for j:=2 to MAXDELS do
begin
for i:=1 to j do y[i,j]:=0;
for i:=j to n do
y[i,j]:=y[i,j-1]-y[i-1,j-1];
end;
for i:=1 to n do
begin
write(x[i]:SIGPLACES:DECPLACES,' ');
for j:=1 to i do
write(y[i,j]:SIGPLACES:DECPLACES,' ');
writeln;
end;
writeln;writeln;
yn:=y[1,1];
r:=(xn-x[1])/h;
for i:=2 to MAXDELS do
yn:=yn+mult(r,i-2)/fact(i-1)*y[i,i];
writeln('answer: ',yn:SIGPLACES:DECPLACES);
readln;
end.
Newton's Backward Difference Interpolation Method
program backward_diff;
uses crt;
const SIGPLACES=4;
DECPLACES=4;
MAXDATA=100;
MAXDELS=10;
var
n,i,j:integer;
x:array[1..MAXDATA] of real;
y:array[1..MAXDATA,1..MAXDELS] of real;
r,xn,yn,h:real;
function fact(num:longint):longint;
begin
if num>0 then fact:=num*fact(num-1) else fact:=1;
end;
function mult(num:real;subtract:integer):real;
begin
if subtract>0 then mult:=(num+subtract)*mult(num,subtract-1) else mult:=num;
end;
begin
clrscr;
write('how many data points? ');
readln(n);
write('enter starting value of x: ');
readln(x[1]);
write('enter the x interval: ');
readln(h);
for i:=1 to n do
begin
x[i]:=x[1]+(i-1)*h;
write('enter y',i,': ');
readln(y[i,1]);
end;
write('what value of x do you want y for? ');
readln(xn);
for j:=2 to MAXDELS do
begin
for i:=1 to j do y[i,j]:=0;
for i:=j to n do
y[i,j]:=y[i,j-1]-y[i-1,j-1];
end;
for i:=1 to n do
begin
write(x[i]:SIGPLACES:DECPLACES,' ');
for j:=1 to i do
write(y[i,j]:SIGPLACES:DECPLACES,' ');
writeln;
end;
writeln;writeln;
yn:=y[n,1];
r:=(xn-x[n])/h;
for i:=2 to MAXDELS do
yn:=yn+mult(r,i-2)/fact(i-1)*y[n,i];
writeln('answer: ',yn:SIGPLACES:DECPLACES);
readln;
end.
Newton's Divided Difference Interpolation Method
program div_diff;
uses crt;
const SIGPLACES=4;
DECPLACES=4;
MAXDATA=100;
MAXDELS=10;
var
n,i,j:integer;
x:array[1..MAXDATA] of real;
y:array[1..MAXDATA,1..MAXDELS] of real;
xn,yn:real;
function mult(k:integer):real;
begin
if k>1 then mult:=mult(k-1)*(xn-x[k]) else mult:=xn-x[1];
end;
begin
clrscr;
write('how many data points? ');
readln(n);
for i:=1 to n do
begin
write('enter x',i,': ');
readln(x[i]);
write('enter y',i,': ');
readln(y[i,1]);
end;
readln;
clrscr;
write('what value of x do you want y for? ');
readln(xn);
{calc the divdiff table}
for j:=2 to MAXDELS do
begin
for i:=1 to j do y[i,j]:=0;
for i:=j to n do
y[i,j]:=(y[i,j-1]-y[i-1,j-1]) / (x[i]-x[i-j+1]);
end;
{show the table}
for i:=1 to n do
begin
write(x[i]:SIGPLACES:DECPLACES,' ');
for j:=1 to i do
write(y[i,j]:SIGPLACES:DECPLACES,' ');
writeln;
end;
writeln;writeln;
yn:=y[1,1];
for i:=2 to MAXDELS do
yn:=yn+mult(i-1)*y[i,i];
writeln('answer: ',yn:SIGPLACES:DECPLACES);
readln;
end.
Euler's Method
program euler; {for diff eqns}
uses crt;
const SIGPLACES=4;
DECPLACES=4;
MAXDATA=100;
var x,y:array[0..MAXDATA] of real;
h,xn:real;
i:integer;
function f(x:real;y:real):real;
begin
f:=x+2*y; {given function}
end;
begin
clrscr;
write('enter x0: '); readln(x[0]);
write('enter y0: '); readln(y[0]);
write('enter the interval h: '); readln(h);
write('enter the value of x for which y is required (xn): ');
readln(xn);
i:=0;
repeat
begin
i:=i+1;
x[i]:=x[0]+i*h;
y[i]:=y[i-1] + (h * f(x[i-1],y[i-1]));
writeln(x[i]:SIGPLACES:DECPLACES,' ',y[i]:SIGPLACES:DECPLACES);
end;
until ((i=MAXDATA) or (x[i]=xn));
if (i=MAXDATA) and (x[i]<>xn) then
writeln('Insufficient storage space, or h is not the correct size to reach xn')
else
writeln('answer is: ',y[i]:SIGPLACES:DECPLACES);
readln;
end.
Euler's Modified Method (Predictor-Corrector formula)
program modified_euler; {for diff eqns}
uses crt;
const SIGPLACES=4;
DECPLACES=4;
MAXDATA=100;
var x,y:array[0..MAXDATA] of real;
h,xn:real;
i:integer;
function f(x:real;y:real):real;
begin
f:=x+3*y; {given function}
end;
begin
clrscr;
write('enter x0: '); readln(x[0]);
write('enter y0: '); readln(y[0]);
write('enter the interval h: '); readln(h);
write('enter the value of x for which y is required (xn): '); readln(xn);
i:=0;
repeat
begin
i:=i+1;
y[i]:=y[i-1] + (h * f(x[i-1],y[i-1]));
x[i]:=x[i-1]+h;
y[i]:=y[i-1] + (h/2) * ( f(x[i-1],y[i-1])+f(x[i],y[i]) );
writeln(xn:SIGPLACES:DECPLACES,' ',y[i]:SIGPLACES:DECPLACES);
end;
until (x[0]+(i*h)=xn);
writeln('answer is: ',y[i]:SIGPLACES:DECPLACES);
readln;
end.
Runge-Kutta 4th Order Method
program runge_kutta_4th;
uses crt;
const SIGPLACES=4;
DECPLACES=4;
MAXDATA=100;
var x,y:array[0..MAXDATA] of real;
i,n:integer;
k1,k2,k3,k4,h,xn:real;
function f(x:real;y:real):real;
begin
f:=-y; { given function y = e^(-x) }
end;
begin
clrscr;
write('enter the first x value (x0): '); readln(x[0]);
write('enter the x value for which y is required (xn): '); readln(xn);
write('enter the step value (h): '); readln(h);
i:=0;
repeat
begin
k1:=h * f(x[i],y[i]);
k2:=h * f(x[i]+h/2,y[i]+k1/2);
k3:=h * f(x[i]+h/2,y[i]+k2/2);
k4:=h * f(x[i]+h,y[i]+k3);
i:=i+1;
x[i]:=x[0]+h*i;
y[i]:=y[i-1] + (k1+2*(k2+k3)+k4)/6;
writeln(x[i]:SIGPLACES:DECPLACES,' ',y[i]:SIGPLACES:DECPLACES);
end;
until ((i=MAXDATA) or (x[i]=xn));
if (i=MAXDATA) and (x[i]<>xn) then
writeln('Insufficient storage space, or h is not the correct size to reach xn')
else
writeln('answer is: ',y[i]:SIGPLACES:DECPLACES);
readln;
end.
Trapezoidal Method
program trapezoidal;
uses crt;
const SIGPLACES=4;
DECPLACES=4;
var x1,xn,h,sum:real;
i:integer;
function f(x:real):real;
begin
f:=1/(1+x);
end;
begin
clrscr;
writeln('Trapezoidal method of numerical integration');
write('what is the starting value of x? ');
readln(x1);
write('what is the ending value of x? ');
readln(xn);
write('enter the interval h: ');
readln(h);
sum:=f(x1)+f(xn);
i:=1;
repeat
begin
i:=i+1;
sum:=sum+2*f(x1+(i-1)*h);
end;
until (x1+i*h = xn);
writeln('answer is ',h*sum/2:SIGPLACES:DECPLACES);
readln;
end.
Simpson's 3/8th Rule Method
program simpson_3_8th;
uses crt;
const SIGPLACES=4;
DECPLACES=4;
MAXDATA=100;
var y:array[0..MAXDATA] of real;
h,sum:real;
i,n:integer;
{i: counter, n:number of ordinates, h:common diff}
begin
clrscr;
write('how many data points? ');
readln(n);
if ((n-1) mod 3)<>0 then {check if numdatapoints=3m+1}
writeln('number of data points must be 3n+1')
else
begin
write('enter the common difference h: ');
readln(h);
for i:=0 to n-1 do
begin
write('enter y',i,': '); readln(y[i]);
end;
sum:=y[0]+y[n-1];
for i:=1 to n-2 do
if i/3=int(i/3) then sum:=sum+2*y[i] else sum:=sum+3*y[i];
writeln('answer is: ',3*h*sum/8:SIGPLACES:DECPLACES);
end;
readln;
end.
