Numerical Methods.

wxmx 1 / 58

Bisection Methods:
bisection(f, a, b, tol) := block(
[c],
if float(f(a)) · float(f(b)) > 0 then error("f(a) and f(b) must have opposite signs"),
while abs(float(b − a)) > tol do (
c: float((a + b) / 2),
print("a: ", float(a), " b: ", float(b), " c: ", c, " f(c): ", float(f(c))),
if float(f(c)) = 0 then return(c),
if float(f(a)) · float(f(c)) < 0 then b: c else a: c
),
return(float((a + b) / 2))
);

bisection ( f , a , b , tol ) := block ( [ c ] , if float ( f ( a ) ) float ( f ( b ) ) > 0

then error ( f(a) and f(b) must have opposite signs ) , while
a+b
float ( b − a ) > tol do ( c : float ,
2

print ( a: , float ( a ) , b: , float ( b ) , c: , c , f(c): , float ( f ( c ) ) ) , if float ( f ( c ) ) =

0 then return ( c ) , if float ( f ( a ) ) float ( f ( c ) ) < 0 then b : c else a : c ) ,

a+b
return float )
2

1 Find the root of sinx in (0,5)

f(x) := sin(x);
wxplot2d(f(x), [x,−5,5], [y,−3,3]);
root: bisection(f, 0, 5, 1e-4);
print("Root: ", root);
f ( x ) := sin ( x )
Numerical [Link] 2 / 58

a: 0.0 b: 5.0 c: 2.5 f(c):0.5984721441039564

a: 2.5 b: 5.0 c: 3.75 f(c):− 0.5715613187423437
a: 2.5 b: 3.75 c: 3.125 f(c):0.016591892229347906
a: 3.125 b: 3.75 c: 3.4375 f(c):− 0.2916078813138529
a: 3.125 b: 3.4375 c: 3.28125 f(c):− 0.13920380539312127
a: 3.125 b: 3.28125 c: 3.203125 f(c):− 0.06149352449458594
a: 3.125 b: 3.203125 c: 3.1640625 f(c):−
0.022467955642842295
a: 3.125 b: 3.1640625 c: 3.14453125 f(c):−
0.0029385921809077266
a: 3.125 b: 3.14453125 c: 3.134765625 f(c):
0.006826975557195355
a: 3.134765625 b: 3.14453125 c: 3.1396484375 f(c):
0.0019442148649450443
a: 3.1396484375 b: 3.14453125 c: 3.14208984375 f(c):−
−4
4.971901397226883 10

a: 3.1396484375 b: 3.14208984375 c: 3.140869140625 f(c):

−4
7.235129016702296 10
a: 3.140869140625 b: 3.14208984375 c: 3.1414794921875 f(c):
−4
1.1316140205172368 10
a: 3.1414794921875 b: 3.14208984375 c: 3.14178466796875
−4
f(c): − 1.920143777768485 10
Numerical [Link] 3 / 58

a: 3.1414794921875 b: 3.14178466796875 c: 3.141632080078125

−5
f(c): − 3.942648832154713 10
a: 3.1414794921875 b: 3.141632080078125 c:
−5
3.1415557861328125 f(c): 3.68674569723867 10
3.1415939331054688
Root: 3.1415939331054688 3.1415939331054688

2 Find the roots of x^3+2x^2-2=0 in

[0,2]
f(x) := x^3+2·x^2−2;
wxplot2d(f(x), [x,0,2]);
root: bisection(f, 0, 2, 1e-4);
print("Root: ", root);
3 2
f ( x ) := x + 2 x − 2

a: 0.0 b: 2.0 c: 1.0 f(c): 1.0

a: 0.0 b: 1.0 c: 0.5 f(c): − 1.375
a: 0.5 b: 1.0 c: 0.75 f(c): − 0.453125
a: 0.75 b: 1.0 c: 0.875 f(c): 0.201171875
a: 0.75 b: 0.875 c: 0.8125 f(c): − 0.143310546875
a: 0.8125 b: 0.875 c: 0.84375 f(c): 0.024505615234375
a: 0.8125 b: 0.84375 c: 0.828125 f(c): − 0.060497283935546875
Numerical [Link] 4 / 58

a: 0.828125 b: 0.84375 c: 0.8359375 f(c): −

0.01827096939086914
a: 0.8359375 b: 0.84375 c: 0.83984375 f(c):
0.0030483603477478027
a: 0.8359375 b: 0.83984375 c: 0.837890625 f(c): −
0.00762852281332016
a: 0.837890625 b: 0.83984375 c: 0.8388671875 f(c): −
0.002294388599693775
a: 0.8388671875 b: 0.83984375 c: 0.83935546875 f(c):
−4
3.7590868305414915 10
a: 0.8388671875 b: 0.83935546875 c: 0.839111328125 f(c): −
−4
9.595092124072835 10
a: 0.839111328125 b: 0.83935546875 c: 0.8392333984375 f(c):
−4
− 2.91867583655403 10
a: 0.8392333984375 b: 0.83935546875 c: 0.83929443359375
−5
f(c): 4.2003719272543094 10 0.839263916015625
Root: 0.839263916015625 0.839263916015625

3 Find the roots of tan(x)-x in [0,10]

f(x) := tan(x)−x;
wxplot2d(f(x), [x,0,10]);
print("Finding First Root");
root1: bisection(f, 0, 2, 1e-4);
print("First Root: ", float(root1));
print("Finding second Root");
root2: bisection(f, 4.5, 6, 1e-4);
print("Second Root: ", float(root2));
print("Finding third Root");
root3: bisection(f, 7.8, 8, 1e-4);
print("Third Root: ", float(root3));
f ( x ) := tan ( x ) − x
Numerical [Link] 5 / 58

Finding First Root Finding First Root

a: 0.0 b: 2.0 c: 1.0 f(c):0.5574077246549023
a: 1.0 b: 2.0 c: 1.5 f(c):12.601419947171719
a: 1.5 b: 2.0 c: 1.75 f(c):− 7.27037992250933
a: 1.5 b: 1.75 c: 1.625 f(c):− 20.05586276236962
a: 1.5 b: 1.625 c: 1.5625 f(c):118.97000572254261
a: 1.5625 b: 1.625 c: 1.59375 f(c):− 45.15211040673973
a: 1.5625 b: 1.59375 c: 1.578125 f(c):− 138.02602884284482
a: 1.5625 b: 1.578125 c: 1.5703125 f(c):2065.2848772466036
a: 1.5703125 b: 1.578125 c: 1.57421875 f(c):−
293.76371015449143
a: 1.5703125 b: 1.57421875 c: 1.572265625 f(c):−
682.1688095492681
a: 1.5703125 b: 1.572265625 c: 1.5712890625 f(c):−
2031.0566880568701
a: 1.5703125 b: 1.5712890625 c: 1.57080078125 f(c):−
224495.92011736153
a: 1.5703125 b: 1.57080078125 c: 1.570556640625 f(c):
4170.551603271635
a: 1.570556640625 b: 1.57080078125 c: 1.5706787109375 f(c):
8500.684183205163
a: 1.5706787109375 b: 1.57080078125 c: 1.57073974609375
f(c):17672.30000889566 1.570770263671875
First Root:1.570770263671875 1.570770263671875
Numerical [Link] 6 / 58

Finding second Root Finding second Root

a: 4.5 b: 6.0 c: 5.25 f(c):− 6.927326406624389
a: 4.5 b: 5.25 c: 4.875 f(c):− 10.970345265174792
a: 4.5 b: 4.875 c: 4.6875 f(c):35.48262705640037
a: 4.6875 b: 4.875 c: 4.78125 f(c):− 19.280293033974072
a: 4.6875 b: 4.78125 c: 4.734375 f(c):− 50.21049500483272
a: 4.6875 b: 4.734375 c: 4.7109375 f(c):684.2403623472053
a: 4.7109375 b: 4.734375 c: 4.72265625 f(c):−
102.11611120711501
a: 4.7109375 b: 4.72265625 c: 4.716796875 f(c):−
231.581005471694
a: 4.7109375 b: 4.716796875 c: 4.7138671875 f(c):−
681.2085621982591
a: 4.7109375 b: 4.7138671875 c: 4.71240234375 f(c):−
74836.16217057766
a: 4.7109375 b: 4.71240234375 c: 4.711669921875 f(c):
1385.9955036611636
a: 4.711669921875 b: 4.71240234375 c: 4.7120361328125 f(c):
2829.3728132917913
a: 4.7120361328125 b: 4.71240234375 c: 4.71221923828125
f(c):5886.577980015014
a: 4.71221923828125 b: 4.71240234375 c: 4.712310791015625
f(c):12784.750117696492 4.7123565673828125
Second Root: 4.7123565673828125 4.7123565673828125
Finding third Root Finding third Root
a: 7.8 b: 8.0 c: 7.9 f(c):− 29.61511268125451
a: 7.8 b: 7.9 c: 7.85 f(c):243.3018441758436
a: 7.85 b: 7.9 c: 7.875 f(c):− 55.44543136522765
a: 7.85 b: 7.875 c: 7.8625 f(c):− 125.25306637947396
a: 7.85 b: 7.8625 c: 7.856249999999999 f(c):−
448.70145469326593
a: 7.85 b: 7.856249999999999 c: 7.8531249999999995 f(c):
1159.5063134293748
a: 7.8531249999999995 b: 7.856249999999999 c:
7.854687499999999 f(c):− 1424.5538878715226
a: 7.8531249999999995 b: 7.854687499999999 c: 7.85390625
f(c):13257.564970347914
a: 7.85390625 b: 7.854687499999999 c: 7.854296874999999
f(c):− 3180.030143076836
a: 7.85390625 b: 7.854296874999999 c: 7.8541015624999995
f(c):− 8346.153858855943
a: 7.85390625 b: 7.8541015624999995 c: 7.85400390625 f(c):−
Numerical [Link] 7 / 58

