initial commit / add all books

16 files changed, 568 insertions, 0 deletions

diff --git a/845/CH7/EX7.1/Ex7_1.sce b/845/CH7/EX7.1/Ex7_1.sce
new file mode 100755
index 000000000..999247d17
--- /dev/null
+++ b/845/CH7/EX7.1/Ex7_1.sce
@@ -0,0 +1,29 @@
+//Example 7.1

+

+clc

+clear

+

+x = 0:0.2:1;

+y = [1 1.16 3.56 13.96 41.96 101];

+

+n = length(x);

+del = %nan*ones(n,6);

+del(:,1) = y';

+for j = 2:6

+    for i = 1:n-j+1

+        del(i,j) = del(i+1,j-1) - del(i,j-1);

+    end

+end

+del = round(del*10^2)/10^2;

+mprintf("%5s %6s %9s %8s %8s %8s %7s",'x','y','dy','d2y','d3y','d4y','d5y')

+disp([x' del])

+

+h = x(2) - x(1);

+del0 = del(1,:);

+del1 = del(2,:);

+

+df1 = (del1(2) - del1(3)/2 + del1(4)/3 - del1(5)/4) / h;

+d2f0 = (del0(2) - del0(3) + del0(4)*11/12 - del0(5)*5/6) / h^2;

+disp(round(d2f0*10^1)/10^1,"f''''(0) = ")

+disp(round(df1*10)/10,"f''(0.2) = ")

+

diff --git a/845/CH7/EX7.10/Ex7_10.sce b/845/CH7/EX7.10/Ex7_10.sce
new file mode 100755
index 000000000..6c7e359e4
--- /dev/null
+++ b/845/CH7/EX7.10/Ex7_10.sce
@@ -0,0 +1,23 @@
+//Etample 7.10

+

+clc

+clear

+

+t = 0:10:80;

+a = [30 31.63 33.34 35.47 37.75 40.33 43.25 46.69 50.67];

+

+h = t(2) - t(1);

+n = length(t);

+

+Is13 = a(1);

+for i = 2:n-1

+    rem2 = i-fix(i./2).*2;

+    if rem2 == 0 then

+        Is13 = Is13 + 4*a(i);

+    else

+        Is13 = Is13 + 2*a(i);

+    end

+end

+Is13 = (Is13 + a(n))/10^3;

+Is13 = round(h/3*Is13*10^4)/10^4;

+disp(strcat(["v = ",string(Is13)," km/s"]))

diff --git a/845/CH7/EX7.11/Ex7_11.sce b/845/CH7/EX7.11/Ex7_11.sce
new file mode 100755
index 000000000..1ff7ec466
--- /dev/null
+++ b/845/CH7/EX7.11/Ex7_11.sce
@@ -0,0 +1,31 @@
+//Example 7.11

+

+clc

+clear

+

+x = 1:0.1:1.8;

+x = round(x*10)/10;

+y = [1.543 1.669 1.811 1.971 2.151 2.352 2.577 2.828 3.107];

+n = length(x);

+x0 = x(1);

+xn = x(n);

+

+N = [1 2 4 8]

+for j = 1:length(N)

+    h = (xn - x0)./N(j);

+    I = y(1);

+    for xx = x0+h:h:xn-h

+        xx = round(xx*10)/10;

+        I = I + 2*y(x==xx);

+    end

+    Itrap(j) = h/2*(I + y(n));

+    IRomb(1) = Itrap(1);

+    if j~=1 then

+        IRomb(j) = (4^(j-1)*Itrap(j)-IRomb(j-1))/(4^(j-1)-1);

+    end

+end

+IRomb = round(IRomb*10^5)/10^5;

+

+disp(Itrap(length(N)),"Integral using Trapezoidal rule:")

+disp(IRomb(length(N)),"Integral using Romberg''s formula:")

+//In third step of computation of integral using Romberg's formula, author mistakenly took the 1.7672 instead of 1.7684 which resulted in a difference

diff --git a/845/CH7/EX7.12/Ex7_12.sce b/845/CH7/EX7.12/Ex7_12.sce
new file mode 100755
index 000000000..bd73dac92
--- /dev/null
+++ b/845/CH7/EX7.12/Ex7_12.sce
@@ -0,0 +1,39 @@
+//Example 7.12

+

+clc

+clear

+

+function [f] = fun1(x,y)

+    f = 1 / (x+y);

+endfunction

+

+x = 1:0.25:2;

+y = x;

+

+m = length(x);

+n = length(y);

+

+del = %nan*ones(m,n);

+for j = 1:n

+    for i = 1:m

+        del(i,j) = fun1(x(i),y(j));

+    end

+end

+

+hx = x(2) - x(1);

+for i = 1:m

+    I = del(i,1);

+    for j = 2:n-1

+        I = I + 2*del(i,j);

+    end

+    Itrap1(i) = hx/2 * (I+del(i,n));

+end

+Itrap1 = round(Itrap1*10^4)/10^4;

+

+hy = y(2) - y(1);

+Itrap2 = Itrap1(1)

+for i = 2:n-1

+    Itrap2 = Itrap2 + 2* Itrap1(i);

+end

+Itrap2 = round(hy/2*(Itrap2+Itrap1(m))*10^4)/10^4;

+disp(Itrap2,"I = ")

diff --git a/845/CH7/EX7.13/Ex7_13.sce b/845/CH7/EX7.13/Ex7_13.sce
new file mode 100755
index 000000000..bfc8ba495
--- /dev/null
+++ b/845/CH7/EX7.13/Ex7_13.sce
@@ -0,0 +1,39 @@
+//Example 7.13

+

+clc

+clear

+

+function [f] = fun1(x,y)

+    f = sqrt(sin(x+y));

+endfunction

+

+x = 0:%pi/8:%pi/2;

+y = x;

+

+m = length(x);

+n = length(y);

+

+del = %nan*ones(m,n);

+for j = 1:n

+    for i = 1:m

+        del(i,j) = fun1(x(i),y(j));

+    end

+end

+

+hx = x(2) - x(1);

+for i = 1:m

+    I = del(i,1);

+    for j = 2:n-1

+        I = I + 2*del(i,j);

+    end

+    Itrap1(i) = hx/2 * (I+del(i,n));

+end

+Itrap1 = round(Itrap1*10^4)/10^4;

+

+hy = y(2) - y(1);

+Itrap2 = Itrap1(1)

+for i = 2:n-1

+    Itrap2 = Itrap2 + 2* Itrap1(i);

+end

+Itrap2 = round(hy/2*(Itrap2+Itrap1(m))*10^4)/10^4;

+disp(Itrap2,"I = ")

diff --git a/845/CH7/EX7.14/Ex7_14.sce b/845/CH7/EX7.14/Ex7_14.sce
new file mode 100755
index 000000000..47da1f96f
--- /dev/null
+++ b/845/CH7/EX7.14/Ex7_14.sce
@@ -0,0 +1,25 @@
+//Example 7.14

+

+clc

+clear

+

+n = 1;

+if n==1 then

+    M = [0 2];

+elseif n==2

+    M = [sqrt(1/3) 1; -sqrt(1/3) 1];

+elseif n==3

+    M = [0 8/9; -0.774597 5/9; 0.774597 5/9];

+elseif n==4

+    M = [-0.339981 0.652145; -0.861136 0.347855; 0.339981 0.652145; 0.861136 0.347855];

+elseif n==5

+    M = [-0 0.568889; -0.538469 0.467914; -0.906180 0.236927; 0 0.568889; 0.538469 0.467914; 0.906180 0.236927];

+elseif n==6

+    M = [-0.238619 0.467914; -0.661209 0.360762; -0.932470 0.171325; 0.238619 0.467914; 0.661209 0.360762; 0.932470 0.171325];

+end

+

+X = M(:,1);

+W = M(:,2);

+

+disp(X,"E1 = ")

+disp(W,"W1 = ")

diff --git a/845/CH7/EX7.15/Ex7_15.sce b/845/CH7/EX7.15/Ex7_15.sce
new file mode 100755
index 000000000..b3ec6213b
--- /dev/null
+++ b/845/CH7/EX7.15/Ex7_15.sce
@@ -0,0 +1,27 @@
+//Example 7.15

+

+clc

+clear

+

+n = 2;

+if n==1 then

+    M = [0 2];

+elseif n==2

+    M = [sqrt(1/3) 1; -sqrt(1/3) 1];

+elseif n==3

+    M = [0 8/9; -0.774597 5/9; 0.774597 5/9];

+elseif n==4

+    M = [-0.339981 0.652145; -0.861136 0.347855; 0.339981 0.652145; 0.861136 0.347855];

+elseif n==5

+    M = [-0 0.568889; -0.538469 0.467914; -0.906180 0.236927; 0 0.568889; 0.538469 0.467914; 0.906180 0.236927];

+elseif n==6

+    M = [-0.238619 0.467914; -0.661209 0.360762; -0.932470 0.171325; 0.238619 0.467914; 0.661209 0.360762; 0.932470 0.171325];

+end

+

+X = M(:,1);

+W = M(:,2);

+

+disp(W(1),"W1 = ")

+disp(W(2),"W2 = ")

+disp(X(1),"E1 = ")

+disp(X(2),"E2 = ")

diff --git a/845/CH7/EX7.16/Ex7_16.sce b/845/CH7/EX7.16/Ex7_16.sce
new file mode 100755
index 000000000..fcd56815d
--- /dev/null
+++ b/845/CH7/EX7.16/Ex7_16.sce
@@ -0,0 +1,31 @@
+//Example 7.16

+

+clc

+clear

+

+function [f] = fun1(x)

+    f = 3*x^2 + x^3;

+endfunction

+

+n = 4;

+if n==1 then

+    M = [0 2];

+elseif n==2

+    M = [sqrt(1/3) 1];

+elseif n==3

+    M = [0 8/9; -0.774597 5/9; 0.774597 5/9];

+elseif n==4

+    M = [-0.339981 0.652145; -0.861136 0.347855; 0.339981 0.652145; 0.861136 0.347855];

+elseif n==5

+    M = [-0 0.568889; -0.538469 0.467914; -0.906180 0.236927; 0 0.568889; 0.538469 0.467914; 0.906180 0.236927];

+elseif n==6

+    M = [-0.238619 0.467914; -0.661209 0.360762; -0.932470 0.171325; 0.238619 0.467914; 0.661209 0.360762; 0.932470 0.171325];

+end

+

+X = M(:,1);

+W = M(:,2);

+I = 0;

+for i = 1:length(X)

+    I = I + W(i)*fun1(X(i));

+end

+disp(I,"I = ")

diff --git a/845/CH7/EX7.2/Ex7_2.sce b/845/CH7/EX7.2/Ex7_2.sce
new file mode 100755
index 000000000..01a902a7d
--- /dev/null
+++ b/845/CH7/EX7.2/Ex7_2.sce
@@ -0,0 +1,28 @@
+//Example 7.2

+

+clc

+clear

+

+x = 1.4:0.2:2.2;

+y = [4.0552 4.953 6.0496 7.3891 9.025];

+

+n = length(x);

+del = %nan*ones(n,5);

+del(:,1) = y';

+for j = 2:5

+    for i = 1:n-j+1

+        del(i+j-1,j) = del(i+j-1,j-1) - del(i+j-2,j-1);

+    end

+end

+mprintf("%5s %6s %10s %10s %9s %9s",'x','y','dy','d2y','d3y','d4y')

+disp([x' del])

+

+h = x(2) - x(1);

+deln = del(5,:);

+

+dfn = (deln(2) + deln(3)/2 + deln(4)/3 + deln(5)/4) / h;

+d2fn = (deln(3) + deln(4) + deln(5)*11/12) / h^2;

+dfn = round(dfn*10^4)/10^4;

+d2fn = round(d2fn*10^4)/10^4;

+disp(dfn,"y''(2.2) = ")

+disp(d2fn,"y''''(2.2) = ")

diff --git a/845/CH7/EX7.3/Ex7_3.sce b/845/CH7/EX7.3/Ex7_3.sce
new file mode 100755
index 000000000..1d02509df
--- /dev/null
+++ b/845/CH7/EX7.3/Ex7_3.sce
@@ -0,0 +1,52 @@
+//Example 7.3

+

+clc

+clear

+

+x = 0:4;

+y = [6.9897 7.4036 7.7815 8.1281 8.451];

+

+n = length(x);

+del = %nan*ones(n,5);

+del(:,1) = y';

+for j = 2:6

+    for i = 1:n-j+1

+        del(i,j) = del(i+1,j-1) - del(i,j-1);

+    end

+end

+del(:,1) = [];

+n0 = length(del(1,:));

+

+X = 2;

+i = find(x==X);

+dowy = 0;

+

+for j = 1:n0

+    if j==2*int(j/2) then

+        add = del(i,j);

+    else

+        add = (del(i-1,j) + del(i,j))/2;

+        i = i-1;

+        if i==0 then

+            break

+        end

+    end

+

+    if add == %nan then

+        break

+    else

+        dowy(j) = add;

+    end

+end

+mprintf("%5s %6s %10s %9s %9s %9s",'x','y','dy','d2y','d3y','d4y')

+disp([x' y' del])

+

+mu = 1;

+h = x(2) - x(1);

+dy2 = mu/h*(dowy(1) - 1/6*dowy(3));

+d2y2 = mu/h^2*(dowy(2)-1/12*dowy(4));

+dy2 = round(dy2*10^4)/10^4;

+d2y2 = round(d2y2*10^4)/10^4;

+

+disp(dy2,"y''(2) = ")

+disp(d2y2,"y''''(2) = ")

diff --git a/845/CH7/EX7.4/Ex7_4.sce b/845/CH7/EX7.4/Ex7_4.sce
new file mode 100755
index 000000000..f35e51791
--- /dev/null
+++ b/845/CH7/EX7.4/Ex7_4.sce
@@ -0,0 +1,43 @@
+//Example 7.4

+

+clc

+clear

+

+x = [0.15 0.21 0.23 0.27 0.32 0.35];

+y = [0.1761 0.3222 0.3617 0.4314 0.5051 0.5441];

+

+n = length(x);

+del = %nan*ones(n,6);

+del(:,1) = y';

+for j = 2:6

+    for i = 1:n-j+1

+        del(i,j) = (del(i+1,j-1) - del(i,j-1)) / (x(i+j-1)-x(i));

+    end

+end

+del(:,1) = [];

+del = round(del*10^3)/10^3;

+mprintf("%5s %6s %10s %10s %8s %9s %9s",'x','y','dy','d2y','d3y','d4y','d5y')

+disp([x' y' del])

+X = poly(0, "X");

+del0 = del(1,:);

+y0 = y(1);

+Y = y0;

+for i = 1:length(del0)

+    p = 1;

+    for j = 1:i

+        p = p*(X-x(j));

+    end

+    Y = Y + p*del0(i);

+end

+

+dydx = derivat(Y);

+d2ydx = derivat(dydx);

+

+XX = 0.25;

+dy = horner(dydx,XX);

+d2y = horner(d2ydx,XX);

+

+disp(round(dy*10^4)/10^4,"y''(0.25) = ")

+disp(d2ydx,"y''''(x) = ")

+disp(d2y,"y''''(0.25) = ")

+//The constant term in y''(x) is incorrectly computed to -91.7 instead of -97.42 in the text.

diff --git a/845/CH7/EX7.5/Ex7_5.sce b/845/CH7/EX7.5/Ex7_5.sce
new file mode 100755
index 000000000..ebc95f32b
--- /dev/null
+++ b/845/CH7/EX7.5/Ex7_5.sce
@@ -0,0 +1,28 @@
+//Example 7.5

+

+clc

+clear

+

+function [f] = y(x)

+    f = -1/x;

+endfunction

+

+H = [0.0128 0.0064 0.0032];

+n = length(H);

+x = 0.05;

+h = H(1);

+Fh = (y(x+h) - y(x-h)) / (2*h);

+Fh2 = (y(x+h/2) - y(x-h/2)) / (h);

+Fh4 = (y(x+h/4) - y(x-h/4)) / (h/2);

+

+F1h2 = (4*Fh2 - Fh) / (4-1);

+F1h4 = (4*Fh4 - Fh2) / (4-1);

+F2h4 = (4^2*F1h4 - F1h2) / (4^2-1);

+del = %nan*ones(n,3);

+del(:,1) = [Fh Fh2 Fh4]';

+del(1:2,2) = [F1h2 F1h4]';

+del(1,3) = F2h4;

+

+disp(del(1,n),"y''(0.05) = ")

+Exact = 1/x^2;

+disp(Exact,"Exact Value:")

diff --git a/845/CH7/EX7.6/Ex7_6.sce b/845/CH7/EX7.6/Ex7_6.sce
new file mode 100755
index 000000000..d8b00fe69
--- /dev/null
+++ b/845/CH7/EX7.6/Ex7_6.sce
@@ -0,0 +1,52 @@
+//Example 7.6

+

+clc

+clear

+

+function [I] = trap (fun,a,b,n)

+// Integrate the function over the interval using Trapezoidal Formula

+// trap (fun,a,b,n)

+// fun - function to be integrated

+// a - lower limit of integration

+// b - upper limit of integration

+// n - No. of times trapezoidal rule needs to be performed

+

+N = n + 1;  // N - total no. of points

+h = (b-a) / (N-1);

+x = linspace(a,b,N);

+y = fun(x);

+

+sum1 = y(1) + 2 * sum(y(2:N-1)) + y(N);

+I = h * sum1 / 2;                   // Trapezoidal Integral Value

+endfunction

+

+function [I] = simp13 (fun,a,b,n)

+// Integrate the function over the interval using Simpson's 1/3rd rule

+// simp13 (fun,a,b,n)

+// fun - function to be integrated

+// a - lower limit of integration

+// b - upper limit of integration

+// n - No. of times simpson's 1/3rd rule needs to be performed

+

+N = 2 * n + 1;      // N - total no. of points

+h = (b-a) / (N-1);

+x = linspace(a,b,N);

+y = fun(x);

+

+sum1 = y(1) + 4 * sum(y(2:2:N-1)) + 2 * sum(y(3:2:N-2)) + y(N);

+I = h* sum1 / 3;                    // Simpson's 1/3rd Integral Value

+endfunction

+

+n = 6;

+ntrap = n;

+ns13 = n/2;

+I = [trap(sin,0,%pi,ntrap); simp13(sin,0,%pi,ns13)];

+I = round(I*10^4)/10^4;

+true = integrate('sin(x)','x',0,%pi);

+err = abs(true - I) / true*100;

+err = round(err*100)/100;

+

+disp(I(1),"y_trap = ")

+disp(I(2),"y_simp13 = ")

+disp(err(1),"error_trap = ")

+disp(err(2),"error_simp13 = ")

diff --git a/845/CH7/EX7.7/Ex7_7.sce b/845/CH7/EX7.7/Ex7_7.sce
new file mode 100755
index 000000000..df21dadbd
--- /dev/null
+++ b/845/CH7/EX7.7/Ex7_7.sce
@@ -0,0 +1,34 @@
+//Example 7.7

+

+clc

+clear

+

+function [I] = simp13 (fun,a,b,n)

+// Integrate the function over the interval using Simpson's 1/3rd rule

+// simp13 (fun,a,b,n)

+// fun - function to be integrated

+// a - lower limit of integration

+// b - upper limit of integration

+// n - No. of times simpson's 1/3rd rule needs to be performed

+

+N = 2 * n + 1;      // N - total no. of points

+h = (b-a) / (N-1);

+x = linspace(a,b,N);

+y = fun(x);

+

+sum1 = y(1) + 4 * sum(y(2:2:N-1)) + 2 * sum(y(3:2:N-2)) + y(N);

+I = h* sum1 / 3;                    // Simpson's 1/3rd Integral Value

+endfunction

+

+n = 8;

+ns13 = n/2;

+I = simp13(log,1,5,ns13);

+I = round(I*10^4)/10^4;

+deff('[y] = true(x)',['y = x * log(x) - x']);

+trueVal = true(5) - true(1);

+err = abs(trueVal - I) / trueVal*100;

+err = round(err*100)/100;

+

+disp(I,"y_simp13 = ")

+disp(trueVal,"Actual Integral = ")

+disp(err,"error_simp13 = ")

diff --git a/845/CH7/EX7.8/Ex7_8.sce b/845/CH7/EX7.8/Ex7_8.sce
new file mode 100755
index 000000000..0318b5ee7
--- /dev/null
+++ b/845/CH7/EX7.8/Ex7_8.sce
@@ -0,0 +1,54 @@
+//Example 7.8

+

+clc

+clear

+

+function [I] = trap (fun,a,b,n)

+// Integrate the function over the interval using Trapezoidal Formula

+// trap (fun,a,b,n)

+// fun - function to be integrated

+// a - lower limit of integration

+// b - upper limit of integration

+// n - No. of times trapezoidal rule needs to be performed

+

+N = n + 1;  // N - total no. of points

+h = (b-a) / (N-1);

+x = linspace(a,b,N);

+y = fun(x);

+

+sum1 = y(1) + 2 * sum(y(2:N-1)) + y(N);

+I = h * sum1 / 2;                   // Trapezoidal Integral Value

+endfunction

+

+function [I] = simp13 (fun,a,b,n)

+// Integrate the function over the interval using Simpson's 1/3rd rule

+// simp13 (fun,a,b,n)

+// fun - function to be integrated

+// a - lower limit of integration

+// b - upper limit of integration

+// n - No. of times simpson's 1/3rd rule needs to be performed

+

+N = 2 * n + 1;      // N - total no. of points

+h = (b-a) / (N-1);

+x = linspace(a,b,N);

+y = fun(x);

+

+sum1 = y(1) + 4 * sum(y(2:2:N-1)) + 2 * sum(y(3:2:N-2)) + y(N);

+I = h* sum1 / 3;                    // Simpson's 1/3rd Integral Value

+endfunction

+

+function [f] = fun1(x)

+    f = 1 ./(1+x^2);

+endfunction

+

+

+n = 4;

+ntrap = n;

+ns13 = n/2;

+I = [trap(fun1,0,1,ntrap); simp13(fun1,0,1,ns13)];

+I = round(I*10^4)/10^4;

+true = intg(0,1,fun1);

+

+disp(I(1),"y_trap = ")

+disp(I(2),"y_simp13 = ")

+disp(I(2)*4,"Approx pi = ")

diff --git a/845/CH7/EX7.9/Ex7_9.sce b/845/CH7/EX7.9/Ex7_9.sce
new file mode 100755
index 000000000..ad4ef22ae
--- /dev/null
+++ b/845/CH7/EX7.9/Ex7_9.sce
@@ -0,0 +1,33 @@
+//Example 7.9

+

+clc

+clear

+

+function [I] = simp13 (fun,a,b,n)

+// Integrate the function over the interval using Simpson's 1/3rd rule

+// simp13 (fun,a,b,n)

+// fun - function to be integrated

+// a - lower limit of integration

+// b - upper limit of integration

+// n - No. of times simpson's 1/3rd rule needs to be performed

+

+N = 2 * n + 1;      // N - total no. of points

+h = (b-a) / (N-1);

+x = linspace(a,b,N);

+y = fun(x);

+

+sum1 = y(1) + 4 * sum(y(2:2:N-1)) + 2 * sum(y(3:2:N-2)) + y(N);

+I = h* sum1 / 3;                    // Simpson's 1/3rd Integral Value

+endfunction

+

+function [f] = fun1(x)

+    f = sqrt(2/%pi)*exp(-x^2/2);

+endfunction

+

+h = 0.125;

+n = (1-0)/h;

+ns13 = n/2;

+I = simp13(fun1,0,1,ns13);

+I = round(I*10^4)/10^4;

+

+disp(I,"Integral value, I = ")