initial commit / add all books

15 files changed, 400 insertions, 0 deletions

diff --git a/75/DEPENDENCIES/aitken1.sce b/75/DEPENDENCIES/aitken1.sce
new file mode 100755
index 000000000..3e50ad87d
--- /dev/null
+++ b/75/DEPENDENCIES/aitken1.sce
@@ -0,0 +1,10 @@
+// this program is exclusively coded to perform one iteration of aitken method,
+
+function x0aa=aitken(x0,x1,x2,g)
+x0a=x0-(x1-x0)^2/(x2-2*x1+x0);
+x1a=g(x0a);
+x2a=g(x1a);
+x0aa=x0a-(x1a-x0a)^2/(x2a-2*x1a+x0a);
+
+endfunction
+
diff --git a/75/DEPENDENCIES/bisection.sce b/75/DEPENDENCIES/bisection.sce
new file mode 100755
index 000000000..23216d2bf
--- /dev/null
+++ b/75/DEPENDENCIES/bisection.sce
@@ -0,0 +1,28 @@
+function x=bisection(a,b,f)
+    N=100;                                      //define max. number of iterations
+    PE=10^-4                                  //define tolerance
+    if (f(a)*f(b) > 0) then
+	     error('no root possible f(a)*f(b) > 0')    // checking if the decided range is containing a root
+	      abort;
+    end;
+    if(abs(f(a)) <PE) then
+	    error('solution at a')                    // seeing if there is an approximate root at a,
+	     abort;
+    end;
+    if(abs(f(b)) < PE) then                       //seeing if there is an approximate root at b,
+	error('solution at b')
+	abort;
+    end;
+    x=(a+b)/2
+    for n=1:1:N                                   //initialising 'for' loop,
+        p=f(a)*f(x)
+        if p<0 then b=x ,x=(a+x)/2;                //checking for the required conditions( f(x)*f(a)<0),
+        else
+             a=x
+            x=(x+b)/2;
+        end
+        if abs(f(x))<=PE then break                // instruction to come out of the loop after the required condition is achived,   
+        end
+    end
+    disp(n," no. of iterations =")                 // display the no. of iterations took to achive required condition,
+endfunction
\ No newline at end of file
diff --git a/75/DEPENDENCIES/bisection1.sce b/75/DEPENDENCIES/bisection1.sce
new file mode 100755
index 000000000..23216d2bf
--- /dev/null
+++ b/75/DEPENDENCIES/bisection1.sce
@@ -0,0 +1,28 @@
+function x=bisection(a,b,f)
+    N=100;                                      //define max. number of iterations
+    PE=10^-4                                  //define tolerance
+    if (f(a)*f(b) > 0) then
+	     error('no root possible f(a)*f(b) > 0')    // checking if the decided range is containing a root
+	      abort;
+    end;
+    if(abs(f(a)) <PE) then
+	    error('solution at a')                    // seeing if there is an approximate root at a,
+	     abort;
+    end;
+    if(abs(f(b)) < PE) then                       //seeing if there is an approximate root at b,
+	error('solution at b')
+	abort;
+    end;
+    x=(a+b)/2
+    for n=1:1:N                                   //initialising 'for' loop,
+        p=f(a)*f(x)
+        if p<0 then b=x ,x=(a+x)/2;                //checking for the required conditions( f(x)*f(a)<0),
+        else
+             a=x
+            x=(x+b)/2;
+        end
+        if abs(f(x))<=PE then break                // instruction to come out of the loop after the required condition is achived,   
+        end
+    end
+    disp(n," no. of iterations =")                 // display the no. of iterations took to achive required condition,
+endfunction
\ No newline at end of file
diff --git a/75/DEPENDENCIES/bvpeigen.sce b/75/DEPENDENCIES/bvpeigen.sce
new file mode 100755
index 000000000..1fcbf02aa
--- /dev/null
+++ b/75/DEPENDENCIES/bvpeigen.sce
@@ -0,0 +1,36 @@
+function [x,y,lam] = BVPeigen1(L,n)
+
+Dx = L/(n-1); 
+x=[0:Dx:L]; 
+a = 1/Dx^2; 
+k  = n-2;  
+   
+A = zeros(k,k);
+for j = 1:k
+    A(j,j) = 2*a;
+end;
+for j = 1:k-1
+    A(j,j+1) = -a;
+    A(j+1,j) = -a;
+end;
+
+exec eigenvectors.sce
+
+[yy,lam]=eigenvectors(A);
+//disp('yy');disp(yy);
+
+y = [zeros(1,k);yy;zeros(1,k)];
+//disp('y');disp(y);
+
+
+xmin=min(x);xmax=max(x);ymin=min(y);ymax=max(y);
+rect = [xmin ymin xmax ymax];
+
+if k>=5 then
+   m = 5;
+else
+   m = k;
+end
+
+
+endfunction
\ No newline at end of file
diff --git a/75/DEPENDENCIES/eigenvectors.sce b/75/DEPENDENCIES/eigenvectors.sce
new file mode 100755
index 000000000..a8ce16cf0
--- /dev/null
+++ b/75/DEPENDENCIES/eigenvectors.sce
@@ -0,0 +1,30 @@
+function [x,lam] = eigenvectors(A)
+
+//Calculates unit eigenvectors of matrix A
+//returning a matrix x whose columns are 
+//the eigenvectors.  The function also 
+//returns the eigenvalues of the matrix.
+
+[n,m] = size(A);
+
+if m<>n then
+	error('eigenvectors - matrix A is not square');
+	abort;
+end;
+
+lam = spec(A)';                   //Eigenvalues of matrix A
+
+x = [];
+
+for k = 1:n
+	B = A - lam(k)*eye(n,n); //Characteristic matrix
+        C = B(1:n-1,1:n-1);	 //Coeff. matrix for reduced system
+	b = -B(1:n-1,n);         //RHS vector for reduced system
+	y = C\b;		 //Solution for reduced system
+	y = [y;1];		 //Complete eigenvector
+	y = y/norm(y);		 //Make unit eigenvector 
+	x = [x y];		 //Add eigenvector to matrix
+end;
+
+endfunction
+//End of function
\ No newline at end of file
diff --git a/75/DEPENDENCIES/euler.sce b/75/DEPENDENCIES/euler.sce
new file mode 100755
index 000000000..7b9f184ee
--- /dev/null
+++ b/75/DEPENDENCIES/euler.sce
@@ -0,0 +1,34 @@
+function [x,y] = Euler1(x0,y0,xn,h,g)
+
+//Euler 1st order method solving ODE
+//  dy/dx = g(x,y), with initial 
+//conditions y=y0 at x = x0.  The 
+//solution is obtained for x = [x0:h:xn]
+//and returned in y
+
+ymaxAllowed = 1e+100
+
+x = [x0:h:xn];
+y = zeros(x);
+n = length(y);
+y(1) = y0;
+
+for j = 1:n-1
+	y(j+1) = y(j) + h*g(x(j),y(j));
+	if y(j+1) > ymaxAllowed then
+           disp('Euler 1 - WARNING: underflow or overflow');
+	   disp('Solution sought in the following range:');
+           disp([x0 h xn]);
+	   disp('Solution evaluated in the following range:');
+	   disp([x0 h x(j)]);
+           n = j;
+           x = x(1,1:n); y = y(1,1:n);
+	   break;
+	end;
+end;
+
+endfunction
+
+//End function Euler1
+
+
diff --git a/75/DEPENDENCIES/gausselim.sce b/75/DEPENDENCIES/gausselim.sce
new file mode 100755
index 000000000..9e1ef25f3
--- /dev/null
+++ b/75/DEPENDENCIES/gausselim.sce
@@ -0,0 +1,43 @@
+function [x] = gausselim(A,b)
+
+    //This function obtains the solution to the system of 
+    //linear equations A*x = b, given the matrix of coefficients A
+    //and the right-hand side vector, b
+
+[nA,mA] = size(A)
+[nb,mb] = size(b)
+
+if nA<>mA then
+	error('gausselim - Matrix A must be square');
+	abort;
+elseif mA<>nb then
+        error('gausselim - incompatible dimensions between A and b');
+	abort;
+end;
+
+a = [A b];
+
+    //Forward elimination
+
+n = nA; 
+for k=1:n-1
+    for i=k+1:n
+	for j=k+1:n+1
+		a(i,j)=a(i,j)-a(k,j)*a(i,k)/a(k,k);
+	end;
+    end; 
+end;
+
+    //Backward substitution
+
+x(n) = a(n,n+1)/a(n,n);
+
+for i = n-1:-1:1
+	sumk=0
+	for k=i+1:n
+		sumk=sumk+a(i,k)*x(k); 
+	end; 
+	x(i)=(a(i,n+1)-sumk)/a(i,i); 
+end;
+
+    //End function
\ No newline at end of file
diff --git a/75/DEPENDENCIES/gaussseidel.sce b/75/DEPENDENCIES/gaussseidel.sce
new file mode 100755
index 000000000..e3ae2f3e6
--- /dev/null
+++ b/75/DEPENDENCIES/gaussseidel.sce
@@ -0,0 +1,18 @@
+function [x]=gaussseidel(A,b,x0)
+[nA,mA]=size(A)
+n=nA
+[L,U] = lu(A)
+d = diag(A)
+a11 = d(1,1)
+a22 = d(2,1)
+a33 = d(3,1)
+D = [a11 0 0;0 a22 0;0 0 a33]
+H = -inv(L+D)*U
+C = inv(L+D)*b
+for m=0:3
+            x = -inv(D)*(L+U)*x + inv(D)*b
+            m=m+1
+            disp(x)
+end
+
+endfunction
\ No newline at end of file
diff --git a/75/DEPENDENCIES/lagrange.sce b/75/DEPENDENCIES/lagrange.sce
new file mode 100755
index 000000000..ec4bc3fc4
--- /dev/null
+++ b/75/DEPENDENCIES/lagrange.sce
@@ -0,0 +1,14 @@
+function[P]=lagrange(X,Y)
+
+    //    X nodes,Y values
+    //    P is the numerical Lagrange polynomial interpolation
+n=length(X)
+    //    n is the number of nodes. (n-1) is the degree
+x=poly(0,"x")
+P=0
+for i=1:n, L=1
+    for j=[1:i-1,i+1:n] L=L*(x-X(j))/(X(i)-X(j))
+        end
+P=P+L*Y(i)
+end
+endfunction
\ No newline at end of file
diff --git a/75/DEPENDENCIES/legendrepol.sce b/75/DEPENDENCIES/legendrepol.sce
new file mode 100755
index 000000000..a1b2931a6
--- /dev/null
+++ b/75/DEPENDENCIES/legendrepol.sce
@@ -0,0 +1,29 @@
+
+function [pL] = legendrepol(n,var)
+
+//    Generates the Legendre polynomial
+//    of order n in variable var
+
+if n == 0 then
+	cc = [1];
+elseif n == 1 then
+	cc = [0 1];
+else
+	if modulo(n,2) == 0 then
+		M = n/2
+	else
+		M = (n-1)/2
+	end;
+
+	cc = zeros(1,M+1);
+	for m = 0:M
+	    k = n-2*m;
+	    cc(k+1)=...
+        (-1)^m*gamma(2*n-2*m+1)/(2^n*gamma(m+1)*gamma(n-m+1)*gamma(n-2*m+1));
+	end;
+end;
+
+pL = poly(cc,var,'coeff');
+
+//    End function legendrepol
+
diff --git a/75/DEPENDENCIES/muller.sce b/75/DEPENDENCIES/muller.sce
new file mode 100755
index 000000000..4faba6b75
--- /dev/null
+++ b/75/DEPENDENCIES/muller.sce
@@ -0,0 +1,32 @@
+function x=muller(x0,x1,x2,f)
+    R=3;
+    PE=10^-8;
+    maxval=10^4;
+     for n=1:1:R
+   
+    La=(x2-x1)/(x1-x0);      
+    Da=1+La;
+    ga=La^2*f(x0)-Da^2*f(x1)+(La+Da)*f(x2);
+    Ca=La*(La*f(x0)-Da*f(x1)+f(x2));
+    
+     q=ga^2-4*Da*Ca*f(x2);
+     if q<0 then q=0;
+     end
+     p= sqrt(q);
+     if ga<0 then p=-p;
+     end
+        La=-2*Da*f(x2)/(ga+p);
+        x=x2+(x2-x1)*La; 
+        if abs(f(x))<=PE then break
+        end
+        if (abs(f(x))>maxval) then error('Solution diverges');
+            abort;
+            break
+        else 
+        x0=x1;
+        x1=x2;
+         x2=x;
+        end 
+    end
+    disp(n," no. of iterations =")
+endfunction
\ No newline at end of file
diff --git a/75/DEPENDENCIES/newton.sce b/75/DEPENDENCIES/newton.sce
new file mode 100755
index 000000000..6f7557e38
--- /dev/null
+++ b/75/DEPENDENCIES/newton.sce
@@ -0,0 +1,16 @@
+function x=newton(x,f,fp)
+    R=100;
+    PE=10^-8;
+    maxval=10^4;
+        
+    for n=1:1:R
+        x=x-f(x)/fp(x); 
+        if abs(f(x))<=PE then break
+        end
+        if (abs(f(x))>maxval) then error('Solution diverges');
+            abort
+            break
+        end
+    end
+    disp(n," no. of iterations =")
+endfunction
\ No newline at end of file
diff --git a/75/DEPENDENCIES/romberg.sce b/75/DEPENDENCIES/romberg.sce
new file mode 100755
index 000000000..fe42bb620
--- /dev/null
+++ b/75/DEPENDENCIES/romberg.sce
@@ -0,0 +1,38 @@
+function [I]=Romberg(a,b,f,h)
+
+//    This function calculates the numerical integral of f(x) between
+//    x = a and x = b, with intervals h.  Intermediate results are obtained
+//    by using SCILAB's own inttrap function
+
+x=(a:h:b)
+x1=x(1,1)
+x2=x(1,2)
+x3=x(1,3)
+x4=x(1,4)
+y1=f(x1)
+y2=f(x2)
+y3=f(x3)
+y4=f(x4)
+y=[y1 y2 y3 y4]
+I1 = inttrap(x,y)
+x=(a:h/2:b)
+x1=x(1,1)
+x2=x(1,2)
+x3=x(1,3)
+x4=x(1,4)
+x5=x(1,5)
+x6=x(1,6)
+x7=x(1,7)
+y1=f(x1)
+y2=f(x2)
+y3=f(x3)
+y4=f(x4)
+y5=f(x5)
+y6=f(x6)
+y7=f(x7)
+y=[y1 y2 y3 y4 y5 y6 y7]
+I2 = inttrap(x,y)
+I = I2 + (1.0/3.0)*(I2-I1)
+
+endfunction
+//end function Romberg
diff --git a/75/DEPENDENCIES/secant.sce b/75/DEPENDENCIES/secant.sce
new file mode 100755
index 000000000..e7d914803
--- /dev/null
+++ b/75/DEPENDENCIES/secant.sce
@@ -0,0 +1,12 @@
+function [x]=secant(a,b,f)
+    N=100;             // define max. no. iterations to be performed
+    PE=10^-4           // define tolerance for convergence
+     for n=1:1:N       // initiating for loop
+        x=a-(a-b)*f(a)/(f(a)-f(b));  
+        if abs(f(x))<=PE then break; //checking for the required condition
+        else a=b;
+            b=x;
+        end     
+     end
+     disp(n," no. of iterations =") //
+endfunction
\ No newline at end of file
diff --git a/75/DEPENDENCIES/trapezoidal.sce b/75/DEPENDENCIES/trapezoidal.sce
new file mode 100755
index 000000000..8c1182a52
--- /dev/null
+++ b/75/DEPENDENCIES/trapezoidal.sce
@@ -0,0 +1,32 @@
+function [x,y] = trapezoidal(x0,y0,xn,h,g)
+
+//Trapezoidal method solving ODE
+//  dy/dx = g(x,y), with initial 
+//conditions y=y0 at x = x0.  The 
+//solution is obtained for x = [x0:h:xn]
+//and returned in y
+
+ymaxAllowed = 1e+100
+
+x = [x0:h:xn];
+y = zeros(x);
+n = length(y);
+y(1) = y0;
+
+for j = 1:n-1
+	y(j+1) = y(j) + h*(g(x(j),y(j))+g(x(j+1),y(j+1)))/2;
+	if y(j+1) > ymaxAllowed then
+           disp('Euler 1 - WARNING: underflow or overflow');
+	   disp('Solution sought in the following range:');
+           disp([x0 h xn]);
+	   disp('Solution evaluated in the following range:');
+	   disp([x0 h x(j)]);
+           n = j;
+           x = x(1,1:n); y = y(1,1:n);
+	   break;
+	end;
+end;
+
+endfunction
+
+//End function trapezoidal
\ No newline at end of file