44906.723860095175 7.853955078125
Third Root: 7.853955078125 7.853955078125

4 Find the root of f(x) = cos(x) - x*e(x)

in [0,2] by using bisection. Also plot
the graph of the given function.
f(x) := cos(x)−x·exp(x);
wxplot2d(f(x), [x,0,2]);
root: bisection(f, 0, 2, 1e-4);
print("Root: ", root);
f ( x ) := cos ( x ) − x exp ( x )

a: 0.0 b: 2.0 c: 1.0 f(c): − 2.1779795225909053

a: 0.0 b: 1.0 c: 0.5 f(c): 0.05322192654030866
a: 0.5 b: 1.0 c: 0.75 f(c): − 0.8560611435856852
a: 0.5 b: 0.75 c: 0.625 f(c): − 0.356690603889921
a: 0.5 b: 0.625 c: 0.5625 f(c): − 0.1412937453091
a: 0.5 b: 0.5625 c: 0.53125 f(c): − 0.04151221167208241
a: 0.5 b: 0.53125 c: 0.515625 f(c): 0.006475340827341247
a: 0.515625 b: 0.53125 c: 0.5234375 f(c): −
0.01736202527677333
a: 0.515625 b: 0.5234375 c: 0.51953125 f(c): −
0.005404401827669103
a: 0.515625 b: 0.51953125 c: 0.517578125 f(c):
Numerical [Link] 8 / 58

−4
5.451844176196374 10

a: 0.517578125 b: 0.51953125 c: 0.5185546875 f(c): −

0.0024271774567166116
a: 0.517578125 b: 0.5185546875 c: 0.51806640625 f(c): −
−4
9.403890225676559 10
a: 0.517578125 b: 0.51806640625 c: 0.517822265625 f(c): −
−4
1.9745046759500617 10
a: 0.517578125 b: 0.517822265625 c: 0.5177001953125 f(c):
−4
1.7390492881264041 10
a: 0.5177001953125 b: 0.517822265625 c: 0.51776123046875
−5
f(c): − 1.1763280326260173 10 0.517730712890625
Root: 0.517730712890625 0.517730712890625

5 Find the root of f(x)=

x^3+2*x^2-3*x-1 by using bisection
method. Also plot the graph of the
given function.
f(x) := x^3+2·x^2−3·x−1;
wxplot2d(f(x), [x,−5,5]);
root1: bisection(f, −4, −2, 1e-4);
root2: bisection(f, −2, 0, 1e-4);
root3: bisection(f, 0, 2, 1e-4);

print("Root 1: ", float(root1));

print("Root 2: ", float(root2));
print("Root 3: ", float(root3));
3 2
f ( x ) := x + 2 x + − 3 x − 1
Numerical [Link] 9 / 58

a: − 4.0 b: − 2.0 c: − 3.0 f(c): − 1.0

a: − 3.0 b: − 2.0 c: − 2.5 f(c): 3.375
a: − 3.0 b: − 2.5 c: − 2.75 f(c): 1.578125
a: − 3.0 b: − 2.75 c: − 2.875 f(c): 0.392578125
a: − 3.0 b: − 2.875 c: − 2.9375 f(c): − 0.277099609375
a: − 2.9375 b: − 2.875 c: − 2.90625 f(c): 0.064300537109375
a: − 2.9375 b: − 2.90625 c: − 2.921875 f(c): −
0.10474777221679688
a: − 2.921875 b: − 2.90625 c: − 2.9140625 f(c): −
0.01981210708618164
a: − 2.9140625 b: − 2.90625 c: − 2.91015625 f(c):
0.022346913814544678
a: − 2.9140625 b: − 2.91015625 c: − 2.912109375 f(c):
0.0012931004166603088
a: − 2.9140625 b: − 2.912109375 c: − 2.9130859375 f(c): −
0.009253076277673244
a: − 2.9130859375 b: − 2.912109375 c: − 2.91259765625 f(c): −
0.003978381515480578
a: − 2.91259765625 b: − 2.912109375 c: − 2.912353515625 f(c):
− 0.0013422389893094078
a: − 2.912353515625 b: − 2.912109375 c: − 2.9122314453125
−5
f(c): − 2.446890175633598 10
a: − 2.9122314453125 b: − 2.912109375 c: − 2.91217041015625
Numerical [Link] 10 / 58

−4
f(c): 6.343408529119188 10 − 2.912200927734375
a: − 2.0 b: 0.0 c: − 1.0 f(c): 3.0
a: − 1.0 b: 0.0 c: − 0.5 f(c): 0.875
a: − 0.5 b: 0.0 c: − 0.25 f(c): − 0.140625
a: − 0.5 b: − 0.25 c: − 0.375 f(c): 0.353515625
a: − 0.375 b: − 0.25 c: − 0.3125 f(c): 0.102294921875
a: − 0.3125 b: − 0.25 c: − 0.28125 f(c): − 0.020294189453125
a: − 0.3125 b: − 0.28125 c: − 0.296875 f(c):
0.040729522705078125
a: − 0.296875 b: − 0.28125 c: − 0.2890625 f(c):
0.01014852523803711
a: − 0.2890625 b: − 0.28125 c: − 0.28515625 f(c): −
0.005090296268463135
a: − 0.2890625 b: − 0.28515625 c: − 0.287109375 f(c):
0.0025247707962989807
a: − 0.287109375 b: − 0.28515625 c: − 0.2861328125 f(c): −
0.0012838514521718025
a: − 0.287109375 b: − 0.2861328125 c: − 0.28662109375 f(c):
−4
6.201878422871232 10
a: − 0.28662109375 b: − 0.2861328125 c: − 0.286376953125
−4
f(c): − 3.318998060422018 10
a: − 0.28662109375 b: − 0.286376953125 c: − 0.2864990234375
−4
f(c): 1.4412702330446336 10
a: − 0.2864990234375 b: − 0.286376953125 c: −
−5
0.28643798828125 f(c): − 9.389064075548958 10 −
0.286468505859375
a: 0.0 b: 2.0 c: 1.0 f(c): − 1.0
a: 1.0 b: 2.0 c: 1.5 f(c): 2.375
a: 1.0 b: 1.5 c: 1.25 f(c): 0.328125
a: 1.0 b: 1.25 c: 1.125 f(c): − 0.419921875
a: 1.125 b: 1.25 c: 1.1875 f(c): − 0.067626953125
a: 1.1875 b: 1.25 c: 1.21875 f(c): 0.124725341796875
a: 1.1875 b: 1.21875 c: 1.203125 f(c): 0.027179718017578125
a: 1.1875 b: 1.203125 c: 1.1953125 f(c): −
0.020564556121826172
a: 1.1953125 b: 1.203125 c: 1.19921875 f(c):
0.003222167491912842
a: 1.1953125 b: 1.19921875 c: 1.197265625 f(c): −
0.008692525327205658
a: 1.197265625 b: 1.19921875 c: 1.1982421875 f(c): −
Numerical [Link] 11 / 58

0.00274051446467638
a: 1.1982421875 b: 1.19921875 c: 1.19873046875 f(c):
−4
2.3949227761477232 10
a: 1.1982421875 b: 1.19873046875 c: 1.198486328125 f(c): −
0.0012508446088759229
a: 1.198486328125 b: 1.19873046875 c: 1.1986083984375 f(c):
−4
− 5.057595499238232 10
a: 1.1986083984375 b: 1.19873046875 c: 1.19866943359375
−4
f(c): − 1.3315448290995846 10 1.198699951171875
Root 1: − 2.912200927734375 − 2.912200927734375
Root 2: − 0.286468505859375 − 0.286468505859375
Root 3: 1.198699951171875 1.198699951171875

6 Find the root of f(x)= sin(2x+9)=0 by

using bisection method. Also plot
the graph of the given function.

f(x) := sin(2·x+9);
wxplot2d(f(x), [x,−2,2]);
root: bisection(f, 0, 1, 1e-4);
print("Root: ", root);
f ( x ) := sin ( 2 x + 9 )
Numerical [Link] 12 / 58

a: 0.0 b: 1.0 c: 0.5 f(c):− 0.5440211108893698

a: 0.0 b: 0.5 c: 0.25 f(c):− 0.0751511204618093
a: 0.0 b: 0.25 c: 0.125 f(c):0.17388948538043356
a: 0.125 b: 0.25 c: 0.1875 f(c):0.04975740630107862
a: 0.1875 b: 0.25 c: 0.21875 f(c):− 0.012721696055790868
a: 0.1875 b: 0.21875 c: 0.203125 f(c):0.01852690072471944
a: 0.203125 b: 0.21875 c: 0.2109375 f(c):
0.0029029566920853495
a: 0.2109375 b: 0.21875 c: 0.21484375 f(c):−
0.004909519507735837
a: 0.2109375 b: 0.21484375 c: 0.212890625 f(c):−
0.0010032890623035957
a: 0.2109375 b: 0.212890625 c: 0.2119140625 f(c):
−4
9.498356265579848 10

a: 0.2119140625 b: 0.212890625 c: 0.21240234375 f(c): −

−5
2.6726730617102714 10
a: 0.2119140625 b: 0.21240234375 c: 0.212158203125 f(c):
−4
4.615545029920243 10
a: 0.212158203125 b: 0.21240234375 c: 0.2122802734375 f(c):
−4
2.174138926668997 10
a: 0.2122802734375 b: 0.21240234375 c: 0.21234130859375
−5
f(c): 9.534358173526353 10 0.212371826171875
Root: 0.212371826171875 0.212371826171875

7 Find the root of f(x)=x^3+2x^2-3x-1

= 0 with error 0.01 by using bisection
method.
f(x) := x^3+2·x^2−3·x−1;
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: bisection(f, 0, 2, 1e-2);
print("Root: ", root);
3 2
f ( x ) := x + 2 x + − 3 x − 1
plot2d: some values will be clipped.
Numerical [Link] 13 / 58

a: 0.0 b: 2.0 c: 1.0 f(c): − 1.0

a: 1.0 b: 2.0 c: 1.5 f(c): 2.375
a: 1.0 b: 1.5 c: 1.25 f(c): 0.328125
a: 1.0 b: 1.25 c: 1.125 f(c): − 0.419921875
a: 1.125 b: 1.25 c: 1.1875 f(c): − 0.067626953125
a: 1.1875 b: 1.25 c: 1.21875 f(c): 0.124725341796875
a: 1.1875 b: 1.21875 c: 1.203125 f(c): 0.027179718017578125
a: 1.1875 b: 1.203125 c: 1.1953125 f(c): −
0.020564556121826172 1.19921875
Root: 1.19921875 1.19921875

Regula Falsi
regula_falsi(f, a, b, tol) := block(
[c],
if float(f(a)) · float(f(b)) > 0 then error("f(a) and f(b) must have opposite signs"),
while abs(float(b − a)) > tol do (
c: float(a − f(a) · (b − a) / (f(b) − f(a))),
print("a: ", float(a), " b: ", float(b), " c: ", float(c), " f(c): ", float(f(c))),
if float(f(c)) = 0 then return(c),
if float(f(a)) · float(f(c)) < 0 then b: c else a: c
),
return(float(c))
);
Numerical [Link] 14 / 58

regula_falsi
(f,a ,b ,tol):=block([c ],if float(f(a )) float(f(b ))>

0 then error(f(a) and f(b) must have opposite signs

) ,while

float(b − a ) > tol do (c:floata −

f(a) (b−a)
,
f ( b )− f ( a )

print ( a: , float ( a ) , b: , float ( b ) , c: , float ( c ) , f(c): , float ( f ( c ) ) ) , if

float ( f ( c ) ) = 0 then return ( c ) , if float ( f ( a ) ) float ( f ( c ) ) < 0 then b : c

else a : c ) , return ( float ( c ) ) )

1 Find the roots of cosx=0 by

regula-falsi method.
f(x) := cos(x);
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: regula_falsi(f, 0, 2, 1e-4);
print("Root: ", root);
f ( x ) := cos ( x )

a: 0.0 b: 2.0 c: 1.412282927437392 f(c): 0.15785041839594702

a: 1.412282927437392 b: 2.0 c: 1.5739063237228794 f(c): −
0.003109991914628247
a: 1.412282927437392 b: 1.5739063237228794 c:
−5
1.570783521943903 f(c): 1.2804850993242101 10
a: 1.570783521943903 b: 1.5739063237228794 c:
Numerical [Link] 15 / 58

−11
1.5707963268154532 f(c): − 2.0556606247014516 10
1.5707963268154532
Root: 1.5707963268154532 1.5707963268154532

2 Find the root of sin(2x+9)=0 by using

regula-falsi method.
f(x) := sin(2·x+9);
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: regula_falsi(f, 0, 1, 1e-4);
print("Root: ", root);
f ( x ) := sin ( 2 x + 9 )

−4
regula_falsi ( f , 0 , 1 , 1.0 10 )
−4
Root: regula_falsi ( f , 0 , 1 , 1.0 10 )
−4
regula_falsi ( f , 0 , 1 , 1.0 10 )

3 Find the root of sin(2x+9)=0 by using

regula-falsi method with error 0.01
f(x) := sin(2·x+9);
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: regula_falsi(f, 0, 1, 1e-2);
print("Root: ", root);
Numerical [Link] 16 / 58

f(x ):=sin(2 x+ 9 )

regula_falsi
(f,0 ,1 ,0.01)

Root: regula_falsi
(f,0 ,1 ,0.01) regula_falsi
(f,0 ,1 ,0.01)

4 Find the positive root of x^3= 2x +5.

(Do only 4 iterations).
regula_falsi_iter(f, a, b, max_iter) := block(
[c, iter: 0],
if float(f(a)) · float(f(b)) > 0 then error("f(a) and f(b) must have opposite signs"),
while iter < max_iter do (
c: float(a − f(a) · (b − a) / (f(b) − f(a))),
print("Iteration: ", iter + 1, " a: ", float(a), " b: ", float(b), " c: ", c, " f(c): ", float(f
if abs(float(f(c))) < 1e-6 then return(c),
if float(f(a)) · float(f(c)) < 0 then b: c else a: c,
iter: iter + 1
),
return(c)
);

f(x) := x^3−2·x−5;
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: regula_falsi_iter(f, 2, 3, 4);
print("Root: ", root);
Numerical [Link] 17 / 58

regula_falsi_iter
(f,a ,b ,max_iter):=block([c,iter
:0 ],if

float(f(a )) float(f(b ))> 0 then

error(f(a) and f(b) must have opposite signs

) ,while iter
< max_iter do

(c:floata −
f(a) (b−a)
, print ( Iteration: , iter + 1 , a: , float ( a ) , b: ,
f ( b )− f ( a )
−6
float ( b ) , c: , c , f(c): , float ( f ( c ) ) ) , if float ( f ( c ) ) < 1.0 10 then

return ( c ) , if float ( f ( a ) ) float ( f ( c ) ) < 0 then b : c else a : c , iter : iter + 1

3
) , return ( c ) ) f ( x ) := x − 2 x − 5
plot2d: some values will be clipped.

Iteration: 1 a: 2.0 b: 3.0 c: 2.0588235294117645 f(c): −

0.39079991858335283
Iteration: 2 a: 2.0588235294117645 b: 3.0 c:
2.081263659845023 f(c): − 0.14720405955375426
Iteration: 3 a: 2.081263659845023 b: 3.0 c:
2.0896392100908474 f(c): − 0.054676503273289434
Iteration: 4 a: 2.0896392100908474 b: 3.0 c:
2.0927395743180055 f(c): − 0.020202866312461154
2.0927395743180055
Root: 2.0927395743180055 2.0927395743180055
Numerical [Link] 18 / 58

5 Find the root of xe^x=3 using

regula-falsi method.

f(x) := x·exp(x)−3;
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: regula_falsi(f, 1, 2, 1e-4);
print("Root: ", root);
f ( x ) := x exp ( x ) − 3
plot2d: some values will be clipped.

−4
regula_falsi ( f , 1 , 2 , 1.0 10 )
−4
Root: regula_falsi ( f , 1 , 2 , 1.0 10 )
−4
regula_falsi ( f , 1 , 2 , 1.0 10 )

Secant Method
Numerical [Link] 19 / 58

secant(f, x0, x1, tol) := block(

[x2],
while abs(float(x1 − x0)) > tol do (
x2: float(x1 − f(x1) · (x1 − x0) / (f(x1) − f(x0))),
print("x0: ", float(x0), " x1: ", float(x1), " x2: ", x2, " f(x2): ", float(f(x2))),
x0: x1,
x1: x2
),
return(x1)
);
secant(f,x0 ,x1 ,tol):=block([x2 ],while float(x1 − x0 ) > tol

do (x2 :floatx1 −
f ( x1 ) ( x1 − x0 )
, print
f ( x1 ) − f ( x0 )

( x0: , float ( x0 ) , x1: , float ( x1 ) , x2: , x2 , f(x2): , float ( f ( x2 ) ) ) , x0 : x1 , x1 :

x2 ) , return ( x1 ) )

1 Find the root of cosx =0 by secant

method.
f(x) := cos(x);
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: secant(f, 0, 2, 1e-4);
print("Root: ", root);
f ( x ) := cos ( x )
Numerical [Link] 20 / 58

x0: 0.0 x1: 2.0 x2: 1.4122829274373918 f(x2):

0.15785041839594724
x0: 2.0 x1: 1.4122829274373918 x2: 1.5739063237228794 f(x2):
− 0.003109991914628247
x0: 1.4122829274373918 x1: 1.5739063237228794 x2:
−5
1.570783521943903 f(x2):1.2804850993242101 10

x0: 1.5739063237228794 x1: 1.570783521943903 x2:

−11
1.5707963268154532 f(x2): − 2.0556606247014516 10
1.5707963268154532
Root: 1.5707963268154532 1.5707963268154532

2 Find the root of sinx = 0 by secant

method.
f(x) := sin(x);
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: secant(f, 3, 4, 1e-4);
print("Root: ", root);
f ( x ) := sin ( x )
Numerical [Link] 21 / 58

x0: 3.0 x1: 4.0 x2: 3.157162792479947 f(x2):−

0.015569509788328599
x0: 4.0 x1: 3.157162792479947 x2: 3.1394590982180786 f(x2):
0.0021335537530368953
x0: 3.157162792479947 x1: 3.1394590982180786 x2:
−8
3.141592727984857 f(x2):− 7.439506376460243 10

x0: 3.1394590982180786 x1: 3.141592727984857 x2:

−14
3.1415926535897367 f(x2): 5.652179433087269 10
3.1415926535897367
Root: 3.1415926535897367 3.1415926535897367

3 Find the root of sin(2x+9)=0 by using

secant method.
f(x) := sin(2·x+9);
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: secant(f, 0, 1, 1e-4);
print("Root: ", root);
f ( x ) := sin ( 2 x + 9 )
Numerical [Link] 22 / 58

x0: 0.0 x1: 1.0 x2: 0.29184615011372395 f(x2):−

0.1582463192264177
x0: 1.0 x1: 0.29184615011372395 x2: 0.1587145148526203
f(x2):0.10714287201900013
x0: 0.29184615011372395 x1: 0.1587145148526203 x2:
−4
0.2124623924188545 f(x2):− 1.4682406780122068 10

x0: 0.1587145148526203 x1: 0.2124623924188545 x2:

−7
0.21238883939085848 f(x2): 2.8198766203292956 10
0.21238883939085848
Root: 0.21238883939085848 0.21238883939085848

4 Find the root of x^3 = 2*x +5 by

using secant method
f(x) :=x^3 − 2·x −5;
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: secant(f, 2, 3, 1e-4);
print("Root: ", root);
3
f ( x ) := x − 2 x − 5
plot2d: some values will be clipped.
Numerical [Link] 23 / 58

x0: 2.0 x1: 3.0 x2: 2.0588235294117645 f(x2): −

0.39079991858335283
x0: 3.0 x1: 2.0588235294117645 x2: 2.081263659845023 f(x2):
− 0.14720405955375426
x0: 2.0588235294117645 x1: 2.081263659845023 x2:
2.0948241460940524 f(x2): 0.003043795598889787
x0: 2.081263659845023 x1: 2.0948241460940524 x2:
−5
2.0945494310352473 f(x2): − 2.2886580653747046 10
x0: 2.0948241460940524 x1: 2.0945494310352473 x2:
−9
2.094551481227599 f(x2): − 3.512811375117053 10
2.094551481227599
Root: 2.094551481227599 2.094551481227599

5 Find the smallest positive root of

x^4-x-10=0 using secant method
f(x) :=x^4−x−10;
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: secant(f, 1.5, 2, 1e-4);
print("Root: ", root);
4
f ( x ) := x − x − 10
plot2d: some values will be clipped.
Numerical [Link] 24 / 58

x0: 1.5 x1: 2.0 x2: 1.8083832335329342 f(x2): −

1.1138486987006448
x0: 2.0 x1: 1.8083832335329342 x2: 1.8501193306593111 f(x2):
− 0.13359056028927263
x0: 1.8083832335329342 x1: 1.8501193306593111 x2:
1.8558071677513204 f(x2): 0.005468272484490555
x0: 1.8501193306593111 x1: 1.8558071677513204 x2:
−5
1.8555835023972678 f(x2): − 2.52009818684229 10
x0: 1.8558071677513204 x1: 1.8555835023972678 x2:
−9
1.855584528448739 f(x2): − 4.7197445951496775 10
1.855584528448739
Root: 1.855584528448739 1.855584528448739

Newton Raphson Method

Numerical [Link] 25 / 58

newton_raphson(f, x0, tol) := block(

[x1, f_prime],
f_prime: diff(f(x), x),
x1: float(x0 − f(x0) / ev(f_prime, x=x0)),
while abs(float(f(x1))) > tol do (
print("x0: ", float(x0), " x1: ", x1, " f(x1): ", float(f(x1))),
x0: x1,
x1: float(x0 − f(x0) / ev(f_prime, x=x0))
),
return(x1)
);
newton_raphson (f,x0 ,tol):=block([x1 ,f_prime],f_prime:

d f ( x0 )
f ( x ) , x1 : float x0 − , while float ( f ( x1 ) ) > tol
dx ev ( f_prime , x = x0 )

do ( print ( x0: , float ( x0 ) , x1: , x1 , f(x1): , float ( f ( x1 ) ) ) , x0 : x1 , x1 :

f ( x0 )
float x0 − ) , return ( x1 ) )
ev ( f_prime , x = x0 )

1 Find the root of sin(x^3 + x - 10).

f(x) :=sin(x^3 + x − 10);
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: newton_raphson(f, 0.5, 1e-4 );
print("Root: ", root);
3
f ( x ) := sin ( x + x − 10 )
Numerical [Link] 26 / 58

x0: 0.5 x1: 0.4715319339485794 f(x1): − 0.0011514197049267549

0.4708412310366765
Root: 0.4708412310366765 0.4708412310366765

2 Find the root of f(x) = x^3-x-1

f(x) :=x^3 − x − 1;
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: newton_raphson(f, 1.5, 1e-4);
print("Root: ", root);
3
f ( x ) := x − x − 1
plot2d: some values will be clipped.
Numerical [Link] 27 / 58

x0: 1.5 x1: 1.3478260869565217 f(x1): 0.10068217309114824

x0: 1.3478260869565217 x1: 1.325200398950907 f(x1):
0.0020583619166634204 1.3247181739990537
Root: 1.3247181739990537 1.3247181739990537

3 Find the root of f(x) = 2x^3-2x-5

f(x) := 2·x^3−2·x−5;
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: newton_raphson(f, 1.6, 1e-4);
print("Root: ", root);
3
f ( x ) := 2 x − 2 x − 5
plot2d: some values will be clipped.
Numerical [Link] 28 / 58

1.6005988023952096
Root: 1.6005988023952096 1.6005988023952096

4 Find the root of f(x) = x^3-x-1 with an

error of 10^(-5).
f(x) :=x^3−x−1;
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: newton_raphson(f, 1.5, 1e-5);
print("Root: ", root);
3
f ( x ) := x − x − 1
plot2d: some values will be clipped.
Numerical [Link] 29 / 58

x0: 1.5 x1: 1.3478260869565217 f(x1): 0.10068217309114824

x0: 1.3478260869565217 x1: 1.325200398950907 f(x1):
0.0020583619166634204 1.3247181739990537
Root: 1.3247181739990537 1.3247181739990537

5 Find the root of f(x) = x^4 + x^2 - 10

with an error of 10^(-5).
f(x) :=x^4 + x^2 − 10;
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: newton_raphson(f, 1.6, 1e-5);
print("Root: ", root);
4 2
f ( x ) := x + x − 10
plot2d: some values will be clipped.
Numerical [Link] 30 / 58

x0: 1.6 x1: 1.6452614379084967 f(x1): 0.034112680006114005

x0: 1.6452614379084967 x1: 1.643645080057159 f(x1):
−5
4.501708426829509 10 1.6436429413740434
Root: 1.6436429413740434 1.6436429413740434

6 Find the root of f(x) = cos x with error

0.00001
f(x) :=cos(x);
wxplot2d(f(x), [x,−5,5],[y,−5,5]);
root: newton_raphson(f, 0.75, 1e-5);
print("Root: ", root);
f ( x ) := cos ( x )
Numerical [Link] 31 / 58

x0: 0.75 x1: 1.8234261485493775 f(x1):− 0.2499511677645151

x0: 1.8234261485493775 x1: 1.5652810537151907 f(x1):
0.005515245118934408 1.5707963827172051
Root: 1.5707963827172051 1.5707963827172051

Trapezoidal Rule
trapezoidal_rule(f, a, b, n) := block(
[h, sum, i, x, I, err],
h: (b − a) / n,
sum: f(a) + f(b),

for i: 1 thru n−1 do (

sum: sum + 2 · f(a + i · h)
),

I: float(integrate(f(x), x, a, b)),
err: float(abs(I − h · sum / 2)),

print("Approximation: ",float(h · sum / 2)),

print("Actual Value: ", I),
print("Error: ", err)
);
Numerical [Link] 32 / 58

(f,a ,b ,n ):=block([h ,sum ,i,x,I,err],h : b − a ,

trapezoidal_rule
n
sum : f ( a ) + f ( b ) , for i thru n − 1 do sum : sum + 2 f ( a + i h ) , I :
b

h sum
float f ( x ) d x , err : float I − ,
2
a

h sum
print Approximation: , float , print ( Actual Value: , I ) ,
2
print ( Error: , err ) )

1 Evaluate f(x)=1/x in the interval [1,2].

f(x) := 1 / x;
trapezoidal_rule(f, 1, 2, 1);
1
f ( x ) :=
x
Approximation: 0.75
Actual Value: 0.6931471805599453
Error: 0.056852819440054714 0.056852819440054714
Result: 0.056852819440054714 0.056852819440054714

2 Evaluate f(x)=e^(x^2) in [0,1].

f(x) := exp(x^2);
trapezoidal_rule(f, 0, 1, 1);
2
f ( x ) := exp ( x )
Approximation: 1.8591409142295225
Actual Value: 1.4626517459071817
Error: 0.3964891683223408 0.3964891683223408
Result: 0.3964891683223408 0.3964891683223408

3 Evaluate f(x)=1/(1+(x^2)) in [0,1]

f(x) :=1/(1+x^2);
trapezoidal_rule(f, 0, 1, 2);
1
f ( x ) :=
2
1+x
Approximation: 0.775
Actual Value: 0.7853981633974483
Error: 0.010398163397448257 0.010398163397448257
Numerical [Link] 33 / 58

4 Evaluate f(x)=1+x/1+x^3 in [0,1].

f(x) := (1+x)/(1+x^3);
trapezoidal_rule(f, 0, 1, 2);
1+x
f ( x ) :=
3
1+x
Approximation: 1.1666666666666667
Actual Value: 1.2091995761561454
Error: 0.04253290948947863 0.04253290948947863

5 Evaluate f(x)=2^x in [0,4].

f(x) := 2^x;
trapezoidal_rule(f, 0, 4, 4);
x
f ( x ) := 2
Approximation: 22.5
Actual Value: 21.64042561333445
Error: 0.8595743866655496 0.8595743866655496

Simpson's 1/3 rule

simpsons_rule(f, a, b, n) := block(
[h, sum_even, sum_odd, i, x, I, err],

h: (b − a) / n,
sum_even: 0,
sum_odd: 0,

for i: 1 step 2 thru n−1 do (

sum_odd: sum_odd + f(a + i · h)
),
for i: 2 step 2 thru n−2 do (
sum_even: sum_even + f(a + i · h)
),

estimated_value: (h / 3) · (f(a) + 4 · sum_odd + 2 · sum_even + f(b)),

I: float(integrate(f(x), x, a, b)),
err: abs(I − estimated_value),

print("Actual Value: ", I),

print("Estimated Value (Simpson's Rule): ", estimated_value),
print("Error: ", err)
);
Numerical [Link] 34 / 58

simpsons_rule(f,a ,b ,n ):=block(

[h ,sum_even ,sum_odd ,i,x,I,err],h : b − a , sum_even : 0 , sum_odd : 0 ,

n
for i step 2 thru n − 1 do sum_odd : sum_odd + f ( a + i h ) , for i from
2 step 2 thru n − 2 do sum_even : sum_even + f ( a + i h ) ,
h
estimated_value : ( f ( a ) + 4 sum_odd + 2 sum_even + f ( b ) ) , I :
3
b

float f ( x ) d x , err : I − estimated_value ,

print ( Actual Value: , I ) ,

print ( Estimated Value (Simpson's Rule): , estimated_value ) ,
print ( Error: , err ) )

1 Evaluate f(x)=1/x for[1,2], using

simpson's 1/3 rule.
f(x) := 1 / x;
simpsons_rule(f, 1, 2, 2);
1
f ( x ) :=
x
Actual Value: 0.6931471805599453
25
Estimated Value (Simpson's Rule):
36
Error: 0.0012972638844991335 0.0012972638844991335

2 Evaluate f(x)=e^(x^2) in [0,1] using

simpson's 1/3 rule.
f(x) := exp(x^2);
simpsons_rule(f, 0, 1, 2);
2
f ( x ) := exp ( x )
Actual Value: 1.4626517459071817
1/4
%e + 4 %e +1
Estimated Value (Simpson's Rule):
6
1/4
%e + 4 %e +1
Error: − 1.4626517459071817
6
Numerical [Link] 35 / 58

1/4
%e + 4 %e +1
− 1.4626517459071817
6

3 Evaluate f(x)=1/(1+(x^2)) in [0,1]

using simpson's 1/3 rule.
f(x) :=1/(1+(x^2));
simpsons_rule(f, 0, 1, 2);
1
f ( x ) :=
2
1+x
Actual Value: 0.7853981633974483
47
Estimated Value (Simpson's Rule):
60
Error: 0.002064830064114953 0.002064830064114953

4 Evaluate f(x)=1/(1+x) in [0,1] using

simpson's 1/3 rule.
f(x) := 1 /(1+ x);
simpsons_rule(f, 0, 1, 2);
1
f ( x ) :=
1+x
Actual Value: 0.6931471805599453
25
Estimated Value (Simpson's Rule):
36
Error: 0.0012972638844991335 0.0012972638844991335

Euler's Method
euler_method(f, x0, y0, h, x_end) := block(
[x: float(x0), y: float(y0), steps: floor((x_end − x0) / h)],
print("Iteration: 0", " x: ", x, " y: ", y),
for i: 1 thru steps do (
y: y + h · f(x, y),
x: x + h,
print("Iteration: ", i, " x: ", x, " y: ", y)
)
);

euler_method ( f , x0 , y0 , h , x ) := block (
end
x − x0
x : float ( x0 ) , y : float ( y0 ) , steps : floor end ,
h
Numerical [Link] 36 / 58

print ( Iteration: 0 , x: , x , y: , y ) , for i thru steps do

( y : y + h f ( x , y ) , x : x + h , print ( Iteration: , i , x: , x , y: , y ) ) )

1 Solve the IVP dy/dx=x^2+y given

x0=0,y0=1 . Find the value at f(0.4)
with step size of h=0.1
f(x, y) := x^2 + y;
euler_method(f1, 0, 1, 0.1, 0.4);
2
f ( x , y ) := x + y
Iteration: 0 x: 0.0 y: 1.0
Iteration: 1 x: 0.1 y: 1.1
Iteration: 2 x: 0.2 y: 1.211
Iteration: 3 x: 0.30000000000000004 y: 1.3361
Iteration: 4 x: 0.4 y: 1.47871 done

2 Solve the IVP dy/dx=x+2*y given

x0=0,y0=0 . Find the value at f(1) with
step size of h=0.25.
f(x, y) := x + 2 · y;
euler_method(f2, 0, 0, 0.25, 1);
f ( x , y ) := x + 2 y
Iteration: 0 x: 0.0 y: 0.0
Iteration: 1 x: 0.25 y: 0.0
Iteration: 2 x: 0.5 y: 0.0625
Iteration: 3 x: 0.75 y: 0.21875
Iteration: 4 x: 1.0 y: 0.515625 done

3 Solve the IVP dy/dx=3*e^(-x)-0.4*y

given x0=0,y0=5 . Find the value at
f(3) with step size of h=1.5
f(x, y) := 3 · exp(−x) − 0.4 · y;
euler_method(f3, 0, 5, 1.5, 3);
f ( x , y ) := 3 exp ( − x ) − 0.4 y
Iteration: 0 x: 0.0 y: 5.0
Iteration: 1 x: 1.5 y: 6.5
Iteration: 2 x: 3.0 y: 3.604085720667934 done

RK 4th Order Method

Numerical [Link] 37 / 58

rk4_method(f, x0, y0, h, x_end) := block(

[x: float(x0), y: float(y0), steps: floor((x_end − x0) / h)],

print("Iteration: 0", " x: ", x, " y: ", y),

for i: 1 thru steps do (

k1: h · f(x, y),
k2: h · f(x + h / 2, y + k1 / 2),
k3: h · f(x + h / 2, y + k2 / 2),
k4: h · f(x + h, y + k3),
y: y + (k1 + 2 · k2 + 2 · k3 + k4) / 6,
x: x + h,
print("Iteration: ", i, " x: ", x, " y: ", y)
)
);
rk4_method (f,x0 ,y0 ,h ,x ):= block (
end
x − x0
x : float ( x0 ) , y : float ( y0 ) , steps : floor end ,
h
print ( Iteration: 0 , x: , x , y: , y ) , for i thru steps do ( k1 : h f ( x , y ) , k2 : h
h k1 h k2
f x+ ,y+ , k3 : h f x + ,y+ , k4 : h f ( x + h , y + k3 ) , y : y +
2 2 2 2
k1 + 2 k2 + 2 k3 + k4
, x : x + h , print ( Iteration: , i , x: , x , y: , y ) ) )
6

1 Solve the IVP dy/dx=x^2+y given

x0=0,y0=1 . Find the value at f(0.4)
with step size of h=0.1
f1(x, y) := x^2 + y;
rk4_method(f1, 0, 1, 0.1, 0.4);
2
f1 ( x , y ) := x + y
Iteration: 0 x: 0.0 y: 1.0
Iteration: 1 x: 0.1 y: 1.240549527280918
Iteration: 2 x: 0.2 y: 1.4450442806532784
Iteration: 3 x: 0.30000000000000004 y: 1.617431472124734
Iteration: 4 x: 0.4 y: 1.7612624484183106 done
Numerical [Link] 38 / 58

2 Solve the IVP dy/dx= x + 2 * y given

x0=0,y0=0 . Find the value at f(0.25)

with step size of

h=1

f2(x, y) := x + 2 · y;
rk4_method(f2, 0, 0, 0.25, 1);
f2 ( x , y ) := x + 2 y
Iteration: 0 x: 0.0 y: 0.0
Iteration: 1 x: 0.25 y: 0.6301814371757238
Iteration: 2 x: 0.5 y: 1.060997592910006
Iteration: 3 x: 0.75 y: 1.3422547723035547
Iteration: 4 x: 1.0 y: 1.5121990856142926 done

Gauss Jacobi Method

Numerical [Link] 39 / 58

gauss_jacobi(A, b, x0, tol, max_iter) := block(

[n, x, x_new, i, j, sum, err, iter],

n: length(b),
x: x0,
iter: 0,

while iter < max_iter do (

x_new: makelist(0, i, 1, n),

for i: 1 thru n do (
sum: 0,

for j: 1 thru n do (
if i # j then sum: sum + A[i][j] · x[j]
),

x_new[i]: (b[i] − sum) / A[i][i]

),

err: max(map(lambda([v, v_new], abs(v − v_new)), x, x_new)),

print("Iteration: ", iter + 1, " x_new: ", x_new, " Error: ", err),

if err < tol then return(x_new),

x: x_new,
iter: iter + 1
),

return("Solution did not converge within the maximum iterations.")

);

gauss_jacobi(A ,b ,x0 ,tol,max_iter):=block(

n ,x,x , i , j , sum , err , iter , n : length ( b ) , x : x0 , iter : 0 , while iter <

new
max_iter do ( x : makelist ( 0 , i , 1 , n ) , for i thru n do ( sum : 0 , for j
new
b − sum
i
thru n do if i ≠ j then sum : sum + [ A , array ] x , x : ),
Numerical [Link] 40 / 58

1 Solve system of linear equation using

Gauss Jacobi method

3x1+x2+x3=1

x1+3x2+x3=1

x1+x2+3x3=1

gauss_jacobi(
[[3, 1, 1],
[1, 3, 1],
[1, 1, 3]],
[1, 1, 1],
[0, 0, 0],
1e-4,
10
);

1 1 1 1 1 1
Iteration: 1 x_new: , , Error: max , ,
3 3 3 3 3 3

1 1 1 2 2 2
Iteration: 2 x_new: , , Error: max , ,
9 9 9 9 9 9

7 7 7 4 4 4
Iteration: 3 x_new: , , Error: max , ,
27 27 27 27 27 27

13 13 13 8 8 8
Iteration: 4 x_new: , , Error: max , ,
81 81 81 81 81 81

55 55 55
Iteration: 5 x_new: , , Error:
243 243 243

16 16 16
max , ,
243 243 243

133 133 133

Iteration: 6 x_new: , , Error:
729 729 729

32 32 32
max , ,
729 729 729

463 463 463

Iteration: 7 x_new: , , Error:
2187 2187 2187

64 64 64
max , ,
2187 2187 2187
Numerical [Link] 41 / 58

1261 1261 1261

Iteration: 8 x_new: , , Error:
6561 6561 6561

128 128 128

max , ,
6561 6561 6561

4039 4039 4039

Iteration: 9 x_new: , , Error:
19683 19683 19683

256 256 256

max , ,
19683 19683 19683

11605 11605 11605

Iteration: 10 x_new: , , Error:
59049 59049 59049

512 512 512

max , ,
59049 59049 59049
Solution did not converge within the maximum iterations.

2 Solve system of linear eqiuation using

Gauss Jacobi Method
10x1+3x2+4x3=17
3x1+6x2+x3=10
x1+x2+x3=3
gauss_jacobi(
[[10, 3, 4],
[3, 6, 1],
[1, 1, 1]],
[17, 10, 3],
[0, 0, 0],
1e-4,
10
);

17 5 17 5
Iteration: 1 x_new: , ,3 Error: max , ,3
10 3 10 3

19 11
Iteration: 2 x_new: 0, ,− Error:
60 30

17 27 101
max , ,
10 20 30

1051 311 161

Iteration: 3 x_new: , , Error:
600 180 60
Numerical [Link] 42 / 58

1051 127 61
max , ,
600 90 20

13 1237 863
Iteration: 4 x_new: , ,− Error:
120 3600 1800

493 1661 5693

max , ,
300 1200 1800

64393 9139 9173

Iteration: 5 x_new: , , Error:
36000 5400 3600

60493 14567 1211

max , ,
36000 10800 400

623 75091 51959

Iteration: 6 x_new: , ,− Error:
3600 216000 108000

58163 96823 327149

max , ,
36000 72000 108000

3862399 1075889 535529

Iteration: 7 x_new: , , Error:
2160000 648000 216000

3488599 106327 213149

max , ,
2160000 81000 72000

90799 4657513 2906087

Iteration: 8 x_new: , ,− Error:
432000 12960000 6480000

852101 1873363 18971957

max , ,
540000 1440000 6480000

229596157 15905033 31498517

Iteration: 9 x_new: , , Error:
129600000 9720000 12960000

202356457 49647593 12436897

max , ,
129600000 38880000 4320000

153529 292226359 158589791

Iteration: 10 x_new: , ,−
648000 777600000 388800000

198890357 326725427 1103545301

Error: max , ,
129600000 259200000 388800000
Solution did not converge within the maximum iterations.
Numerical [Link] 43 / 58

3 Solve system of linear equation using

Gauss-Jacobi method.

8x1 + x2 -x3 = 2

-x1 + 7x2 - 2x3 = 4

2x1 + x2 + 9x3 = 12

gauss_jacobi(
[[8, 1, −1],
[−1, 7, −2],
[2, 1, 9]],
[2, 4, 12],
[0, 0, 0],
1e-4,
10
);

1 4 4 1 4 4
Iteration: 1 x_new: , , Error: max , ,
4 7 3 4 7 3

29 83 17 2 5 5
Iteration: 2 x_new: , , Error: max , ,
84 84 14 21 12 42

187 569 289

Iteration: 3 x_new: , , Error:
672 588 252

15 1 17
max , ,
224 49 252

961 13249 8213

Iteration: 4 x_new: , , Error:
3528 14112 7056

83 407 121
max , ,
14112 14112 7056

3489 3881 7067

Iteration: 5 x_new: , , Error:
12544 4116 6048

649 401 191

max , ,
112896 98784 42336

659555 2240459 1383083

Iteration: 6 x_new: , , Error:
2370816 2370816 1185408

67 5003 683
max , ,
1185408 2370816 395136
Numerical [Link] 44 / 58

752477 15675151 8296741

Iteration: 7 x_new: , , Error:
2709504 16595712 7112448

9101 4031 251

max , ,
18966528 8297856 1016064

55313003 376148209 25815637

Iteration: 8 x_new: , , Error:
199148544 398297088 22127616

11887 55415 31985

max , ,
398297088 398297088 199148544

885127433 1316588645 1394054945

Iteration: 9 x_new: , ,
3186376704 1394039808 1194891264

1895 46609 10547

Error: max , ,
50577408 929359872 1194891264

18587330441 63196781867 39033073459

Iteration: 10 x_new: , ,
66913910784 66913910784 33456955392

86413 526907 465001

Error: max , ,
16728477696 66913910784 33456955392
Solution did not converge within the maximum iterations.

4 Solve system of linear equation using

Gauss-Jacobi method
5x1 - 2x2 + 3x3 = -1
-3x1 + 9x2 + x3 = 2
2x1 - x2 - 7x3 = 3
gauss_jacobi(
[[5, −2, 3],
[−3, 9, 1],
[2, −1, −7]],
[−1, 2, 3],
[0, 0, 0],
1e-4,
10
);

1 2 3 1 2 3
Iteration: 1 x_new: − , ,− Error: max , ,
5 9 7 5 9 7
Numerical [Link] 45 / 58

46 64 163
Iteration: 2 x_new: , ,− Error:
315 315 315

109 2 4
max , ,
315 105 45

302 133 131

Iteration: 3 x_new: , ,− Error:
1575 405 315

8 71 32
max , ,
175 567 315

2564 673 41744

Iteration: 4 x_new: , ,− Error:
14175 2025 99225

22 8 479
max , ,
2025 2025 99225

91961 294038 14036

Iteration: 5 x_new: , ,− Error:
496125 893025 33075

2221 551 52
max , ,
496125 178605 14175

831967 3353 1887181

Iteration: 6 x_new: , ,− Error:
4465125 10125 4465125

4318 8483 1097

max , ,
4465125 4465125 637875

1384588 13313332 13210114

Iteration: 7 x_new: , ,− Error:
7441875 40186125 31255875

6071 211 17
max , ,
22325625 1607445 3472875

261756851 465838364 594590783

Iteration: 8 x_new: , ,−
1406514375 1406514375 1406514375

69719 42752 19379

Error: max , ,
1406514375 468838125 200930625

1308934702 598984298 594552541

Iteration: 9 x_new: , ,−
7032571875 1808375625 1406514375

50149 68962 38242

Error: max , ,
2344190625 2531725875 1406514375
Iteration: 10 x_new:
Numerical [Link] 46 / 58

11780069404 20964710561 187283066419

, ,− Error:
63293146875 63293146875 443052028125

342914 260131 983996

max , ,
63293146875 63293146875 443052028125
Solution did not converge within the maximum iterations.

Gauss Seidel Method

gauss_seidel(A, b, x0, tol, max_iter) := block(
[n, x, i, j, sum, err, iter],

n: length(b),
x: x0,
iter: 0,

while iter < max_iter do (

for i: 1 thru n do (
sum: 0,
for j: 1 thru n do (
if i # j then sum: sum + A[i][j] · x[j]
),
x[i]: (b[i] − sum) / A[i][i]
),
err: max(map(lambda([v, v_new], abs(v − v_new)), x0, x)),
print("Iteration: ", iter + 1, " x: ", x, " Error: ", err),

if err < tol then return(x),

x0: x,
iter: iter + 1
),

return("Solution did not converge within the maximum iterations.")

);
gauss_seidel ( A , b , x0 , tol , max_iter ) := block (
[ n , x , i , j , sum , err , iter ] , n : length ( b ) , x : x0 , iter : 0 , while iter < max_iter
do ( for i thru n do ( sum : 0 , for j thru n do if i ≠ j then sum : sum
b − sum
i
+ [ A , array ] x , x : ) , err :
j j i [ A , array ]
i

max map lambda v , vnew , v − v , x0 , x ,

new
Numerical [Link] 47 / 58

print ( Iteration: , iter + 1 , x: , x , Error: , err ) , if err < tol then return ( x ) ,
x0 : x , iter : iter + 1 ) , return
( Solution did not converge within the maximum iterations. ) )

1 Solve system of linear eqiuation using

Gauss Seidel Method
5x1+x2+2x3=10
-3x1+9x2+4x3=14
-x1-x2+7x3=33
gauss_seidel(
[[5, 1, 2],
[−3, 9, 4],
[−1, −1, 7]],
[10, 14, 33],
[0, 0, 0],
1e-4,
10
);

20 335
Iteration: 1 x: 2, , Error: max ( [ 0 , 0 , 0 ] )
9 63

4 566 17821
Iteration: 2 x: − ,− , Error: max ( [ 0 , 0 , 0 ] )
7 567 3969

178 10912 1182299

Iteration: 3 x: ,− , Error: max ( [ 0 , 0 , 0 ] )
441 35721 250047

23584 5505922 367724177

Iteration: 4 x: ,− , Error:
138915 11252115 78764805

max ( [ 0 , 0 , 0 ] )

403294 313744748 23242736407

Iteration: 5 x: ,− , Error:
1750329 708883245 4962182715

max ( [ 0 , 0 , 0 ] )

592507508 101504235374 7315330204549

Iteration: 6 x: ,− ,
2756768175 223298222175 1563087555225

Error: max ( [ 0 , 0 , 0 ] )
Iteration: 7 x:
38016550906 6351648170584 460964696354627
,− , Error:
173676395025 14067787997025 98474515979175
Numerical [Link] 48 / 58

max ( [ 0 , 0 , 0 ] )
11919700475176 2004245895720538
Iteration: 8 x: [ ,− ,
54708064432875 4431353219062875
145195906071843593
] Error: max ( [ 0 , 0 , 0 ] )
31019472533440125
751836321350174 126211428142757924
Iteration: 9 x: [ ,− ,
3446608059271125 279175252800961125
1829494131263664659
] Error: max ( [ 0 , 0 , 0 ] )
390845353921345575
236756264492215292
Iteration: 10 x: [ ,−
1085681538670404375
39761120112605701766 2881442890177254631261
, ] Error:
87940204632302754375 615581432426119280625
max ( [ 0 , 0 , 0 ] )
Solution did not converge within the maximum iterations.

2 Solve system of linear eqiuation using

Gauss Seidel Method
10x1+3x2+4x3=17
3x1+6x2+x3=10
x1+x2+x3=3
gauss_seidel(
[[10, 3, 4],
[3, 6, 1],
[1, 1, 1]],
[17, 10, 3],
[0, 0, 0],
1e-4,
10
);

17 49 29
Iteration: 1 x: , , Error: max ( [ 0 , 0 , 0 ] )
10 60 60
757 3439 2819
Iteration: 2 x: , , Error: max ( [ 0 , 0 , 0 ] )
600 3600 3600
39607 212989 197369
Iteration: 3 x: , , Error: max ( [ 0 , 0 , 0 ] )
36000 216000 216000
Numerical [Link] 49 / 58

2243557 12895639 12523019

Iteration: 4 x: , , Error:
2160000 12960000 12960000
max ( [ 0 , 0 , 0 ] )
131541007 776146789 767407169
Iteration: 5 x: , , Error:
129600000 777600000 777600000
max ( [ 0 , 0 , 0 ] )
7821130957 46622535439 46418678819
Iteration: 6 x: , , Error:
7776000000 46656000000 46656000000
max ( [ 0 , 0 , 0 ] )
467609678407 2798584176589 2793837752969
Iteration: 7 x: , ,
466560000000 2799360000000 2799360000000
Error: max ( [ 0 , 0 , 0 ] )
Iteration: 8 x:
28018016458357 167943573095239 167833128154619
, , Error:
27993600000000 167961600000000 167961600000000
max ( [ 0 , 0 , 0 ] )
1680183968095807 10077276814166389
Iteration: 9 x: [ , ,
1679616000000000 10077696000000000
10074707377258769
] Error: max ( [ 0 , 0 , 0 ] )
10077696000000000
100790172048465757 604652010082015039
Iteration: 10 x: [ , ,
100776960000000000 604661760000000000
604592237627190419
] Error: max ( [ 0 , 0 , 0 ] )
604661760000000000
Solution did not converge within the maximum iterations.

3 Solve system of linear eqiuation using

Gauss Seidel Method
3x1+x2+x3=1
x1+3x2+x3=1
x1+x2+3x3=1
Numerical [Link] 50 / 58

gauss_seidel(
[[3, 1, 1],
[1, 3, 1],
[1, 1, 3]],
[1, 1, 1],
[0, 0, 0],
1e-4,
10
);

Iteration:1 x: 1 2 4
, , Error: max ( [ 0 , 0 , 0 ] )
3 9 27
17 52 140
Iteration: 2 x: , , Error: max ( [ 0 , 0 , 0 ] )
81 243 729
433 1334 3928
Iteration: 3 x: , , Error: max ( [ 0 , 0 , 0 ] )
2187 6561 19683
11753 35512 106376
Iteration: 4 x: , , Error: max ( [ 0 , 0 , 0 ] )
59049 177147 531441
318529 956666 2870716
Iteration: 5 x: , , Error:
1594323 4782969 14348907
max ( [ 0 , 0 , 0 ] )
8608193 25826380 77489204
Iteration: 6 x: , , Error:
43046721 129140163 387420489
max ( [ 0 , 0 , 0 ] )
232452145 697341710 2092086256
Iteration: 7 x: , , Error:
1162261467 3486784401 10460353203
max ( [ 0 , 0 , 0 ] )
6276241817 18828559024 56485894352
Iteration: 8 x: , ,
31381059609 94143178827 282429536481
Error: max ( [ 0 , 0 , 0 ] )
169457965057 508372961330 1525118971828
Iteration: 9 x: , ,
847288609443 2541865828329 7625597484987
Error: max ( [ 0 , 0 , 0 ] )
Iteration: 10 x:
4575359629169 13726075910308 41178222567068
, , Error:
22876792454961 68630377364883 205891132094649
max ( [ 0 , 0 , 0 ] )
Solution did not converge within the maximum iterations.
Numerical [Link] 51 / 58

4 Solve system of linear equation using

Gauss-Seidel method

5x1 - 2x2 + 3x3 = -1

-3x1 + 9x2 + x3 = 2

2x1 - x2 - 7x3 = 3

gauss_seidel(
[[5, −2, 3],
[−3, 9, 1],
[2, −1, −7]],
[−1, 2, 3],
[0, 0, 0],
1e-4,
10
);

1 7 32
Iteration: 1 x: − , ,− Error: max ( [ 0 , 0 , 0 ] )
5 45 63

263 677 8506

Iteration: 2 x: , ,− Error: max ( [ 0 , 0 , 0 ] )
1575 2025 19845

94711 212719 2635922

Iteration: 3 x: , ,− Error:
496125 637875 6251175

max ( [ 0 , 0 , 0 ] )

29129417 66549293 832211734

Iteration: 4 x: , ,− Error:
156279375 200930625 1969120125

max ( [ 0 , 0 , 0 ] )

9159406099 20962788271 262205258498

Iteration: 5 x: , ,−
49228003125 63293146875 620272839375

Error: max ( [ 0 , 0 , 0 ] )
Iteration: 6 x:
2886067931153 6603853889237 82592384208406
, ,− Error:
15506820984375 19937341265625 195385944403125

max ( [ 0 , 0 , 0 ] )
909133722255691 2080215427447639
Iteration: 7 x: [ , ,−
4884648610078125 6280262498671875
26016522692728082
] Error: max ( [ 0 , 0 , 0 ] )
61546572486984375
Numerical [Link] 52 / 58

286376089845867977 655267137315003533
Iteration: 8 x: [ , ,−
1538664312174609375 1978282687081640625
8195207354538770854
] Error: max ( [ 0 , 0 , 0 ] )
19387170333400078125
90208438107951518419 206409144979618247551
Iteration: 9 x: [ ,
484679258335001953125 623159046430716796875
2581490420791850627138
,− ] Error: max ( [ 0 , 0 , 0 ] )
6106958655021024609375
28415659244775224620193
Iteration: 10 x: [ ,
152673966375525615234375
65018881572167595715397 813169479347682147634486
,− ]
196295099625675791015625 1923691976331622751953125
Error: max ( [ 0 , 0 , 0 ] )
Solution did not converge within the maximum iterations.

5 Solve system of linear equation using

Gauss-Seidel method.
8x1 + x2 -x3 = 2
-x1 + 7x2 - 2x3 = 4
2x1 + x2 + 9x3 = 12
gauss_seidel(
[[8, 1, −1],
[−1, 7, −2],
[2, 1, 9]],
[2, 4, 12],
[0, 0, 0],
1e-4,
10
);

1 17 305
Iteration: 1 x: , , Error: max ( [ 0 , 0 , 0 ] )
4 28 252
41 425 4580
Iteration: 2 x: , , Error: max ( [ 0 , 0 , 0 ] )
126 441 3969
8693 208981 2336485
Iteration: 3 x: , , Error: max ( [ 0 , 0 , 0 ] )
31752 222264 2000376
557051 13231525 147001345
Iteration: 4 x: , , Error:
2000376 14002632 126023688
Numerical [Link] 53 / 58

max ( [ 0 , 0 , 0 ] )
69991249 1666186133 18525915965
Iteration: 5 x: , , Error:
252047376 1764331632 15878984688
max ( [ 0 , 0 , 0 ] )
1102756567 52489398475 583550774545
Iteration: 6 x: , ,
3969746172 55576446408 500188017672
Error: max ( [ 0 , 0 , 0 ] )
Iteration: 7 x:
555761111807 13227175590919 147055342781575
, , Error:
2000752070688 14005264494816 126047380453344
max ( [ 0 , 0 , 0 ] )
Iteration: 8 x:
35013190421249 833313397797775 9264481894385635
, , Error:
126047380453344 882331663173408 7940984968560672
max ( [ 0 , 0 , 0 ] )
4411657812831751 104997465138859667
Iteration: 9 x: [ , ,
15881969937121344 111173789559849408
1167324800199688715
] Error: max ( [ 0 , 0 , 0 ] )
1000564106038644672
69483619563351283 3307420251403680625
Iteration: 10 x: [ , ,
250141026509661168 3501974371135256352
36770730854445559675
] Error: max ( [ 0 , 0 , 0 ] )
31517769340217307168
Solution did not converge within the maximum iterations.

Langrange Interpolation
Numerical [Link] 54 / 58

lagrange_interpolation(x_vals, y_vals, eval_point) := block(

[n, L, P, x, i, j, result],
n: length(x_vals),
P: 0,

for i: 1 thru n do (
L: 1,
for j: 1 thru n do (
if i # j then L: expand(L · (x − x_vals[j]) / (x_vals[i] − x_vals[j]))
),
P: expand(P + y_vals[i] · L)
),

if not is(eval_point = false) then (

result: ev(P, x = eval_point),
return([P, result])
) else (
return(P)
)
);

lagrange_interpolation
(x , y , eval_point ):= block (
vals vals
[ n , L , P , x , i , j , result ] , n : length ( x ) , P : 0 , for i thru n do ( L : 1 , for j
vals

L x − xvals
j
thru n do if i ≠ j then L : expand ,P:
x −x
vals vals
i j

expand P + yvals L ) , if not is ( eval_point = false ) then

i

( result : ev ( P , x = eval_point ) , return ( [ P , result ] ) ) else return ( P ) )

1 Determine the polynomial

f(1)=5
f(2)=10
f(3)=15
x_vals: [1, 2, 3];
y_vals: [5, 10,15];
P: lagrange_interpolation(x_vals, y_vals,1);
[1,2,3] [ 5 , 10 , 15 ] [5 x,5]
Numerical [Link] 55 / 58

2 Determine the polynomial

f(1)=1

f(3)=27

f(5)=55

Also find the value of f(2.5).

x_vals: [1, 3, 5];

y_vals: [1, 27,55];
P: lagrange_interpolation(x_vals, y_vals,2.5);
2
[1,3,5] [ 1 , 27 , 55 ] x 45
+ 12 x − , 20.3125
4 4

3 Find the interpolating polynomial

using Lagrange interpolation given
x=[0 2 4] and y=[2.7 3.4 4.2].
EVALUATE f(3.0).
x_vals: [0, 2,4];
y_vals: [2.7,3.4,4.2];
P: lagrange_interpolation(x_vals, y_vals,3.0);
[0,2,4] [ 2.7 , 3.4 , 4.2 ] [ 0.012500000000000067
2
x + 0.3249999999999995 x + 2.7 , 3.7874999999999996 ]

Newton Divided Difference

Interpolation
Numerical [Link] 56 / 58

divided_difference(x_values, y_values) := block(

[n: length(x_values), dd: []],
/* Initialize the divided difference table */
dd: makelist(y_values[i], i, 1, n),
for i: 2 thru n do (
for j: n step −1 thru i do (
dd[j]: (dd[j] − dd[j−1]) / (x_values[j] − x_values[j−i+1])
)
),
return(dd)
);

/* Define function for Newton Polynomial */

newton_polynomial(x_values, dd, x) := block(
[n: length(x_values), poly: dd[1], term: 1],
for i: 2 thru n do (
term: term · (x − x_values[i−1]),
poly: poly + dd[i] · term
),
return(expand(poly))
);
divided_difference
(x ,y ):= block (
values values

n : length ( x ) , dd : [ ] , dd : makelist yvalues , i , 1 , n , for i from 2

values i
dd − dd
j j−1
thru n do for j from n step − 1 thru i do dd :
j x −x
values values
j j−i+1

, return ( dd ) ) newton_polynomial ( x , dd , x ) := block (

values
n : length ( x ) , poly : dd , term : 1 , for i from 2 thru n do
values 1

term : term x − xvalues , poly : poly + dd term ,

i−1 i

return ( expand ( poly ) ) )

1 Determine the polynomial

f(1)=1
f(3=27
f(5)=55
Also find the value of f(2.5).
Numerical [Link] 57 / 58

x_vals: [1, 3, 5];

y_vals: [1, 27, 55];
dd: divided_difference(x_vals, y_vals);
polynomial: newton_polynomial(x_vals, dd, x);
value_at_2_5: ev(polynomial, x = 2.5);
print("Interpolating Polynomial: ", polynomial);
print("f(2.5): ", value_at_2_5);

2
[1 ,3 ,5 ] [1 ,27 ,55 ] 1 ,13 , 1 x
+ 12
4 4
45
x− 20.3125
4
2 2
x 45 x 45
Interpolating Polynomial: + 12 x − + 12 x −
4 4 4 4
f(2.5): 20.3125 20.3125

2 Find the interpolating polynomial

using divided difference interpolation
given x=[1 2 3] and y=[5 10 15].
EVALUATE f(1).
x_vals: [1, 2, 3];
y_vals: [5, 10, 15];
dd: divided_difference(x_vals, y_vals);
polynomial: newton_polynomial(x_vals, dd, x);
value_at_1: ev(polynomial, x = 1);
print("Interpolating Polynomial: ", polynomial);
print("f(1): ", value_at_1);

[1,2,3] [ 5 , 10 , 15 ] [5,5,0] 5x
5
Interpolating Polynomial: 5 x 5x
f(1): 5 5
Numerical [Link] 58 / 58

3 Find the interpolating polynomial

using divided difference

interpolation

given x=[0 1 2 3] and y=[1 2 4 8].

EVALUATE f(3.5).

x_vals: [0, 1, 2, 3];

y_vals: [1, 2, 4, 8];
dd: divided_difference(x_vals, y_vals);
polynomial: newton_polynomial(x_vals, dd, x);
value_at_3_5: ev(polynomial, x = 3.5);
print("Interpolating Polynomial: ", polynomial);
print("f(3.5): ", value_at_3_5);

3
1 1 x
[0,1,2,3] [1,2,4,8] 1,1, ,
2 6 6
5x
+ +1 11.0625
6
3 3
x 5x x 5x
Interpolating Polynomial: + +1 + +1
6 6 6 6
f(3.5): 11.0625 11.0625