diff options
| author | priyanka | 2015-06-24 15:03:17 +0530 |
|---|---|---|
| committer | priyanka | 2015-06-24 15:03:17 +0530 |
| commit | b1f5c3f8d6671b4331cef1dcebdf63b7a43a3a2b (patch) | |
| tree | ab291cffc65280e58ac82470ba63fbcca7805165 /50 | |
| download | Scilab-TBC-Uploads-b1f5c3f8d6671b4331cef1dcebdf63b7a43a3a2b.tar.gz Scilab-TBC-Uploads-b1f5c3f8d6671b4331cef1dcebdf63b7a43a3a2b.tar.bz2 Scilab-TBC-Uploads-b1f5c3f8d6671b4331cef1dcebdf63b7a43a3a2b.zip | |
initial commit / add all books
Diffstat (limited to '50')
181 files changed, 3609 insertions, 0 deletions
diff --git a/50/CH2/EX2.1/2_1.PNG b/50/CH2/EX2.1/2_1.PNG Binary files differnew file mode 100755 index 000000000..a113c4de6 --- /dev/null +++ b/50/CH2/EX2.1/2_1.PNG diff --git a/50/CH2/EX2.1/ex_1.sce b/50/CH2/EX2.1/ex_1.sce new file mode 100755 index 000000000..fe2fd52ce --- /dev/null +++ b/50/CH2/EX2.1/ex_1.sce @@ -0,0 +1,18 @@ + // The equation 8*x^3-12*x^2-2*x+3==0 has three real roots.
+ // the graph of this function can be observed here.
+xset('window',0);
+x=-1:.01:2.5; // defining the range of x.
+deff('[y]=f(x)','y=8*x^3-12*x^2-2*x+3'); //defining the cunction
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+
+title(' y = 8*x^3-12*x^2-2*x+3')
+
+// from the above plot we can infre that the function has roots between
+// the intervals (-1,0),(0,1),(1,2).
diff --git a/50/CH2/EX2.11/ex_11.sce b/50/CH2/EX2.11/ex_11.sce new file mode 100755 index 000000000..92b538cc7 --- /dev/null +++ b/50/CH2/EX2.11/ex_11.sce @@ -0,0 +1,25 @@ + // The equation x^3-5*x+1==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',10);
+x=-2:.01:4; // defining the range of x.
+deff('[y]=f(x)','y=x^3-5*x+1'); //defining the cunction.
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = x^3-5*x+1')
+
+// from the above plot we can infre that the function has roots between
+// the intervals (0,1),(2,3).
+// since we have been asked for the smallest positive root of the equation,
+// we are intrested on the interval (0,1)
+
+
+ //sollution by muller method to 3 iterations.
+
+muller3(0,.5,1,f)
+
diff --git a/50/CH2/EX2.12/2_12.PNG b/50/CH2/EX2.12/2_12.PNG Binary files differnew file mode 100755 index 000000000..4c46e2515 --- /dev/null +++ b/50/CH2/EX2.12/2_12.PNG diff --git a/50/CH2/EX2.12/ex_12.sce b/50/CH2/EX2.12/ex_12.sce new file mode 100755 index 000000000..42225d058 --- /dev/null +++ b/50/CH2/EX2.12/ex_12.sce @@ -0,0 +1,24 @@ + // The equation cos(x)-x*%e^x==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',8);
+x=-1:.001:2; // defining the range of x.
+deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the cunction.
+deff('[y]=fp(x)','y=-sin(x)-x*%e^x-%e^x');
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = cos(x)-x*%e^x')
+
+// from the above plot we can infre that the function has root between
+// the interval (0,1)
+
+
+ //sollution by muller method to 5 iterations.
+
+
+muller5(-1,0,1,f)
diff --git a/50/CH2/EX2.13/2_13.PNG b/50/CH2/EX2.13/2_13.PNG Binary files differnew file mode 100755 index 000000000..0dca7c750 --- /dev/null +++ b/50/CH2/EX2.13/2_13.PNG diff --git a/50/CH2/EX2.13/ex_13.sce b/50/CH2/EX2.13/ex_13.sce new file mode 100755 index 000000000..1a87dd5e8 --- /dev/null +++ b/50/CH2/EX2.13/ex_13.sce @@ -0,0 +1,33 @@ + // The equation x^3-5*x+1==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',12);
+x=-2:.01:4; // defining the range of x.
+deff('[y]=f(x)','y=x^3-5*x+1'); //defining the function.
+deff('[y]=fp(x)','y=3*x^2-5');
+deff('[y]=fpp(x)','y=6*x');
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = x^3-5*x+1')
+
+// from the above plot we can infre that the function has roots between
+// the intervals (0,1),(2,3).
+// since we have been asked for the smallest positive root of the equation,
+// we are intrested on the interval (0,1)
+// a=0;b=1,
+
+
+// solution by chebyshev method
+
+// the approximate root after 4 iterations can be observed.
+
+
+chebyshev(0.5,f,fp)
+
+
+// hence the approximate root witin the permissible error of 10^-15 is .2016402,
\ No newline at end of file diff --git a/50/CH2/EX2.14/2_14.PNG b/50/CH2/EX2.14/2_14.PNG Binary files differnew file mode 100755 index 000000000..884daebb0 --- /dev/null +++ b/50/CH2/EX2.14/2_14.PNG diff --git a/50/CH2/EX2.14/ex_14.sce b/50/CH2/EX2.14/ex_14.sce new file mode 100755 index 000000000..54e521762 --- /dev/null +++ b/50/CH2/EX2.14/ex_14.sce @@ -0,0 +1,30 @@ +
+
+ // The equation 1/x-7==0 has a real root.
+ // the graph of this function can be observed here.
+xset('window',13);
+x=0.001:.001:.25; // defining the range of x.
+deff('[y]=f(x)','y=1/x-7'); //defining the function.
+deff('[y]=fp(x)','y=-1/x^2');
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y =1/x-7')
+
+// from the above plot we can infre that the function has roots between
+// the interval (0,2/7)
+
+
+ //solution by chebyshev method
+
+
+ chebyshev(0.1,f,fp) //calling the pre-defined function 'chebyshev' to find the approximate root in the range of (0,2/7).
+
+
+
+
\ No newline at end of file diff --git a/50/CH2/EX2.15/2_15.PNG b/50/CH2/EX2.15/2_15.PNG Binary files differnew file mode 100755 index 000000000..5c4de8ecc --- /dev/null +++ b/50/CH2/EX2.15/2_15.PNG diff --git a/50/CH2/EX2.15/ex_15.sce b/50/CH2/EX2.15/ex_15.sce new file mode 100755 index 000000000..b8bfdab9a --- /dev/null +++ b/50/CH2/EX2.15/ex_15.sce @@ -0,0 +1,35 @@ + // The equation cos(x)-x*%e^x==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',8);
+x=-1:.001:2; // defining the range of x.
+deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the cunction.
+deff('[y]=fp(x)','y=-sin(x)-x*%e^x-%e^x');
+deff('[y]=fpp(x)','y=-cos(x)-x*%e^x-2*%e^x');
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = cos(x)-x*%e^x')
+
+// from the above plot we can infre that the function has root between
+// the interval (0,1)
+
+
+// a=0;b=1,
+
+
+
+ // solution by chebyshev with a permissible error of 10^-15.
+
+// we call a user-defined function 'chebyshev' so as to find the approximate
+// root of the equation within the defined permissible error limit.
+
+chebyshev(1,f,fp)
+
+
+
+// hence the approximate root witin the permissible error of 10^-15 is
\ No newline at end of file diff --git a/50/CH2/EX2.16/2_16.PNG b/50/CH2/EX2.16/2_16.PNG Binary files differnew file mode 100755 index 000000000..dd5cb960c --- /dev/null +++ b/50/CH2/EX2.16/2_16.PNG diff --git a/50/CH2/EX2.16/ex_16.sce b/50/CH2/EX2.16/ex_16.sce new file mode 100755 index 000000000..d1a65a8ca --- /dev/null +++ b/50/CH2/EX2.16/ex_16.sce @@ -0,0 +1,38 @@ + // The equation x^3-5*x+1==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',15);
+x=-2:.01:4; // defining the range of x.
+deff('[y]=f(x)','y=x^3-5*x+1'); //defining the function.
+deff('[y]=fp(x)','y=3*x^2-5');
+deff('[y]=fpp(x)','y=6*x');
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = x^3-5*x+1')
+
+// from the above plot we can infre that the function has roots between
+// the intervals (0,1),(2,3).
+// since we have been asked for the smallest positive root of the equation,
+// we are intrested on the interval (0,1)
+// a=0;b=1,
+
+
+// solution by multipoint iteration method
+
+// the approximate root after 3 iterations can be observed.
+
+
+multipoint_iteration31(0.5,f,fp)
+
+// hence the approximate root witin the permissible error of 10^-15 is .201640,
+
+
+
+multipoint_iteration33(0.5,f,fp)
+
+// hence the approximate root witin the permissible error of 10^-15 is .201640,
\ No newline at end of file diff --git a/50/CH2/EX2.17/ex_17.sce b/50/CH2/EX2.17/ex_17.sce new file mode 100755 index 000000000..f818a0af8 --- /dev/null +++ b/50/CH2/EX2.17/ex_17.sce @@ -0,0 +1,34 @@ + // The equation cos(x)-x*%e^x==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',8);
+x=-1:.001:2; // defining the range of x.
+deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the function.
+deff('[y]=fp(x)','y=-sin(x)-x*%e^x-%e^x');
+deff('[y]=fpp(x)','y=-cos(x)-x*%e^x-2*%e^x');
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = cos(x)-x*%e^x')
+
+// from the above plot we can infre that the function has root between
+// the interval (0,1)
+
+
+// a=0;b=1,
+
+
+
+ // solution by multipoint_iteration method using the formula given in equation no.2.33.
+
+// we call a user-defined function 'multipoint_iteration33' so as to find the approximate
+// root of the equation within the defined permissible error limit.
+
+multipoint_iteration33(1,f,fp)
+
+
+// hence the approximate root witin the permissible error of 10^-5 is 0.5177574.
\ No newline at end of file diff --git a/50/CH2/EX2.2/2_2.PNG b/50/CH2/EX2.2/2_2.PNG Binary files differnew file mode 100755 index 000000000..c38387579 --- /dev/null +++ b/50/CH2/EX2.2/2_2.PNG diff --git a/50/CH2/EX2.2/ex_2.sce b/50/CH2/EX2.2/ex_2.sce new file mode 100755 index 000000000..3fa2a9a1b --- /dev/null +++ b/50/CH2/EX2.2/ex_2.sce @@ -0,0 +1,17 @@ + // The equation cos(x)-x*%e^x==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',1);
+x=0:.01:2; // defining the range of x.
+deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the cunction.
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = cos(x)-x*%e^x')
+
+// from the above plot we can infre that the function has root between
+// the interval (0,1)
diff --git a/50/CH2/EX2.23/ex_23.sce b/50/CH2/EX2.23/ex_23.sce new file mode 100755 index 000000000..3a72496d0 --- /dev/null +++ b/50/CH2/EX2.23/ex_23.sce @@ -0,0 +1,48 @@ + // The equation 3*x^3+4*x^2+4*x+1==0 has three real roots. + // the graph of this function can be observed here. +xset('window',22); +x=-1.5:.001:1.5; // defining the range of x. +deff('[y]=f(x)','y=3*x^3+4*x^2+4*x+1'); //defining the cunction +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph + +title(' y =3*x^3+4*x^2+4*x+1') + +// from the above plot we can infre that the function has root between +// the interval (-1,0), + +x0=-.5; // initial approximation + + +// let the iterative function g(x) be x+A*(3*x^3+4*x^2+4*x+1) =g(x); + +// gp(x)=(1+A*(9*x^2+8*x+4) ) +// we need to choose a value for A , which makes abs(gp(x0))<1 + +// hence abs(gp(x0))=abs(1+9*A/4) + +A=-1:.1:1; + +abs(1+9*A/4) // tryin to check the values of abs(gp(x0)) for different values of A. + + +// from the above values of 'A' and the values of 'abs(gp(x0))', +// we can infer that for the vales of 'A 'in the range (-.8,0) g(x ) will be giving a converging solution , + +// hence deliberatele we choose a to be -0.5, + +A=-0.5; + +deff('[y]=g(x)','y= x-0.5*(3*x^3+4*x^2+4*x+1)'); +deff('[y]=gp(x)','y= 1-0.5*(9*x^2+8*x+4)'); // hence defining g(x) and gp(x), +generaliteration(x0,g,gp) + + + + diff --git a/50/CH2/EX2.24.1/2_24_1.PNG b/50/CH2/EX2.24.1/2_24_1.PNG Binary files differnew file mode 100755 index 000000000..8a9e4a241 --- /dev/null +++ b/50/CH2/EX2.24.1/2_24_1.PNG diff --git a/50/CH2/EX2.24.1/ex_24_1.sce b/50/CH2/EX2.24.1/ex_24_1.sce new file mode 100755 index 000000000..5745698d7 --- /dev/null +++ b/50/CH2/EX2.24.1/ex_24_1.sce @@ -0,0 +1,39 @@ + // The equation x^3-5*x+1==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',2);
+x=-2:.01:4; // defining the range of x.
+deff('[y]=f(x)','y=x^3-5*x+1'); //defining the function.
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = x^3-5*x+1')
+
+// from the above plot we can infre that the function has roots between
+// the intervals (0,1),(2,3).
+// since we have been asked for the smallest positive root of the equation,
+// we are intrested on the interval (0,1)
+
+x0=.5;
+
+ //solution using linear iteration method for the first two iterations and aitken's process two times for the third iteration.
+
+ deff('[y]=g(x)','y=1/5*(x^3+1)');
+ deff('[y]=gp(x)','y=1/5*(3*x^2)');
+
+
+generaliteration2(x0,g,gp)
+
+
+// from the above iterations performed we can infer that-
+x1=0.225;
+x2=0.202278;
+
+
+
+
+aitken(x0,x1,x2,g) // calling the aitken method for one iteration
\ No newline at end of file diff --git a/50/CH2/EX2.24.2/2_24_2.PNG b/50/CH2/EX2.24.2/2_24_2.PNG Binary files differnew file mode 100755 index 000000000..1b89face5 --- /dev/null +++ b/50/CH2/EX2.24.2/2_24_2.PNG diff --git a/50/CH2/EX2.24.2/ex_24_2.sce b/50/CH2/EX2.24.2/ex_24_2.sce new file mode 100755 index 000000000..9ff6fe046 --- /dev/null +++ b/50/CH2/EX2.24.2/ex_24_2.sce @@ -0,0 +1,39 @@ + // The equation x-%e^-x==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',24);
+x=-3:.01:4; // defining the range of x.
+deff('[y]=f(x)','y=x-%e^-x'); //defining the cunction.
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = x-%e^-x')
+
+// from the above plot we can infre that the function has root between
+// the interval (0,1)
+
+x0=1;
+
+ //solution using linear iteration method for the first two iterations and aitken's process two times for the third iteration.
+
+
+
+ deff('[y]=g(x)','y=%e^-x');
+ deff('[y]=gp(x)','y=-%e^-x');
+
+
+generaliteration2(x0,g,gp)
+
+
+// from the above iterations performed we can infer that-
+x1=0.367879;
+x2=0.692201;
+
+
+
+
+aitken(x0,x1,x2,g) // calling the aitken method for one iteration
\ No newline at end of file diff --git a/50/CH2/EX2.25/2_25.PNG b/50/CH2/EX2.25/2_25.PNG Binary files differnew file mode 100755 index 000000000..1478a1d27 --- /dev/null +++ b/50/CH2/EX2.25/2_25.PNG diff --git a/50/CH2/EX2.25/ex_25.sce b/50/CH2/EX2.25/ex_25.sce new file mode 100755 index 000000000..a9c110187 --- /dev/null +++ b/50/CH2/EX2.25/ex_25.sce @@ -0,0 +1,36 @@ + // The equation cos(x)-x*%e^x==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',25);
+x=0:.01:2; // defining the range of x.
+deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the cunction.
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = cos(x)-x*%e^x')
+
+// from the above plot we can infre that the function has root between
+// the interval (0,1)
+
+x0=0;
+
+ //solution using linear iteration method for the first two iterations and aitken's process two times for the third iteration.
+
+ deff('[y]=g(x)','y=x+1/2*(cos(x)-x*%e^x)');
+ deff('[y]=gp(x)','y=1+1/2*(-sin(x)-x*%e^x-%e^x)');
+
+
+generaliteration2(x0,g,gp)
+
+
+// from the above iterations performed we can infer that-
+x1=0.50000000;
+x2=0.5266110;
+
+
+
+aitken(x0,x1,x2,g) // calling the aitken method for one iteration
diff --git a/50/CH2/EX2.26/2_26.PNG b/50/CH2/EX2.26/2_26.PNG Binary files differnew file mode 100755 index 000000000..d303cd438 --- /dev/null +++ b/50/CH2/EX2.26/2_26.PNG diff --git a/50/CH2/EX2.26/ex_26.sce b/50/CH2/EX2.26/ex_26.sce new file mode 100755 index 000000000..2e5761bf3 --- /dev/null +++ b/50/CH2/EX2.26/ex_26.sce @@ -0,0 +1,36 @@ + // The equation x^3-7*x^2+16*x-12==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',25);
+x=0:.001:4; // defining the range of x.
+deff('[y]=f(x)','y=x^3-7*x^2+16*x-12'); //defining the cunction.
+deff('[y]=fp(x)','y=3*x^2-14*x+16');
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = x^3-7*x^2+16*x-12')
+
+
+
+
+// given that the equation has double roots at x=2 hence m=2;
+
+m=2;
+
+ // solution by newton raphson method
+
+
+newton(1,f,fp) // calling the user defined function
+
+
+
+
+ //solution by modified newton raphsons mathod
+
+
+
+modified_newton(1,f,fp)
\ No newline at end of file diff --git a/50/CH2/EX2.27/2_27.PNG b/50/CH2/EX2.27/2_27.PNG Binary files differnew file mode 100755 index 000000000..8cef9c44f --- /dev/null +++ b/50/CH2/EX2.27/2_27.PNG diff --git a/50/CH2/EX2.27/ex_27.sce b/50/CH2/EX2.27/ex_27.sce new file mode 100755 index 000000000..2d97164c5 --- /dev/null +++ b/50/CH2/EX2.27/ex_27.sce @@ -0,0 +1,46 @@ + // The equation 27*x^5+27*x^4+36*x^3+28*x^2+9*x+1==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',26);
+x=-2:.001:3; // defining the range of x.
+deff('[y]=f(x)','y=27*x^5+27*x^4+36*x^3+28*x^2+9*x+1'); //defining the cunction.
+deff('[y]=fp(x)','y=27*5*x^4+27*4*x^3+36*3*x^2+28*2*x+9');
+deff('[y]=fpp(x)','y=27*5*4*x^3+27*4*3*x^2+36*3*2*x+28*2');
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = 27*x^5+27*x^4+36*x^3+28*x^2+9*x+1')
+
+
+
+
+ // solution by newton raphson method as per the equation no. 2.14
+
+
+newton(-1,f,fp) // calling the user defined function
+
+
+newton4(-1,f,fp)
+
+
+ // solution by newton raphson method as per the equation no. 2.63
+
+newton63(-1,f,fp,fpp) // calling the user defined function
+ // solution by the secant method defined to satisfy the equation no.2.64.
+
+
+secant64(0,-1,f,fp)
+
+
+
+
+
+
+ // solution by the secant method defined to satisfy the equation no.2.65.
+
+
+secant65(0,-.5,f)
diff --git a/50/CH2/EX2.3/2_3.PNG b/50/CH2/EX2.3/2_3.PNG Binary files differnew file mode 100755 index 000000000..1d812544f --- /dev/null +++ b/50/CH2/EX2.3/2_3.PNG diff --git a/50/CH2/EX2.3/ex_3.sce b/50/CH2/EX2.3/ex_3.sce new file mode 100755 index 000000000..7add08896 --- /dev/null +++ b/50/CH2/EX2.3/ex_3.sce @@ -0,0 +1,35 @@ + // The equation x^3-5*x+1==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',2);
+x=-2:.01:4; // defining the range of x.
+deff('[y]=f(x)','y=x^3-5*x+1'); //defining the cunction.
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = x^3-5*x+1')
+
+// from the above plot we can infre that the function has roots between
+// the intervals (0,1),(2,3).
+// since we have been asked for the smallest positive root of the equation,
+// we are intrested on the interval (0,1)
+// a=0;b=1,
+
+// we call a user-defined function 'bisection' so as to find the approximate
+// root of the equation with a defined permissible error.
+
+bisection(0,1,f)
+
+// since in the example 2.3 we have been asked to perform 5 itterations
+// the approximate root after 5 iterations can be observed.
+
+
+
+bisection5(0,1,f)
+
+
+// hence the approximate root after 5 iterations is 0.203125 witin the permissible error of 10^-4,
\ No newline at end of file diff --git a/50/CH2/EX2.4/2_4.PNG b/50/CH2/EX2.4/2_4.PNG Binary files differnew file mode 100755 index 000000000..e111288e6 --- /dev/null +++ b/50/CH2/EX2.4/2_4.PNG diff --git a/50/CH2/EX2.4/ex_4.sce b/50/CH2/EX2.4/ex_4.sce new file mode 100755 index 000000000..751cc0f5c --- /dev/null +++ b/50/CH2/EX2.4/ex_4.sce @@ -0,0 +1,32 @@ + // The equation cos(x)-x*%e^x==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',3);
+x=0:.01:2; // defining the range of x.
+deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the cunction.
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = cos(x)-x*%e^x')
+
+// from the above plot we can infre that the function has root between
+// the interval (0,1)
+
+
+// a=0;b=1,
+
+// we call a user-defined function 'bisection' so as to find the approximate
+// root of the equation with a defined permissible error.
+
+bisection(0,1,f)
+
+// since in the example 2.4 we have been asked to perform 5 itterations ,
+
+bisection5(0,1,f)
+
+
+// hence the approximate root after 5 iterations is 0.515625 witin the permissible error of 10^-4,
\ No newline at end of file diff --git a/50/CH2/EX2.5/ex_5.sce b/50/CH2/EX2.5/ex_5.sce new file mode 100755 index 000000000..0a6ff1b8e --- /dev/null +++ b/50/CH2/EX2.5/ex_5.sce @@ -0,0 +1,48 @@ + // The equation x^3-5*x+1==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',4);
+x=-2:.01:4; // defining the range of x.
+deff('[y]=f(x)','y=x^3-5*x+1'); //defining the cunction.
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = x^3-5*x+1')
+
+// from the above plot we can infre that the function has roots between
+// the intervals (0,1),(2,3).
+// since we have been given the interval to be considered as (0,1)
+// a=0;b=1,
+
+
+ // Solution by secant method
+
+
+
+
+
+// since in the example 2.5 we have been asked to perform 4 itterations ,
+secant4(0,1,f) // we call a user-defined function 'bisection' so as to find the approximate
+// root of the equation with a defined permissible error.
+
+
+
+// hence the approximate root occured in secant method after 4 iterations is 0.201640 witin the permissible error of 10^-4,
+
+
+
+ // solution by regular falsi method
+
+
+// since in the example 2.5 we have been asked to perform 4 itterations ,
+
+regulafalsi4(0,1,f) // we call a user-defined function 'regularfalsi4' so as to find the approximate
+// root of the equation with a defined permissible error.
+
+
+
+// hence the approximate root occured in regularfalsi method after 4 iterations is 0.201640 witin the permissible error of 10^-4,
\ No newline at end of file diff --git a/50/CH2/EX2.6/2_6.PNG b/50/CH2/EX2.6/2_6.PNG Binary files differnew file mode 100755 index 000000000..3fafbbe75 --- /dev/null +++ b/50/CH2/EX2.6/2_6.PNG diff --git a/50/CH2/EX2.6/ex_6.sce b/50/CH2/EX2.6/ex_6.sce new file mode 100755 index 000000000..0aa95ce3a --- /dev/null +++ b/50/CH2/EX2.6/ex_6.sce @@ -0,0 +1,55 @@ + // The equation cos(x)-x*%e^x==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',3);
+x=0:.01:2; // defining the range of x.
+deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the cunction.
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = cos(x)-x*%e^x')
+
+// from the above plot we can infre that the function has root between
+// the interval (0,1)
+
+
+// a=0;b=1,
+
+
+ // Solution by secant method
+
+
+
+
+
+// since in the example 2.6 we have no specification of the no. of itterations ,
+// we define a function 'secant' and execute it.
+
+
+
+secant(0,1,f) // we call a user-defined function 'secant' so as to find the approximate
+// root of the equation with a defined permissible error.
+
+
+
+// hence the approximate root occured in secant method witin the permissible error of 10^-5 is ,
+
+
+
+ // solution by regular falsi method
+
+
+
+
+
+// since in the example 2.6 we have no specification of the no. of itterations ,
+
+
+regulafalsi(0,1,f) // we call a user-defined function 'regularfalsi' so as to find the approximate
+// root of the equation with a defined permissible error.
+
+
diff --git a/50/CH2/EX2.7/2_7.PNG b/50/CH2/EX2.7/2_7.PNG Binary files differnew file mode 100755 index 000000000..dcb0e7198 --- /dev/null +++ b/50/CH2/EX2.7/2_7.PNG diff --git a/50/CH2/EX2.7/ex_7.sce b/50/CH2/EX2.7/ex_7.sce new file mode 100755 index 000000000..ad7b083d7 --- /dev/null +++ b/50/CH2/EX2.7/ex_7.sce @@ -0,0 +1,30 @@ + // The equation x^3-5*x+1==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',6);
+x=-2:.01:4; // defining the range of x.
+deff('[y]=f(x)','y=x^3-5*x+1'); //defining the function.
+deff('[y]=fp(x)','y=3*x^2-5');
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = x^3-5*x+1')
+
+// from the above plot we can infre that the function has roots between
+// the intervals (0,1),(2,3).
+// since we have been asked for the smallest positive root of the equation,
+// we are intrested on the interval (0,1)
+// a=0;b=1,
+
+// since in the example 2.7 we have been asked to perform 4 itterations ,
+// the approximate root after 4 iterations can be observed.
+
+
+newton4(0.5,f,fp)
+
+
+// hence the approximate root after 4 iterations is 0.201640 witin the permissible error of 10^-15,
\ No newline at end of file diff --git a/50/CH2/EX2.8/2_8.PNG b/50/CH2/EX2.8/2_8.PNG Binary files differnew file mode 100755 index 000000000..bb4e1f32b --- /dev/null +++ b/50/CH2/EX2.8/2_8.PNG diff --git a/50/CH2/EX2.8/ex_8.sce b/50/CH2/EX2.8/ex_8.sce new file mode 100755 index 000000000..44ac39cef --- /dev/null +++ b/50/CH2/EX2.8/ex_8.sce @@ -0,0 +1,31 @@ + // The equation x^3-17==0 has three real roots.
+ // the graph of this function can be observed here.
+xset('window',7);
+x=-5:.001:5; // defining the range of x.
+deff('[y]=f(x)','y=x^3-17'); //defining the cunction.
+deff('[y]=fp(x)','y=3*x^2');
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = x^3-17')
+
+// from the above plot we can infre that the function has root between
+// the interval (2,3).
+
+
+
+ //solution by newton raphson's method
+
+
+
+ // since in example no.2.8 we have been asked to perform 4 iterations ,we define a fuction newton4'' which does newton raphson's method of finding approximate root upto 4 iterations,
+
+
+
+ newton4(2,f,fp) //calling the pre-defined function 'newton4'.
+
\ No newline at end of file diff --git a/50/CH2/EX2.9/2_9.PNG b/50/CH2/EX2.9/2_9.PNG Binary files differnew file mode 100755 index 000000000..a172ae0c6 --- /dev/null +++ b/50/CH2/EX2.9/2_9.PNG diff --git a/50/CH2/EX2.9/ex_9.sce b/50/CH2/EX2.9/ex_9.sce new file mode 100755 index 000000000..5fdd426b7 --- /dev/null +++ b/50/CH2/EX2.9/ex_9.sce @@ -0,0 +1,37 @@ + // The equation cos(x)-x*%e^x==0 has real roots.
+ // the graph of this function can be observed here.
+xset('window',8);
+x=-1:.001:2; // defining the range of x.
+deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the cunction.
+deff('[y]=fp(x)','y=-sin(x)-x*%e^x-%e^x');
+y=feval(x,f);
+
+a=gca();
+
+a.y_location = "origin";
+
+a.x_location = "origin";
+plot(x,y) // instruction to plot the graph
+title(' y = cos(x)-x*%e^x')
+
+// from the above plot we can infre that the function has root between
+// the interval (0,1)
+
+
+// a=0;b=1,
+
+
+
+ // solution by newton raphson's method with a permissible error of 10^-8.
+
+
+// we call a user-defined function 'newton' so as to find the approximate
+// root of the equation within the defined permissible error limit.
+
+newton(1,f,fp)
+
+
+
+
+
+// hence the approximate root witin the permissible error of 10^-8 is 0.5177574.
\ No newline at end of file diff --git a/50/CH3/EX3.1/ex_3_1.sce b/50/CH3/EX3.1/ex_3_1.sce new file mode 100755 index 000000000..473e588b6 --- /dev/null +++ b/50/CH3/EX3.1/ex_3_1.sce @@ -0,0 +1,14 @@ +// example 3.1
+//Positive definite
+
+A=[12 4 -1;4 7 1;-1 1 6];
+//check if the determinant of the leading minors is a positive value-
+
+A1=A(1);
+det(A1)
+A2=A(1:2,1:2);
+det(A2)
+A3=A(1:3,1:3);
+det(A3)
+
+//we observe that the determinant of the leading minors is a positive ,hence the given matrix A is a positive definite.
\ No newline at end of file diff --git a/50/CH3/EX3.10/ex_3_10.sce b/50/CH3/EX3.10/ex_3_10.sce new file mode 100755 index 000000000..05f721b06 --- /dev/null +++ b/50/CH3/EX3.10/ex_3_10.sce @@ -0,0 +1,13 @@ +//example no. 3.10
+//solve system by decomposition method
+
+A=[1 1 1;4 3 -1;3 5 3]
+n=3;
+
+b=[1;6;4]
+
+[U,L]=LandU(A,3)
+
+Z=fore(L,b)
+
+X=back(U,Z)
\ No newline at end of file diff --git a/50/CH3/EX3.11/ex_3_11.sce b/50/CH3/EX3.11/ex_3_11.sce new file mode 100755 index 000000000..fdf0b3b7c --- /dev/null +++ b/50/CH3/EX3.11/ex_3_11.sce @@ -0,0 +1,12 @@ +//example no.3.11
+//caption: Inverse using LU decomposition
+
+A=[3 2 1;2 3 2;1 2 2]
+
+[U,L]=LandU(A,3) // call LandU function to evaluate U,L of A,
+
+//since A=L*U ,
+// inv(A)=inv(U)*inv(L)
+// let inv(A)=AI
+
+AI=U^-1*L^-1
\ No newline at end of file diff --git a/50/CH3/EX3.12/ex_3_12.sce b/50/CH3/EX3.12/ex_3_12.sce new file mode 100755 index 000000000..4599b1c23 --- /dev/null +++ b/50/CH3/EX3.12/ex_3_12.sce @@ -0,0 +1,28 @@ +//example no. 3.12
+//solve system by decomposition method
+
+A=[1 1 -1;2 2 5;3 2 -3]
+
+b=[2;-3;6]
+
+
+
+ // hence we can observe that LU decomposition method fails to solve this system since the pivot L(2,2)=0;
+
+
+ //we note that the coefficient matrix is not a positive definite matrix and hence its LU decomposition is not guaranteed,
+
+
+ //if we interchange the rows of A as shown below the LU decomposition would work,
+
+ A=[3 2 -3;2 2 5;1 1 -1]
+
+ b=[6;-3;2]
+
+ [U,L]=LandU(A,3) // call LandU function to evaluate U,L of A,
+
+n=3;
+Z=fore(L,b);
+
+X=back(U,Z)
+
\ No newline at end of file diff --git a/50/CH3/EX3.13/ex_3_13.sce b/50/CH3/EX3.13/ex_3_13.sce new file mode 100755 index 000000000..319071f40 --- /dev/null +++ b/50/CH3/EX3.13/ex_3_13.sce @@ -0,0 +1,18 @@ +//example no. 3.13
+//solve system by LU decomposition method
+
+A=[2 1 1 -2;4 0 2 1;3 2 2 0;1 3 2 -1]
+
+b=[-10;8;7;-5]
+
+[U,L]=LandU(A,4)
+n=4;
+Z=fore(L,b);
+
+X=back(U,Z)
+
+//since A=L*U ,
+// inv(A)=inv(U)*inv(L)
+// let inv(A)=AI
+
+AI=U^-1*L^-1
diff --git a/50/CH3/EX3.14/ex_3_14.sce b/50/CH3/EX3.14/ex_3_14.sce new file mode 100755 index 000000000..3590d9c29 --- /dev/null +++ b/50/CH3/EX3.14/ex_3_14.sce @@ -0,0 +1,12 @@ +//example no. 3.14
+//solve system by cholesky method
+
+A=[1 2 3;2 8 22;3 22 82]
+
+b=[5;6;-10]
+
+L=cholesky (A,3) //call cholesky function to evaluate the root of the system
+n=3;
+Z=fore(L,b);
+
+X=back(L',Z)
\ No newline at end of file diff --git a/50/CH3/EX3.15/ex_3_15.sce b/50/CH3/EX3.15/ex_3_15.sce new file mode 100755 index 000000000..902ce414b --- /dev/null +++ b/50/CH3/EX3.15/ex_3_15.sce @@ -0,0 +1,19 @@ +//example no. 3.15
+//solve system by cholesky method
+
+A=[4 -1 0 0;-1 4 -1 0;0 -1 4 -1;0 0 -1 4]
+
+b=[1;0;0;0]
+
+L=cholesky (A,4) //call cholesky function to evaluate the root of the system
+
+n=4;
+Z=fore(L,b);
+
+X=back(L',Z)
+
+//since A=L*L' ,
+// inv(A)=inv(L')*inv(L)
+// let inv(A)=AI
+
+AI=L'^-1*L^-1
\ No newline at end of file diff --git a/50/CH3/EX3.2/ex_3_2.sce b/50/CH3/EX3.2/ex_3_2.sce new file mode 100755 index 000000000..bc98fc202 --- /dev/null +++ b/50/CH3/EX3.2/ex_3_2.sce @@ -0,0 +1,13 @@ +//example no.3.2 +// show if a given matrix 'A' possesses 'property A' + +A=[2 -1 0 -1;-1 2 -1 0;0 -1 2 0;0 0 -1 2 ] + +P=[0 0 0 1;0 1 0 0;0 0 1 0;1 0 0 0] // let us take the pemutation matrix as P + +// then we find that + +P*A*P' + +// is of the form (3.2) . hence the matrix 'A' has property A + diff --git a/50/CH3/EX3.21/ex_3_21.sce b/50/CH3/EX3.21/ex_3_21.sce new file mode 100755 index 000000000..890f06c4d --- /dev/null +++ b/50/CH3/EX3.21/ex_3_21.sce @@ -0,0 +1,14 @@ +//example no. 3.21
+//solve the system by jacobi iteration method
+
+A=[4 1 1;1 5 2;1 2 3]
+
+b=[2 ;-6;-4]
+
+N=3; //no. of ierations
+n=3; // order of the matrix is n*n
+
+X=[.5;-.5;-.5] //initial approximation
+
+
+jacobiiteration(A,n,N,X,b) //call the function which performs jacobi iteration method to solve the system
\ No newline at end of file diff --git a/50/CH3/EX3.22/ex_3_22.sce b/50/CH3/EX3.22/ex_3_22.sce new file mode 100755 index 000000000..4f6e78114 --- /dev/null +++ b/50/CH3/EX3.22/ex_3_22.sce @@ -0,0 +1,14 @@ +//example no. 3.22
+//solve the system by gauss seidel method
+
+A=[2 -1 0;-1 2 -1;0 -1 2]
+
+b=[7;1;1]
+
+N=3; //no. of ierations
+n=3; // order of the matrix is n*n
+
+X=[0;0;0] //initial approximation
+
+
+gaussseidel(A,n,N,X,b) //call the function which performs gauss seidel method to solve the system
\ No newline at end of file diff --git a/50/CH3/EX3.27/ex_3_27.sce b/50/CH3/EX3.27/ex_3_27.sce new file mode 100755 index 000000000..242a16db7 --- /dev/null +++ b/50/CH3/EX3.27/ex_3_27.sce @@ -0,0 +1,24 @@ +// example 3.27
+// a) find eigenvalue and eigen vector;
+// b) verify inv(S)*A*S is a diagonal matrix;
+
+// 1)
+A=[1 2 -2 ;1 1 1;1 3 -1];
+
+B=[1 0 0;0 1 0; 0 0 1];
+
+ [x,lam] = geigenvectors(A,B);
+
+ inv(x)*A*x
+
+ // 2)
+A=[3 2 2;2 5 2;2 2 3];
+
+B=[1 0 0;0 1 0; 0 0 1];
+
+
+ [x,lam] = geigenvectors(A,B);
+
+ inv(x)*A*x
+
+
\ No newline at end of file diff --git a/50/CH3/EX3.4/ex_3_4.sce b/50/CH3/EX3.4/ex_3_4.sce new file mode 100755 index 000000000..9cbb29d89 --- /dev/null +++ b/50/CH3/EX3.4/ex_3_4.sce @@ -0,0 +1,30 @@ +// example no.3.4
+// solve the system of equations
+
+//(a) . by using cramer's rule,
+
+A=[1 2 -1;3 6 1;3 3 2]
+
+B1=[2 2 -1;1 6 1;3 3 2]
+
+B2=[1 2 -1;3 1 1;3 3 2]
+
+B3=[1 2 2;3 6 1; 3 3 3]
+
+
+//we know;
+
+X1=det(B1)/det(A)
+X2=det(B2)/det(A)
+X3=det(B3)/det(A)
+
+
+//(b). by detemining the inverse of the coefficient matrix
+
+A=[1 2 -1;3 6 1;3 3 2]
+
+b=[2 ;1 ;3]
+
+//we know;
+
+X=inv(A)*b
\ No newline at end of file diff --git a/50/CH3/EX3.5/ex_3_5.sce b/50/CH3/EX3.5/ex_3_5.sce new file mode 100755 index 000000000..7f23bc40e --- /dev/null +++ b/50/CH3/EX3.5/ex_3_5.sce @@ -0,0 +1,9 @@ +//example-3.5
+//caption-solution by gauss elimibnation method
+
+A=[10 -1 2;1 10 -1;2 3 20] //matrices A and b from the above
+ //system of equations
+b=[4;3;7]
+
+gausselim(A,b) //call gauss elimination function to solve the
+ //matrices A and b
\ No newline at end of file diff --git a/50/CH3/EX3.6/ex_3_6.sce b/50/CH3/EX3.6/ex_3_6.sce new file mode 100755 index 000000000..47ca03a96 --- /dev/null +++ b/50/CH3/EX3.6/ex_3_6.sce @@ -0,0 +1,10 @@ +//example-3.6
+//caption-solution by gauss elimibnation method
+
+A=[1 1 1;3 3 4;2 1 3] //matrices A and b from the above
+ //system of equations
+
+b=[6;20;13]
+
+pivotgausselim(A,b) //call gauss elimination function to solve the
+ //matrices A and b
\ No newline at end of file diff --git a/50/CH3/EX3.8/ex_3_8.sce b/50/CH3/EX3.8/ex_3_8.sce new file mode 100755 index 000000000..88c13b1ea --- /dev/null +++ b/50/CH3/EX3.8/ex_3_8.sce @@ -0,0 +1,8 @@ +//example no. 3.8
+// solving the matrix equation with partial pivoting in gauss elimination
+
+A=[2 1 1 -2;4 0 2 1;3 2 2 0;1 3 2 -1]
+
+b=[-10;8;7;-5]
+
+pivotgausselim(A,b)
\ No newline at end of file diff --git a/50/CH3/EX3.9/ex_3_9.sce b/50/CH3/EX3.9/ex_3_9.sce new file mode 100755 index 000000000..34c7609a3 --- /dev/null +++ b/50/CH3/EX3.9/ex_3_9.sce @@ -0,0 +1,15 @@ +//example no.3.9
+//solving the system using inverse of the cofficient matrix
+
+A=[1 1 1;4 3 -1;3 5 3]
+
+I=[1 0 0;0 1 0;0 0 1]
+
+b=[1 ;6 ;4]
+
+M=jorden(A,I)
+
+IA=M(1:3,4:6)
+
+X=IA*b
+
diff --git a/50/CH4/EX4.15/ex_4_15.sce b/50/CH4/EX4.15/ex_4_15.sce new file mode 100755 index 000000000..be4304926 --- /dev/null +++ b/50/CH4/EX4.15/ex_4_15.sce @@ -0,0 +1,26 @@ +// example 4.15:
+// obtain the interpolate using backward differences polinomial
+
+
+xL=[.1 .2 .3 .4 .5 ]'
+
+f=[1.4 1.56 1.76 2 2.28]'
+n=2;
+
+
+// hence;
+disp('P=1.4+(x-.1)*(.16/.1)+(x-.5)*(x-.4)*(.04/.02)')
+disp('P=2x^2+x+1.28');
+
+// 1) obtain the interpolate at x=0.25;
+x=0.25;
+[P]=NBDP(x,n,xL,f);
+P
+disp('f(.25)=1.655');
+
+
+// 2) obtain the interpolate at x=0.35;
+x=0.35;
+[P]=NBDP(x,n,xL,f);
+P
+disp('f(.35)=1.875');
\ No newline at end of file diff --git a/50/CH4/EX4.20/ex_4_20.sce b/50/CH4/EX4.20/ex_4_20.sce new file mode 100755 index 000000000..77cba39f6 --- /dev/null +++ b/50/CH4/EX4.20/ex_4_20.sce @@ -0,0 +1,14 @@ +// example 4.20;
+// hermite interpolation:
+
+x=[-1 0 1];
+
+f=[1 1 3];
+
+fp=[-5 1 7];
+
+ P= hermiteinterpol(x,f,fp);
+
+// hence;
+disp('f(-0.5)=3/8');
+disp('f(0.5)=11/8');
\ No newline at end of file diff --git a/50/CH4/EX4.21/ex_4_21.sce b/50/CH4/EX4.21/ex_4_21.sce new file mode 100755 index 000000000..288cd93e5 --- /dev/null +++ b/50/CH4/EX4.21/ex_4_21.sce @@ -0,0 +1,17 @@ +// example: 4.21;
+// piecewise linear interpolating polinomials:
+
+x1=1;x2=2; x3=4;x4=8;
+f1=3; f2=7; f3=21; f4=73;
+// we need to apply legranges interpolation in sub-ranges [1,2];[2,4],[4,8];
+
+ x=poly(0,"x");
+
+P1=legrangeinterpol (x1,x2,f1,f2); // in the range [1,2]
+P1
+
+P1=legrangeinterpol (x2,x3,f2,f3); // in the range [2,4]
+P1
+
+P1=legrangeinterpol (x3,x4,f3,f4); // in the range [4,8]
+P1
\ No newline at end of file diff --git a/50/CH4/EX4.22/ex_4_22.sce b/50/CH4/EX4.22/ex_4_22.sce new file mode 100755 index 000000000..972fd5058 --- /dev/null +++ b/50/CH4/EX4.22/ex_4_22.sce @@ -0,0 +1,30 @@ +// example: 4.22;
+// piecewise quadratic interpolating polinomials:
+
+X=[-3 -2 -1 1 3 6 7];
+F=[369 222 171 165 207 990 1779];
+// we need to apply legranges interpolation in sub-ranges [-3 ,-1];[-1,3],[3,7];
+
+ x=poly(0,"x");
+
+ // 1) in the range [-3,-1]
+ x=[-3 -2 -1];
+ f=[369 222 171];
+ n=2;
+P2=lagrangefundamentalpoly(x,f,n);
+
+ // 2) in the range [-1,3]
+ x=[-1 1 3];
+ f=[171 165 207];
+n=2;
+P2=lagrangefundamentalpoly(x,f,n)
+
+ // 3) in the range [3,7]
+ x=[3 6 7];
+ f=[207 990 1779];
+n=2;
+P2=lagrangefundamentalpoly(x,f,n)
+
+
+
+// hence, we obtain the values of f(-2.5)=48; f(6.5)=1351.5;
\ No newline at end of file diff --git a/50/CH4/EX4.23/ex_4_23.sce b/50/CH4/EX4.23/ex_4_23.sce new file mode 100755 index 000000000..e06601bef --- /dev/null +++ b/50/CH4/EX4.23/ex_4_23.sce @@ -0,0 +1,26 @@ +// example: 4.23;
+// piecewise cubical interpolating polinomials:
+
+X=[-3 -2 -1 1 3 6 7];
+F=[369 222 171 165 207 990 1779];
+// we need to apply legranges interpolation in sub-ranges [-3 ,1];[1,7];
+
+
+ x=poly(0,"x");
+
+ // 1) in the range [-3,1]
+ x=[-3 -2 -1 1];
+ f=[369 222 171 165];
+ n=3;
+P2=lagrangefundamentalpoly(x,f,n);
+
+ // 2) in the range [1,7]
+ x=[1 3 6 7];
+ f=[165 207 990 1779];
+n=3;
+P2=lagrangefundamentalpoly(x,f,n)
+
+
+
+ // hence,
+ disp('f(6.5)=1339.25');
\ No newline at end of file diff --git a/50/CH4/EX4.3/ex_4_3.sce b/50/CH4/EX4.3/ex_4_3.sce new file mode 100755 index 000000000..fb99fe436 --- /dev/null +++ b/50/CH4/EX4.3/ex_4_3.sce @@ -0,0 +1,19 @@ + // example 4.3
+ // find the linear interpolation polinomial
+ // using
+
+ disp('f(2)=4');
+ disp('f(2.5)=5.5');
+ // 1)lagrange interpolation,
+
+ P1=legrangeinterpol (2,2.5,4,5.5)
+ // 2)aitken's iterated interpolation,
+
+ P1=aitkeninterpol (2,2.5,4,5.5)
+
+ // 3) newton devided differance interpolation,
+
+ P1=NDDinterpol (2,2.5,4,5.5)
+
+ // hence approximate value of f(2.2)= 4.6;
+
\ No newline at end of file diff --git a/50/CH4/EX4.31/ex_4_31.sce b/50/CH4/EX4.31/ex_4_31.sce new file mode 100755 index 000000000..117c4dfb9 --- /dev/null +++ b/50/CH4/EX4.31/ex_4_31.sce @@ -0,0 +1,23 @@ +// example 4.31
+// obtain the linear polinomial approximation to the function f(x)=x^3
+
+// let P(x)=a0*x+a1
+
+// hence I(a0,a1)= integral (x^3-(a0*x+a1))^2 in the interval [0,1]
+
+printf('I=1/7-2*(a0/5+a1/4)+a0^2/3+a0*a1+a1^2')
+printf('dI/da0 = -2/5+2/3*a0+a1=0')
+
+printf('dI/da1 = -1/2+a0+2*a1=0')
+
+// hence
+
+printf('[2/3 1;1 2]*[a0 ;a1]=[2/5; 1/2]')
+
+// solving for a0 and a1;
+
+a0=9/10;
+a1=-1/5;
+// hence considering the polinomial with intercept P1(x)=(9*x-2)/10;
+
+// considering the polinomial approximation through origin P2(x)=3*x/5;
\ No newline at end of file diff --git a/50/CH4/EX4.32/ex_4_32.sce b/50/CH4/EX4.32/ex_4_32.sce new file mode 100755 index 000000000..81c43db58 --- /dev/null +++ b/50/CH4/EX4.32/ex_4_32.sce @@ -0,0 +1,45 @@ +// example 4.32
+// obtain the linear polinomial approximation to the function f(x)=x^1/2
+
+// let P(x)=a0*x+a1
+
+
+// for n=1;
+// hence I(c0,c1)= integral (x^1/2-(c1*x+c0))^2 in the interval [0,1]
+
+
+printf('dI/dc0 = -2*(2/3-c0-c1/2)=0')
+
+printf('dI/dc1 =-2*(2/5-c0/2-c1/3) =0')
+
+// hence
+
+printf('[1 1/2;1/2 1/3]*[c0 ;c1]=[-4/3; -4/5]')
+
+// hence solving for c0 and c1;
+
+
+// the first degree square approximation P(x)=4*(1+3*x)/15;
+
+// for n=2;
+
+// hence I(c0,c1,c2)= integral (x^1/2-(c2*x^2+c1*x+c0))^2 in the interval [0,1]
+
+
+printf('dI/dc0 = (2/3-c0-c1/2-c2/2)=0')
+
+printf('dI/dc1 =(2/5-c0/2-c1/3-c2/4) =0')
+
+printf('dI/dc2 =(2/7-c0/3-c1/4-c2/5) =0')
+
+
+// hence
+
+printf('[1 1/2 1/2;1/2 1/3 1/4;1/3 1/4 1/5]*[c0 ;c1;c2]=[-2/3; -2/5;-2/7]')
+
+// hence solving for c0,c1 and c2;
+
+
+// the first degree square approximation P(x)=(6+48*x-20*x^2)/35;
+
+
diff --git a/50/CH4/EX4.34/ex_4_34.sce b/50/CH4/EX4.34/ex_4_34.sce new file mode 100755 index 000000000..ed9752978 --- /dev/null +++ b/50/CH4/EX4.34/ex_4_34.sce @@ -0,0 +1,9 @@ +// example 4.34
+
+// obtain least square straight line fit
+
+x=[.2 .4 .6 .8 1];
+f=[.447 .632 .775 .894 1];
+
+
+[P]=straightlineapprox(x,f) // call of the function to get the desired solution
\ No newline at end of file diff --git a/50/CH4/EX4.35/ex_4_35.sce b/50/CH4/EX4.35/ex_4_35.sce new file mode 100755 index 000000000..d7dceb169 --- /dev/null +++ b/50/CH4/EX4.35/ex_4_35.sce @@ -0,0 +1,7 @@ +// example 4.35
+
+// obtain least square approximation of second degree;
+x=[-2 -1 0 1 2];
+f=[15 1 1 3 19];
+
+[P]=quadraticapprox(x,f) // call of the function to get the desired solution
\ No newline at end of file diff --git a/50/CH4/EX4.36/ex_4_36.sce b/50/CH4/EX4.36/ex_4_36.sce new file mode 100755 index 000000000..afc372683 --- /dev/null +++ b/50/CH4/EX4.36/ex_4_36.sce @@ -0,0 +1,31 @@ +// example 4.36
+// method of least squares to fit the data to the curve P(x)=c0*X+c1/squrt(X)
+
+x=[.2 .3 .5 1 2];
+f=[16 14 11 6 3];
+
+// I(c0,c1)= summation of (f(x)-(c0*X+c1/sqrt(X)))^2
+
+// hence on parcially derivating the summation,
+
+n=length(x);m=length(f);
+if m<>n then
+ error('linreg - Vectors x and f are not of the same length.');
+ abort;
+end;
+
+s1=0; // s1= summation of x(i)*f(i)
+s2=0; // s2= summation of f(i)/sqrt(x(i))
+s3=0;
+for i=1:n
+ s1=s1+x(i)*f(i);
+ s2=s2+f(i)/sqrt(x(i));
+ s3=s3+1/x(i);
+end
+
+c0=det([s1 sum(sqrt(x));s2 s3])/det([sum(x^2) sum(sqrt(x));sum(sqrt(x)) s3])
+
+c1=det([sum(x^2) s1;sum(sqrt(x)) s2])/det([sum(x^2) sum(sqrt(x));sum(sqrt(x)) s3])
+X=poly(0,"X");
+P=c0*X+c1/X^1/2
+// hence considering the polinomial P(x)=7.5961*X^1/2-1.1836*X
\ No newline at end of file diff --git a/50/CH4/EX4.37/ex_4_37.sce b/50/CH4/EX4.37/ex_4_37.sce new file mode 100755 index 000000000..067c73c44 --- /dev/null +++ b/50/CH4/EX4.37/ex_4_37.sce @@ -0,0 +1,30 @@ +// example 4.37
+// method of least squares to fit the data to the curve P(x)=a*%e^(-3*t)+b*%e^(-2*t);
+
+t=[.1 .2 .3 .4];
+f=[.76 .58 .44 .35];
+
+// I(c0,c1)= summation of (f(x)-a*%e^(-3*t)+b*%e^(-2*t))
+
+// hence on parcially derivating the summation,
+
+n=length(t);m=length(f);
+if m<>n then
+ error('linreg - Vectors t and f are not of the same length.');
+ abort;
+end;
+
+s1=0; // s1= summation of f(i)*%e^(-3*t(i));
+s2=0; // s2= summation of f(i)*%e^(-2*t(i));
+
+for i=1:n
+ s1=s1+f(i)*%e^(-3*t(i));
+ s2=s2+f(i)*%e^(-2*t(i));
+
+end
+
+a=det([s1 sum(%e^(-5*t)); s2 sum(%e^(-4*t))])/det([sum(%e^(-6*t)) sum(%e^(-5*t)); sum(%e^(-5*t)) sum(%e^(-4*t))])
+
+b=det([sum(%e^(-6*t)) s1; sum(%e^(-5*t)) s2])/det([sum(%e^(-6*t)) sum(%e^(-5*t)); sum(%e^(-5*t)) sum(%e^(-4*t))])
+
+// hence considering the polinomial P(t)=.06853*%e^(-3*t)+0.3058*%e^(-2*t)
\ No newline at end of file diff --git a/50/CH4/EX4.38/ex_4_38.sce b/50/CH4/EX4.38/ex_4_38.sce new file mode 100755 index 000000000..96d5857a9 --- /dev/null +++ b/50/CH4/EX4.38/ex_4_38.sce @@ -0,0 +1,26 @@ +// example 4.38
+// gram schmidt orthogonalisation
+
+W=1;
+ x=poly(0,"x");
+ P0=1;
+ phi0=P0;
+ a10=integrate('W*x*phi0','x',0,1)/integrate('W*1*phi0','x',0,1)
+ P1=x-a10*phi0
+ phi1=P1;
+
+ a20=integrate('W*x^2*phi0','x',0,1)/integrate('W*1*phi0','x',0,1)
+
+ a21=integrate('(x^2)*(x-1/2)','x',0,1)/integrate('(x-1/2)^2','x',0,1)
+
+ P2=x^2-a20*x-a21*phi1
+
+// since ,I= intgral [x^(1/2)-c0*P0-c1*P1-c2*P2]^2 inthe range [0,1]
+
+// hence partially derivating I
+
+c0=integrate('x^(1/2)','x',0,1)/integrate('1','x',0,1)
+c1=integrate('(x^(1/2))*(x-(1/2))','x',0,1)/integrate('(x-(1/2))^2','x',0,1)
+c1=integrate('(x^(1/2))*(x^2-4*x/3+1/2)','x',0,1)/integrate('(x^2-4*x/3+1/2)^2','x',0,1)
+
+
diff --git a/50/CH4/EX4.39/ex_4_39.sce b/50/CH4/EX4.39/ex_4_39.sce new file mode 100755 index 000000000..3b231135d --- /dev/null +++ b/50/CH4/EX4.39/ex_4_39.sce @@ -0,0 +1,34 @@ +// example 4.39
+// gram schmidt orthogonalisation
+
+// 1)
+ W=1;
+ x=poly(0,"x");
+ P0=1
+ phi0=P0;
+ a10=0;
+ P1=x-a10*phi0
+ phi1=P1;
+
+ a20=integrate('x^2','x',-1,1)/integrate('W*1*phi0','x',-1,1);
+
+ a21=integrate('(x^3)','x',-1,1)/integrate('(x)^2','x',-1,1);
+
+ P2=x^2-a20*x-a21*phi1
+
+
+// 2)
+disp(' W=1/(1-x^2)^(1/2)');
+ x=poly(0,"x");
+ P0=1
+ phi0=P0;
+ a10=0;
+ P1=x-a10*phi0
+ phi1=P1;
+
+ a20=integrate('x^2/(1-x^2)^(1/2)','x',-1,1)/integrate('1/(1-x^2)^(1/2)','x',-1,1);
+
+ a21=0; // since x^3 is an odd function;
+
+ P2=x^2-a20*x-a21*phi1
+
diff --git a/50/CH4/EX4.4/ex_4_4.sce b/50/CH4/EX4.4/ex_4_4.sce new file mode 100755 index 000000000..c31365e30 --- /dev/null +++ b/50/CH4/EX4.4/ex_4_4.sce @@ -0,0 +1,12 @@ + // example 4.4
+ // find the linear interpolation polinomial
+ // using lagrange interpolation,
+
+ disp('sin(.1)=.09983; sin(.2)=.19867');
+
+
+ P1=legrangeinterpol (.1,.2,.09983,.19867)
+
+// hence;
+disp('P(.15)=.00099+.9884*.15')
+disp('P(0.15)=0.14925');
\ No newline at end of file diff --git a/50/CH4/EX4.41/ex_4_41.sce b/50/CH4/EX4.41/ex_4_41.sce new file mode 100755 index 000000000..1ce2864ac --- /dev/null +++ b/50/CH4/EX4.41/ex_4_41.sce @@ -0,0 +1,20 @@ +// example 4.41
+// using chebyshev polinomials obtain least squares approximation of second degree;
+
+// the chebeshev polinomials;
+x=poly(0,"x");
+T0=1;
+T1=x;
+T2=2*x^2-1;
+
+
+// I=integrate ('1/(1-x^2)^(1/2)*(x^4-c0*T0-c1*T1-c2*T2)^2','x',-1,1)
+
+// since;
+c0=integrate('(1/3.14)*(x^4)/(1-x^2)^(1/2)','x',-1,1)
+
+c1=integrate('(2/3.14)*(x^5)/(1-x^2)^(1/2)','x',-1,1)
+
+c2=integrate('(2/3.14)*(x^4)*(2*x^2-1)/(1-x^2)^(1/2)','x',-1,1)
+
+ f=(3/8)*T0+(1/2)*T2;
\ No newline at end of file diff --git a/50/CH4/EX4.6/ex_4_6.sce b/50/CH4/EX4.6/ex_4_6.sce new file mode 100755 index 000000000..869b1f772 --- /dev/null +++ b/50/CH4/EX4.6/ex_4_6.sce @@ -0,0 +1,18 @@ +// example :4.6
+// caption : obtain the legrange linear interpolating polinomial
+
+
+// 1) obtain the legrange linear interpolating polinomial in the interval [1,3] and obtain approximate value of f(1.5),f(2.5);
+x0=1;x1=2;x2=3;f0=.8415;f1=.9093;f2=.1411;
+
+P13=legrangeinterpol (x0,x2,f0,f2) //in the range [1 ,3]
+
+
+P12=legrangeinterpol (x0,x1,f0,f1) // in the range [1,2]
+
+P23=legrangeinterpol (x1,x2,f1,f2) // inthe range [2,3]
+
+// from P23 we find that ; where as exact value is sin(2.5)=0.5985;
+disp('P(1.5)=0.8754');
+disp('exact value of sin(1.5)=.9975')
+disp('P(2.5)=0.5252');
\ No newline at end of file diff --git a/50/CH4/EX4.7/ex_4_7.sce b/50/CH4/EX4.7/ex_4_7.sce new file mode 100755 index 000000000..cbc217059 --- /dev/null +++ b/50/CH4/EX4.7/ex_4_7.sce @@ -0,0 +1,15 @@ + // example 4.7
+ // polinomial of degree 2;
+
+ // f(0)=1;f(1)=3;f(3)=55;
+
+ // using legrange fundamental polinomial rule,
+
+ x=[0 1 3]; // arrainging the inputs of the function as elements of a row,
+ f=[1 3 55]; // arrainging the outputs of the function as elements of a row,
+ n=2; // degree of the polinomial;
+
+
+ P2=lagrangefundamentalpoly(x,f,n)
+
+
\ No newline at end of file diff --git a/50/CH4/EX4.8/ex_4_8.sce b/50/CH4/EX4.8/ex_4_8.sce new file mode 100755 index 000000000..6bfdd8e78 --- /dev/null +++ b/50/CH4/EX4.8/ex_4_8.sce @@ -0,0 +1,15 @@ +
+ // example 4.8
+ // caption: solution by quadratic interpolation;
+
+ // x-degrees:[10 20 30]
+ // hence x in radians is
+ x=[3.14/18 3.14/9 3.14/6];
+ f=[1.1585 1.2817 1.3660];
+ n=2;
+
+
+ P2=lagrangefundamentalpoly(x,f,n)
+
+ // hence from P2 ,the exact value of f(3.14/12) is 1.2246;
+ // where as exact value is 1.2247;
\ No newline at end of file diff --git a/50/CH4/EX4.9/ex_4_9.sce b/50/CH4/EX4.9/ex_4_9.sce new file mode 100755 index 000000000..726da7ed5 --- /dev/null +++ b/50/CH4/EX4.9/ex_4_9.sce @@ -0,0 +1,16 @@ +// example 4.9
+// caption :obtain the polinomial of degree 2
+
+x=[0 1 3];
+f=[1 3 55];
+n=2;
+
+// 1) iterated interpolation;
+
+
+[L012,L02,L01]=iteratedinterpol (x,f,n)
+
+// 2) newton divided diffrences interpolation;
+
+
+ P2=NDDinterpol2 (x,f)
\ No newline at end of file diff --git a/50/CH5/EX5.1/ex_5_1.sce b/50/CH5/EX5.1/ex_5_1.sce new file mode 100755 index 000000000..4ffe53138 --- /dev/null +++ b/50/CH5/EX5.1/ex_5_1.sce @@ -0,0 +1,12 @@ +// example: 5.1
+// linear and quadratic interpolation:
+
+// f(x)=ln x;
+
+xL=[2 2.2 2.6];
+f=[.69315 .78846 .95551];
+
+// 1) fp(2) with linear interpolation;
+
+fp=linearinterpol(xL,f);
+disp(fp);
diff --git a/50/CH5/EX5.10/ex_5_10.sce b/50/CH5/EX5.10/ex_5_10.sce new file mode 100755 index 000000000..919949c50 --- /dev/null +++ b/50/CH5/EX5.10/ex_5_10.sce @@ -0,0 +1,15 @@ +// example 5.10;
+// find the jacobian matrix;
+
+
+// given two functions in x,y;
+// and the point at which the jacobian has to be found out;
+
+deff('[w]=f1(x,y)','w=x^2+y^2-x');
+
+deff('[q]=f2(x,y)','q=x^2-y^2-y');
+
+h=1;k=1;
+
+J= jacobianmat (f1,f2,h,k);
+disp(J);
\ No newline at end of file diff --git a/50/CH5/EX5.11/ex_5_11.sce b/50/CH5/EX5.11/ex_5_11.sce new file mode 100755 index 000000000..3f861a3ea --- /dev/null +++ b/50/CH5/EX5.11/ex_5_11.sce @@ -0,0 +1,18 @@ +// example : 5.11
+// solve the definite integral by 1) trapezoidal rule, 2)simpsons rule
+// exact value of the integral is ln 2= 0.693147,
+
+deff('[y]=F(x)','y=1/(1+x)')
+
+//1)trapezoidal rule,
+
+a=0;
+b=1;
+I =trapezoidal(0,1,F)
+disp(error =.75-.693147)
+
+// simpson's rule
+
+I=simpson(a,b,F)
+
+disp(error =.694444-.693147)
\ No newline at end of file diff --git a/50/CH5/EX5.12/ex_5_12.sce b/50/CH5/EX5.12/ex_5_12.sce new file mode 100755 index 000000000..2e08e16f7 --- /dev/null +++ b/50/CH5/EX5.12/ex_5_12.sce @@ -0,0 +1,17 @@ +// example 5.12
+// caption: solve the integral by 1)mid-point rule,2)two-point open type rule
+
+
+// let integration of f(x)=sin(x)/(x) in the range [0,1] is equal to I1 and I2
+// 1)mid -point rule;
+a=0;b=1;
+h=(b-a)/2;
+
+x=0:h:1;
+deff('[y]=f(x)','y=sin(x)/x')
+I1=2*h*f(x(1)+h)
+
+
+//2) two-point open type rule
+h=(b-a)/3;
+I2=(3/2)*h*(f(x(1)+h)+f(x(1)+2*h))
\ No newline at end of file diff --git a/50/CH5/EX5.13/ex_5_13.sce b/50/CH5/EX5.13/ex_5_13.sce new file mode 100755 index 000000000..d0ada71c8 --- /dev/null +++ b/50/CH5/EX5.13/ex_5_13.sce @@ -0,0 +1,12 @@ +// example 5.13
+// caption: simpson 3-8 rule
+
+
+// let integration of f(x)=1/(1+x) in the range [0,1] by simpson 3-8 rule is equal to I
+
+x=0:1/3:1;
+deff('[y]=f(x)','y=1/(1+x)')
+
+[I] = simpson38(x,f)
+
+
\ No newline at end of file diff --git a/50/CH5/EX5.15/ex_5_15.sce b/50/CH5/EX5.15/ex_5_15.sce new file mode 100755 index 000000000..b56793e06 --- /dev/null +++ b/50/CH5/EX5.15/ex_5_15.sce @@ -0,0 +1,26 @@ +// example :5.15
+// find the quadrature formula of
+// integral of f(x)*(1/sqrt(x(1+x))) in the range [0,1]= a1*f(0)+a2*f(1/2)+a3*f(1)=I
+// hence find integral 1/sqrt(x-x^3) in the range [0,1]
+
+// making the method exact for polinomials of degree upto 2,
+// I=I1=a1+a2+a3
+// I=I2=(1/2)*a2+a3
+// I=I3=(1/4)*a2+a3
+
+// A=[a1 a2 a3]'
+
+I1=integrate('1/sqrt(x*(1-x))','x',0,1)
+I2=integrate('x/sqrt(x*(1-x))','x',0,1)
+I3=integrate('x^2/sqrt(x*(1-x))','x',0,1)
+
+//hence
+// [1 1 1;0 1/2 1 ;0 1/4 1]*A=[I1 I2 I3]'
+
+A=inv([1 1 1;0 1/2 1 ;0 1/4 1])*[I1 I2 I3]'
+// I=(3.14/4)*(f(0)+2*f(1/2)+f(1));
+
+// hence, for solving the integral 1/sqrt(x-x^3) in the range [0,1]=I
+
+deff('[y]=f(x)','y=1/sqrt(1+x)');
+I=(3.14/4)*[1+2*sqrt(2/3)+sqrt(2)/2]
\ No newline at end of file diff --git a/50/CH5/EX5.16/ex_5_16.sce b/50/CH5/EX5.16/ex_5_16.sce new file mode 100755 index 000000000..f3b260f89 --- /dev/null +++ b/50/CH5/EX5.16/ex_5_16.sce @@ -0,0 +1,16 @@ +// example 5.16
+// caption: gauss-legendre three point method
+// I= integral 1/(1+x) in the range [0,1];
+// first we need ti transform the interval [0,1 ] to [-1,1], since gauss-legendre three point method is applicable in the range[-1,1],
+
+// let t=ax+b;
+// solving for a,b from the two ranges, we get a=2; b=-1; t=2x-1;
+
+// hence I=integral 1/(1+x) in the range [0,1]= integral 1/(t+3) in the range [-1,1];
+
+
+deff('[y]=f(t)','y=1/(t+3)');
+// since , from gauss legendre three point rule(n=2);
+I=(1/9)*(5*f(-sqrt(3/5))+8*f(0)+5*f(sqrt(3/5)))
+
+// we know , exact solution is ln 2=0.693147;
\ No newline at end of file diff --git a/50/CH5/EX5.17/ex_5_17.sce b/50/CH5/EX5.17/ex_5_17.sce new file mode 100755 index 000000000..f0e6d0b6f --- /dev/null +++ b/50/CH5/EX5.17/ex_5_17.sce @@ -0,0 +1,24 @@ +// example 5.17
+// caption: gauss-legendre method
+// I= integral 2*x/(1+x^4) in the range [1,2];
+// first we need ti transform the interval [1,2 ] to [-1,1], since gauss-legendre three point method is applicable in the range[-1,1],
+
+// let t=ax+b;
+// solving for a,b from the two ranges, we get a=1/2; b=3/2; x=(t+3)/2;
+
+// hence I=integral 2*x/(1+x^4) in the range [0,1]= integral 8*(t+3)/16+(t+3)^4 in the range [-1,1];
+
+
+deff('[y]=f(t)','y=8*(t+3)/(16+(t+3)^4) ');
+
+// 1) since , from gauss legendre one point rule;
+I1=2*f(0)
+
+// 2) since , from gauss legendre two point rule;
+I2=f(-1/sqrt(3))+f(1/sqrt(3))
+
+// 3) since , from gauss legendre three point rule;
+I=(1/9)*(5*f(-sqrt(3/5))+8*f(0)+5*f(sqrt(3/5)))
+
+
+// we know , exact solution is 0.5404;
\ No newline at end of file diff --git a/50/CH5/EX5.18/ex_5_18.sce b/50/CH5/EX5.18/ex_5_18.sce new file mode 100755 index 000000000..08ce9df7a --- /dev/null +++ b/50/CH5/EX5.18/ex_5_18.sce @@ -0,0 +1,23 @@ +// example 5.18
+// caption: gauss-chebyshev method
+
+// we write the integral as I=integral f(x)/sqrt(1-x^2) in the range [-1,1];
+// where f(x)=(1-x^2)^2*cos(x)
+
+deff('[y]=f(x)','y=(1-x^2)^2*cos(x)');
+
+// 1) since , from gauss chebyshev one point rule;
+I1=(3.14)*f(0)
+
+// 2) since , from gauss chebyshev two point rule;
+I2=(3.14/2)*f(-1/sqrt(2))+f(1/sqrt(2))
+
+// 3) since , from gauss chebyshev three point rule;
+I=(3.14/3)*(f(-sqrt(3)/2)+f(0)+f(sqrt(3)/2))
+
+
+// and 4) since , from gauss legendre three point rule;
+I=(1/9)*(5*f(-sqrt(3/5))+8*f(0)+5*f(sqrt(3/5)))
+
+
+
diff --git a/50/CH5/EX5.2/ex_5_2.sce b/50/CH5/EX5.2/ex_5_2.sce new file mode 100755 index 000000000..cd5e0032e --- /dev/null +++ b/50/CH5/EX5.2/ex_5_2.sce @@ -0,0 +1,9 @@ +// example 5.2
+// evaluate fp(.8) and fpp(.8) with quadratic interpolation;
+
+xL=[.4 .6 .8];
+f=[.0256 .1296 .4096];
+h=.2;
+
+fp=(1/2*h)*(f(1)-4*f(2)+3*f(3))
+fpp=(1/h^2)*(f(1)-2*f(2)+f(3)) // from equation 5.22c and 5.24c in the book;
diff --git a/50/CH5/EX5.20/ex_5_20.sce b/50/CH5/EX5.20/ex_5_20.sce new file mode 100755 index 000000000..3e0f84277 --- /dev/null +++ b/50/CH5/EX5.20/ex_5_20.sce @@ -0,0 +1,16 @@ +// example 5.20
+// caption: gauss-leguerre method
+// I= integral e^-x/(1+x^2) in the range [0,~];
+
+// observing the integral we can inffer that f(x)=1/(1+x^2)
+
+deff('[y]=f(x)','y=1/(1+x^2) ');
+
+
+// 1) since , from gauss leguerre two point rule;
+I2=(1/4)*[(2+sqrt(2))*f(2-sqrt(2))+(2-sqrt(2))*f(2+sqrt(2))]
+
+// 3) since , from gauss leguerre three point rule;
+I=(0.71109*f(0.41577)+0.27852*f(2.29428)+0.01039*f(6.28995))
+
+
diff --git a/50/CH5/EX5.21/ex_5_21.sce b/50/CH5/EX5.21/ex_5_21.sce new file mode 100755 index 000000000..64303a5b5 --- /dev/null +++ b/50/CH5/EX5.21/ex_5_21.sce @@ -0,0 +1,16 @@ +// example 5.21
+// caption: gauss-leguerre method
+// I= integral e^-x*(3*x^3-5*x+1) in the range [0,~];
+
+// observing the integral we can inffer that f(x)=(3*x^3-5*x+1)
+
+deff('[y]=f(x)','y=(3*x^3-5*x+1) ');
+
+
+// 1) since , from gauss leguerre two point rule;
+I2=(1/4)*[(2+sqrt(2))*f(2-sqrt(2))+(2-sqrt(2))*f(2+sqrt(2))]
+
+// 3) since , from gauss leguerre three point rule;
+I3=(0.71109*f(0.41577)+0.27852*f(2.29428)+0.01039*f(6.28995))
+
+
diff --git a/50/CH5/EX5.22/ex_5_22.sce b/50/CH5/EX5.22/ex_5_22.sce new file mode 100755 index 000000000..7afdb7036 --- /dev/null +++ b/50/CH5/EX5.22/ex_5_22.sce @@ -0,0 +1,20 @@ +// example 5.22
+// caption: gauss-leguerre method
+// I= integral 1/(x^2+2*x+2) in the range [0,~];
+
+// since in the gauss-leguerre method the integral would be of the form e^x*f(x);
+
+// observing the integral we can inffer that f(x)=%e^x/(x^2+2*x+2)
+deff('[y]=f(x)','y=%e^x/(x^2+2*x+2) ');
+
+
+// 1) since , from gauss leguerre two point rule;
+I2=(1/4)*[(2+sqrt(2))*f(2-sqrt(2))+(2-sqrt(2))*f(2+sqrt(2))]
+
+// 3) since , from gauss leguerre three point rule;
+I=(0.71109*f(0.41577)+0.27852*f(2.29428)+0.01039*f(6.28995))
+
+
+// the exact solution is given by,
+
+I=integrate('1/((x+1)^2+1)','x',0,1000) // 1000 ~infinite;
diff --git a/50/CH5/EX5.26/ex_5_26.sce b/50/CH5/EX5.26/ex_5_26.sce new file mode 100755 index 000000000..402f55cb0 --- /dev/null +++ b/50/CH5/EX5.26/ex_5_26.sce @@ -0,0 +1,42 @@ +// Example 5.26
+// caption: 1) composite trapizoidal rule, 2) composite simpsons rule with 2,4 ,8 equal sub-intervals,
+
+// I=integral 1/(1+x) in the range [0,1]
+
+deff('[y]=f(x)','y=1/(1+x)')
+
+// when N=2;
+// 1)composite trapizoidal rule
+h=1/2;
+x=0:h:1;
+
+ IT=comptrapezoidal(x,h,f)
+
+ // 2)composite simpsons rule
+
+ [I] = simpson13(x,h,f)
+
+
+ // when N=4
+ // 1)composite trapizoidal rule
+h=1/4;
+x=0:h:1;
+
+ IT=comptrapezoidal(x,h,f)
+
+ // 2)composite simpsons rule
+
+ [I] = simpson13(x,h,f)
+
+
+
+ // when N=8
+ // 1)composite trapizoidal rule
+h=1/8;
+x=0:h:1;
+
+ IT=comptrapezoidal(x,h,f)
+
+ // 2)composite simpsons rule
+
+ [I] = simpson13(x,h,f)
\ No newline at end of file diff --git a/50/CH5/EX5.27/ex_5_27.sce b/50/CH5/EX5.27/ex_5_27.sce new file mode 100755 index 000000000..cd5800369 --- /dev/null +++ b/50/CH5/EX5.27/ex_5_27.sce @@ -0,0 +1,23 @@ +// example 5.27
+// caption: gauss-legendre three point method
+// I= integral 1/(1+x) in the range [0,1];
+
+// we are asked to subdivide the range into two,
+// first we need to sub-divide the interval [0,1 ] to [0,1/2] and [1/2,1] and then transform both to [-1,1], since gauss-legendre three point method is applicable in the range[-1,1],
+
+// t=4x-1 and y=4x-3;
+
+// hence I=integral 1/(1+x) in the range [0,1]= integral 1/(t+5) in the range [-1,1]+ integral 1/(t+7) in the range [-1,1]
+
+
+deff('[y1]=f1(t)','y1=1/(t+5)');
+// since , from gauss legendre three point rule(n=2);
+I1=(1/9)*(5*f1(-sqrt(3/5))+8*f1(0)+5*f1(sqrt(3/5)))
+
+deff('[y2]=f2(t)','y2=1/(t+7)');
+// since , from gauss legendre three point rule(n=2);
+I2=(1/9)*(5*f2(-sqrt(3/5))+8*f2(0)+5*f2(sqrt(3/5)))
+
+I=I1+I2
+
+// we know , exact solution is .693147;
\ No newline at end of file diff --git a/50/CH5/EX5.29/ex_5_29.sce b/50/CH5/EX5.29/ex_5_29.sce new file mode 100755 index 000000000..7e88c2c8c --- /dev/null +++ b/50/CH5/EX5.29/ex_5_29.sce @@ -0,0 +1,12 @@ +// example 5.29
+// evaluate the given double integral using the simpsons rule;
+
+// I= double integral f(x)=1/(x+y) in the range x=[1,2],y=[1,1.5];
+
+h=.5;
+k=.25;
+deff('[w]=f(x,y)','w=1/(x+y)')
+
+I=(.125/9)*[{f(1,1)+f(2,1)+f(1,1.5)+f(2,1.5)}+4*{f(1.5,1)+f(1,1.25)+f(1.5,1.5)+f(2,1.25)}+16*f(1.5,1.25)];
+disp(I);
+
diff --git a/50/CH5/EX5.30/ex_5_30.sce b/50/CH5/EX5.30/ex_5_30.sce new file mode 100755 index 000000000..17f88a29d --- /dev/null +++ b/50/CH5/EX5.30/ex_5_30.sce @@ -0,0 +1,18 @@ +// example 5.30
+// evaluate the given double integral using the simpsons rule;
+
+// I= double integral f(x)=1/(x+y) in the range x=[1,2],y=[1,2];
+// 1)
+h=.5;
+k=.5;
+deff('[w]=f(x,y)','w=1/(x+y)')
+
+I=(1/16)*[{f(1,1)+f(2,1)+f(1,2)+f(2,2)}+2*{f(1.5,1)+f(1,1.5)+f(2,1.5)+f(1.5,2)}+4*f(1.5,1.5)]
+
+// 2)
+h=.25;
+k=.25;
+deff('[w]=f(x,y)','w=1/(x+y)')
+
+I=(1/64)*[{f(1,1)+f(2,1)+f(1,2)+f(2,2)}+2*{f(5/4,1)+f(3/2,1)+f(7/4,1)+f(1,5/4)+f(1,3/2)+f(1,7/4)+f(2,5/4)+f(2,3/4)+f(2,7/4)+f(5/4,2)+f(3/2,2)+f(7/4,2)}+4*{f(5/4,5/4)+f(5/4,3/2)+f(5/4,7/4)+f(3/2,5/4)+f(3/2,3/2)+f(3/2,7/4)+f(7/4,5/4)+f(7/4,3/2)+f(7/4,7/4)}]
+
diff --git a/50/CH6/EX6.12/ex_6_12.sce b/50/CH6/EX6.12/ex_6_12.sce new file mode 100755 index 000000000..045ec5415 --- /dev/null +++ b/50/CH6/EX6.12/ex_6_12.sce @@ -0,0 +1,9 @@ +// example 6.12,
+// caption: solve the IVP by backward euler method,
+// with h=0.2,
+// u'=f(t,u)
+deff('[z]=f(t,u)','z=-2*t*u^2');
+
+
+[u] = backeuler(1,0,0.4,.2,f) // h=0.2;
+
diff --git a/50/CH6/EX6.13/ex_6_13.sce b/50/CH6/EX6.13/ex_6_13.sce new file mode 100755 index 000000000..7fc402959 --- /dev/null +++ b/50/CH6/EX6.13/ex_6_13.sce @@ -0,0 +1,11 @@ +// example 6.13,
+// caption: solve the IVP by euler midpoint method,
+// with h=0.2,
+// u'=f(t,u)
+deff('[z]=f(t,u)','z=-2*t*u^2');
+deff('[w]=fp(t,u)','w=-2*u^2-4*u*t');
+
+
+
+[u] = eulermidpoint(1,0,1,.2,f,fp) // h=0.2;
+
diff --git a/50/CH6/EX6.15/ex_6_15.sce b/50/CH6/EX6.15/ex_6_15.sce new file mode 100755 index 000000000..30daf63a0 --- /dev/null +++ b/50/CH6/EX6.15/ex_6_15.sce @@ -0,0 +1,20 @@ +// example 6.15
+//caption: solving ODE by tailor series method
+// u'=t^2+ u^2, u0=0;
+
+t=0;U=0; //at t=0, the value of u is 0
+U1=0; // u1 is the 1st derivatove of the funtion u
+U2=2*t+2*U*U1 // U2 ---- 2nd derivative
+U3=2+2*(U*U2+U1^2)
+U4=2*(U*U3+3*U1*U2)
+U5=2*(U*U4+4*U1*U3+3*U2^2)
+U6=2*(U*U5+5*U1*U4+10*U2*U3)
+U7=2*(U*U6+6*U1*U5+15*U2*U4+10*U3^2)
+U8=0;
+U9=0;
+U10=0;
+U11=2*(U*U10+10*U1*U9+45*U2*U8+120*U3*U7+210*U4*U6+126*U5^2)
+ // U11 is the 11th derivative of u
+
+
+taylor(1)
\ No newline at end of file diff --git a/50/CH6/EX6.17/ex_6_17.sce b/50/CH6/EX6.17/ex_6_17.sce new file mode 100755 index 000000000..9de434d2a --- /dev/null +++ b/50/CH6/EX6.17/ex_6_17.sce @@ -0,0 +1,17 @@ +// example 6.17,
+// caption: solve by 1)modified euler cauchy, 2)heun method
+
+// h=0.2
+// 1) modified euler cauchy method,
+
+// u'=f(t,u)
+// u'=-2tu^2
+deff('[z]=f(t,u)','z=-2*t*u^2');
+
+modifiedeuler(1,0,.4,.2,f) // calling the function,
+
+// 2) heun method,
+deff('[z]=f(t,u)','z=-2*t*u^2');
+
+
+heun(1,0,.4,.2,f) // calling the function,
diff --git a/50/CH6/EX6.18/ex_6_18.sce b/50/CH6/EX6.18/ex_6_18.sce new file mode 100755 index 000000000..5a8165e09 --- /dev/null +++ b/50/CH6/EX6.18/ex_6_18.sce @@ -0,0 +1,8 @@ +// example 6.18,
+// caption: use of 4th order runge kutta method,
+
+// u'=f(t,u)
+// u'=-2tu^2
+deff('[z]=f(t,u)','z=-2*t*u^2');
+
+RK4(1,0,.4,.2,f) // calling the function,
\ No newline at end of file diff --git a/50/CH6/EX6.20/ex_6_20.sce b/50/CH6/EX6.20/ex_6_20.sce new file mode 100755 index 000000000..5430115f1 --- /dev/null +++ b/50/CH6/EX6.20/ex_6_20.sce @@ -0,0 +1,18 @@ +// example no. 6.20,
+// caption: solve the system of equations
+
+//1) eulercauchy method solving simultanious ODE
+
+deff('[z]=f1(t,u,v)','z=-3*u+2*v');
+deff('[w]=f2(t,u,v)','w=3*u-4*v');
+
+
+[u,v,t] = simeulercauchy(0,.5,0,.4,.2,f1,f2)
+
+// 2) RK4 method solving simultanious ODE
+
+
+
+[u,v,t]=simRK4(0,.5,0,.4,.2,f1,f2)
+
+
diff --git a/50/CH6/EX6.21/ex_6_21.sce b/50/CH6/EX6.21/ex_6_21.sce new file mode 100755 index 000000000..32b9e0d95 --- /dev/null +++ b/50/CH6/EX6.21/ex_6_21.sce @@ -0,0 +1,27 @@ +// example no. 6.21,
+// caption: solving the IVP by implicit RK2 method
+
+// u'=f(t,u)
+// u'=-2tu^2
+//u(0)=1,h=0.2;
+t0=0;h=0.2;tn=.4;u0=1;
+deff('[z]=f(t,u)','z=-2*t*u^2');
+umaxAllowed = 1e+100;
+
+t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0;
+
+for j = 1:n-1
+ // k1=h*f(t(j)+h/2,u(j)+k1/2);
+ // conidering the IVP we can infer that the above expression in non linear in k1,
+// hence we use newton rapson method to solve for k1;
+deff('[w]=F(u2)','w=k1+h*(2*t(j)+h)*(u(j)+k1/2)^2') // u2=u(2)
+deff('[x]=Fp(u2)','x=1+h*(2*t(j)+h)*(u(j)+k1/2)')
+
+k1=h*f(t(j),u(j));
+
+newton(k1,F,Fp);
+ u(j+1) = u(j) +k1
+ disp(u(j+1))
+
+end;
+
diff --git a/50/CH6/EX6.25/ex_6_25.sce b/50/CH6/EX6.25/ex_6_25.sce new file mode 100755 index 000000000..d9a405146 --- /dev/null +++ b/50/CH6/EX6.25/ex_6_25.sce @@ -0,0 +1,8 @@ +// example 6.25
+// caption: solving the IVP by adams-bashforth 3rd order method.
+// u'=f(t,u)
+// u'=-2tu^2
+//u(0)=1,h=0.2;
+deff('[z]=f(t,u)','z=-2*t*u^2');
+
+adamsbashforth3(1,0,1,.2,f) // calling the function,
diff --git a/50/CH6/EX6.27/ex_6_27.sce b/50/CH6/EX6.27/ex_6_27.sce new file mode 100755 index 000000000..41081a0a9 --- /dev/null +++ b/50/CH6/EX6.27/ex_6_27.sce @@ -0,0 +1,32 @@ +// example 6.27
+// solving IVP by 3rd order adams moulton
+// u'=t^2+u^2, u(1)=2,
+// h=0.1, [1,1.2]
+deff('[z]=f(t,u)','z=t^2+u^2');
+t0=1;u0=2;h=0.1;tn=1.2;
+// third order adams moulton method,
+// u(j+2)=u(j+1)+(h/12)*(5*f(t(j+2),u(j+2))+8*f(t(j+1),u(j+1))-f(t(j),u(j)));- is the expression for adamsbas-moulton3
+
+
+// on observing the IVP we can inffer that this would be a non linear equation,
+// u(j+2)=u(j+1)+(h/12)*(5*((t(j+2))^2+(u(j+2))^2)+8*((t(j+1))^2+(u(j+1))^2)-((t(j))^2+(u(j))^2))
+
+t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0;
+for j = 1:n-2
+ if j==1 then
+ k1=h*f(t(j),u(j));
+ k2=h*f(t(j)+h,u(j)+k1);
+ u(j+1) = u(j) + (k2+k1)/2;
+ disp(u(j+1))
+ end;
+end;
+
+// hence the third order adams moulton expression turns to be,
+// u(2)= 0.041667*(u(2))^2+3.194629
+// let us use newton raphsom method to solve this,
+deff('[w]=F(u2)','w=-u2+ 0.041667*(u2)^2+3.194629') // u2=u(2)
+deff('[x]=Fp(u2)','x=-1+ 0.041667*2*u2')
+
+// let us assume the initial guess of u(2)=u(1);
+
+newton(2.633333,F,Fp)
\ No newline at end of file diff --git a/50/CH6/EX6.3/ex_6_3.sce b/50/CH6/EX6.3/ex_6_3.sce new file mode 100755 index 000000000..ea77769e4 --- /dev/null +++ b/50/CH6/EX6.3/ex_6_3.sce @@ -0,0 +1,17 @@ +// example 6.3
+// solution to the given IVP
+
+disp('du/dt= A*u');
+// u=[u1 u2]';
+A=[-3 4 ;-2 3]; // given
+B=[1 0;0 1]; // identity matrix;
+
+
+
+
+[x,lam] = geigenvectors(A,B);
+
+// hence;
+disp('u=c1*%e^t*x(:,1)+c2*%e^-t*x(:,2)');
+disp('u1=c1*%e^t+c2*%e^-t*2')
+disp('u2=c1*%e^t+c2*%e^-t')
\ No newline at end of file diff --git a/50/CH6/EX6.32/ex_6_32.sce b/50/CH6/EX6.32/ex_6_32.sce new file mode 100755 index 000000000..eca7419c7 --- /dev/null +++ b/50/CH6/EX6.32/ex_6_32.sce @@ -0,0 +1,25 @@ +// example 6.32
+// caption: solving the IVP by numerov method
+// u''=(1+t^2)*u
+// u(0)=1, u'(0)=0 ,[0,1]
+// h=0.2,
+
+// expression for numerov method is
+//u(j+1)-2*u(j)+u(j-1)=(h^2/12)*(u''(j+1)+10*u''(j)+u''(j-1));
+
+// observing the IVP we can reduce the numerov method to
+//u(2)=2*u(1)-u(0)+(.2^2/12)*(1.16*u(2)+10.4*u(1)+1); for j=1
+// u(3)=2*u(2)-u(1)+(.2^2/12)*(1.36*u(3)+11.6*u(2)+1.04*u(1)); for j=2
+// u(4)=2*u(3)-u(2)+(.2^2/12)*(1.64*u(4)+13.6*u(3)+1.16*u(2)); for j=3
+// u(5)=2*u(4)-u(3)+(.2^2/12)*(2*u(5)+16.4*u(4)+1.36*u(3)); for j=4
+
+// from taylor series expansion we observe that
+u1=1.0202;u0=1;
+//u2-(.2^2/12)*(1.16*u2)=2*u1-u0+(.2^2/12)*(10.4*u1+1);
+u2=(1/.9961333)*2*u1-u0+(.2^2/12)*(10.4*u1+1)
+
+u3=(1/.995467)*(2.038667*u2-.996533*u1)
+
+u4=(1/.994533)*(2.045333*u3-.996133*u2)
+
+u5=(1/.993333)*(2.054667*u4-.995467*u3)
\ No newline at end of file diff --git a/50/CH6/EX6.4/ex_6_4.sce b/50/CH6/EX6.4/ex_6_4.sce new file mode 100755 index 000000000..f24a42d57 --- /dev/null +++ b/50/CH6/EX6.4/ex_6_4.sce @@ -0,0 +1,16 @@ +// example 6.4
+// solution to the given IVP
+
+disp('du/dt= A*u');
+// u=[u1 u2]';
+A=[-2 1;1 -20]; // given
+B=[1 0;0 1]; // identity matrix;
+
+
+
+
+
+[x,lam] = geigenvectors(A,B);
+
+// hence;
+disp('u=c1*%e^(lam(1)*t)*x(:,1)+c2*%e^-(lam(2)*t)*x(:,2)');
diff --git a/50/CH6/EX6.9/ex_6_9.sce b/50/CH6/EX6.9/ex_6_9.sce new file mode 100755 index 000000000..b3b8ddf45 --- /dev/null +++ b/50/CH6/EX6.9/ex_6_9.sce @@ -0,0 +1,14 @@ +// example 6.9
+// solve the IVP by euler method,
+// with h=0.2, 0.1, 0.05;
+// u'=f(t,u)
+// u'=-2tu^2
+deff('[z]=f(t,u)','z=-2*t*u^2');
+
+
+[u,t] = Euler1(1,0,1,.2,f) // h=0.2;
+
+[u,t] = Euler1(1,0,1,0.1,f) // h=0.1;
+
+
+[u,t] = Euler1(1,0,1,0.05,f) // h=0.05;
\ No newline at end of file diff --git a/50/CH7/EX7.1/ex_7_1.sce b/50/CH7/EX7.1/ex_7_1.sce new file mode 100755 index 000000000..d857f64d6 --- /dev/null +++ b/50/CH7/EX7.1/ex_7_1.sce @@ -0,0 +1,27 @@ +// example 7.1
+// solve by shooting method;
+
+// u''=u+1;
+// u(0)=0; u(1)=%e-1;
+
+// let -> U1(x)=du/dx;
+// U2(x)=d2u/dx2;
+
+// U(x)=[U1(x);U2(x)]
+
+// hence ;
+// dU/dx=f(x,U);
+
+
+
+deff('[w]=f(x,U)','w=[U(2); U(1)+1]')
+
+h=0.25;
+x=[0:h:1];
+ub=[0,%e-1];
+up=[0:1:10];
+
+
+[U] = shooting(ub,up,x,f);
+
+// the solution obtained would show the values of u and their derivatives at various x taken in regular intervals of h;
\ No newline at end of file diff --git a/50/CH7/EX7.11/ex_7_11.sce b/50/CH7/EX7.11/ex_7_11.sce new file mode 100755 index 000000000..2e785541b --- /dev/null +++ b/50/CH7/EX7.11/ex_7_11.sce @@ -0,0 +1,32 @@ +// example 7.11
+// solve the boundary value problem u''=u'+1;
+// u(0)=1; u(x=1)=2(%e-1); h=1/3;
+
+
+// we know; u''=(u(j-1)-2*u(j)+u(j+1))/h^2;
+// we know; u'=(u(j+1)-u(j-1))/2h;
+
+// 1) second order method;
+ x=0:1/3:1;
+
+ u=[u0 u1 u2 u3 ];
+// hence;
+disp('(u(j-1)-2*u(j)+u(j+1))/h^2=((u(j+1)-u(j-1))/2h)+1') // for j=1,2;
+
+
+disp('for j=1 (7/6)*u0-2*u1+(5/6)*u2=(1/9)')
+
+disp('for j=2 (7/6)*u1-2*u2+(5/6)*u3=(1/9)')
+
+
+// hence eliminating u1!
+// solving for u1,u2,
+u0=1;
+u3=2*(%e-1);
+u1=1.454869;
+u2=2.225019;
+
+disp(x);
+disp(u);
+
+
diff --git a/50/CH7/EX7.3/ex_7_3.sce b/50/CH7/EX7.3/ex_7_3.sce new file mode 100755 index 000000000..6c8f5c1c5 --- /dev/null +++ b/50/CH7/EX7.3/ex_7_3.sce @@ -0,0 +1,29 @@ +// example 7.3
+// solve by shooting method;
+
+// u''=2*u*u';
+// u(0)=0.5; u(1)=1;
+
+// let -> U1(x)=du/dx;
+// U2(x)=d2u/dx2;
+
+// U(x)=[U1(x);U2(x)]
+
+// hence ;
+// dU/dx=f(x,U);
+
+h=.25;
+
+ub=[.5,1];
+
+up=[0:.1:1];
+
+x=0:h:1;
+
+deff('[w]=f(x,U)','w=[U(2); 2*U(1)*U(2)]')
+
+
+
+[U] = shooting(ub,up,x,f);
+
+// the solution obtained would show the values of u in the first collumn and their corresponding derivatives in the second collumn ;
\ No newline at end of file diff --git a/50/CH7/EX7.4/ex_7_4.sce b/50/CH7/EX7.4/ex_7_4.sce new file mode 100755 index 000000000..2536351fd --- /dev/null +++ b/50/CH7/EX7.4/ex_7_4.sce @@ -0,0 +1,28 @@ +// example 7.4
+// solve by shooting method;
+
+// u''=2*u*u';
+// u(0)=0.5; u(1)=1;
+
+// let -> U1(x)=du/dx;
+// U2(x)=d2u/dx2;
+
+// U(x)=[U1(x);U2(x)]
+
+// hence ;
+// dU/dx=f(x,U);
+
+h=.25;
+
+ub=[.5,1];
+
+up=[0:.1:1];
+
+x=0:h:1;
+
+deff('[w]=f(x,U)','w=[U(2); 2*U(1)*U(2)]')
+
+
+[U] = shooting(ub,up,x,f);
+
+// the solution obtained would show the values of u in the first collumn and their corresponding derivatives in the second collumn ;
\ No newline at end of file diff --git a/50/CH7/EX7.5/ex_7_5.sce b/50/CH7/EX7.5/ex_7_5.sce new file mode 100755 index 000000000..8d6eed3f7 --- /dev/null +++ b/50/CH7/EX7.5/ex_7_5.sce @@ -0,0 +1,52 @@ +// example 7.5
+// solve the boundary value problem u''=u+x;
+// u(x=0)=u(0)=0; u(x=1)=u(4)=0; h=1/4;
+
+
+// we know; u''=(u(j-1)-2*u(j)+u(j+1))/h^2;
+
+// 1) second order method;
+ x=0:1/4:1;
+ u0=0;
+ u4=0;
+ u=[u0 u1 u2 u3 u4];
+// hence;
+disp('(u(j-1)-2*u(j)+u(j+1))/h^2=u(j)+x(j)') // for j=1,2,3;
+
+disp('for j=1 -16*u0+33*u1-16*u2=-.25')
+
+disp('for j=2 -16*u1+33*u2-16*u3=-.50')
+
+disp('for j=3 -16*u2+33*u3-16*u4=-.75')
+
+// hence solving for u1,u2,u3) ,
+u1=-.034885;
+u2=-.056326;
+u3=-.050037;
+
+disp(x);
+disp(u);
+
+// 2) numerov method;
+ x=0:1/4:1;
+ u0=0;
+ u4=0;
+ u=[u0 u1 u2 u3 u4];
+// since according to numerov method we get the following system of equations;
+disp('(191*u(j-1)-394*u(j)+191*u(j+1)=x(j-1)+10*x(j)+x(j+1)') // for j=1,2,3;
+
+disp('for j=1 191*u0-394*u1+191*u2=3')
+
+disp('for j=2 191*u1-394*u2+191*u3=6')
+
+disp('for j=3 191*u2-394*u3+191*u4=9')
+
+// hence solving for u1,u2,u3 ,
+u1=-.034885
+u2=-.056326
+u3=-.050037
+
+
+disp(x);
+disp(u);
+
diff --git a/50/CH7/EX7.6/ex_7_6.sce b/50/CH7/EX7.6/ex_7_6.sce new file mode 100755 index 000000000..04a28e398 --- /dev/null +++ b/50/CH7/EX7.6/ex_7_6.sce @@ -0,0 +1,32 @@ +// example 7.6
+// solve the boundary value problem u''=u*x;
+// u(0)+u'(0)=1; u(x=1)=0; h=1/3;
+
+
+// we know; u''=(u(j-1)-2*u(j)+u(j+1))/h^2;
+
+// 1) second order method;
+ x=0:1/3:1;
+
+ u3=1;
+ u=[u0 u1 u2 u3 ];
+// hence;
+disp('(u(j-1)-2*u(j)+u(j+1))/h^2=u(j)*x(j)') // for j=0,1,2,3;
+
+disp('for j=0 u1!-2*u0+u1=0') // u1!=u(-1)
+
+disp('for j=1 u0-2*u1+u2=(1/27)u1')
+
+disp('for j=2 u1-2*u2+u3=(2/27)u2')
+
+// we know; u'=(u(j+1)-u(j-1))/2h
+// hence eliminating u1!
+// solving for u0,u1,u2,u3 ,
+u0=-.9879518;
+u1=-.3253012;
+u2=-.3253012;
+
+disp(x);
+disp(u);
+
+
diff --git a/50/DEPENDENCIES/Euler1.sce b/50/DEPENDENCIES/Euler1.sce new file mode 100755 index 000000000..67fa05b6d --- /dev/null +++ b/50/DEPENDENCIES/Euler1.sce @@ -0,0 +1,26 @@ +function [u,t] = Euler1(u0,t0,tn,h,f)
+
+//Euler 1st order method solving ODE
+// du/dt = f(u,t), with initial
+//conditions u=u0 at t=t0. The
+//solution is obtained for t = [t0:h:tn]
+//and returned in u
+
+umaxAllowed = 1e+100;
+
+t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0;
+
+for j = 1:n-1
+ u(j+1) = u(j) + h*f(t(j),u(j));
+ if u(j+1) > umaxAllowed then
+ disp('Euler 1 - WARNING: underflow or overflow');
+ disp('Solution sought in the following range:');
+ disp([t0 h tn]);
+ disp('Solution evaluated in the following range:');
+ disp([t0 h t(j)]);
+ n = j; t = t(1,1:n); u = u(1,1:n);
+ break;
+ end;
+end;
+
+endfunction
diff --git a/50/DEPENDENCIES/F28.sce b/50/DEPENDENCIES/F28.sce new file mode 100755 index 000000000..e79b88348 --- /dev/null +++ b/50/DEPENDENCIES/F28.sce @@ -0,0 +1,8 @@ +function [F] = f28(x)
+
+//f2(x) = 0, with x = [x(1);x(2)]
+//represents a system of 2 non-linear equations
+F(1) = x(1)^2+x(1)*x(2) + x(2)^2 - 7;
+F(2) = x(1)^3+x(2)^3 -9;
+
+endfunction
diff --git a/50/DEPENDENCIES/LandU.sce b/50/DEPENDENCIES/LandU.sce new file mode 100755 index 000000000..fee5d09ea --- /dev/null +++ b/50/DEPENDENCIES/LandU.sce @@ -0,0 +1,12 @@ +function [U,L]=LandU(A,n)
+ U=A
+ L=eye(n,n)
+ for p=1:1:n-1
+ for i=p+1:1:n
+ m=A(i,p)/A(p,p);
+ L(i,p)=m;
+ A(i,:)=A(i,:)-m*A(p,:);
+ U=A;
+ end
+ end
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/NBDP.sci b/50/DEPENDENCIES/NBDP.sci new file mode 100755 index 000000000..535583f62 --- /dev/null +++ b/50/DEPENDENCIES/NBDP.sci @@ -0,0 +1,47 @@ +function [P]=NBDP(x,n,xL,f)
+//This function calculates a Newton Forward-Difference Polynomial of
+//order n, evaluated at x, using column vectors xL, f as the reference
+//table. The first value of xL and of f, represent, respectively,
+//xo and fo in the equation for the polynomial.
+[m,nc]=size(f)
+//check that it is indeed a column vector
+if (nc<>1) then
+ error('f is not a column vector.');
+ abort
+end;
+//check the difference order
+if (n >= m) then
+ disp(n,"n=");
+ disp(m,"m=");
+ error('n must be less than or equal to m-1');
+ abort
+end;
+//
+xo = xL(m,1);
+delx = mtlb_diff(xL);
+h = delx(1,1);
+s = (x-xo)/h;
+P = f(m,1);
+delf = f;
+disp(delf);
+for i = 1:n
+ delf = mtlb_diff(delf);
+ [m,nc] = size(delf);
+ disp(delf);
+ P = P + Binomial(s+i-1,i)*delf(m,1)
+end;
+endfunction
+
+function[C]=Binomial(s,i)
+ C = 1.0;
+ for k = 0:i-1
+ C = C*(s-k);
+ end;
+ C = C/factorial(i)
+endfunction
+function[fact]=factorial(nn)
+ fact = 1.0
+ for k = nn:-1:1
+ fact=fact*k
+ end;
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/NDDinterpol.sci b/50/DEPENDENCIES/NDDinterpol.sci new file mode 100755 index 000000000..5b49ebc9c --- /dev/null +++ b/50/DEPENDENCIES/NDDinterpol.sci @@ -0,0 +1,17 @@ +function P1=NDDinterpol (x0,x1,f0,f1)
+ x=poly(0,"x");
+ f01=(f1-f0)/(x1-x0);
+ P1=f0+(x-x0)*f01;
+endfunction
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file diff --git a/50/DEPENDENCIES/NDDinterpol2.sci b/50/DEPENDENCIES/NDDinterpol2.sci new file mode 100755 index 000000000..9618c8e9e --- /dev/null +++ b/50/DEPENDENCIES/NDDinterpol2.sci @@ -0,0 +1,19 @@ +function P2=NDDinterpol2 (x,f)
+ X=poly(0,"X");
+ f01=(f(2)-f(1))/(x(2)-x(1));
+ f13=(f(3)-f(2))/(x(3)-x(2));
+ f013=(f13-f01)/(x(3)-x(1));
+ P2=f(1)+(X-x(1))*f01+(X-x(1))*(X-x(2))*f013;
+endfunction
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file diff --git a/50/DEPENDENCIES/RK4.sci b/50/DEPENDENCIES/RK4.sci new file mode 100755 index 000000000..7539a4b15 --- /dev/null +++ b/50/DEPENDENCIES/RK4.sci @@ -0,0 +1,30 @@ +function [u,t] = RK4(u0,t0,tn,h,f)
+
+// RK4 method solving ODE
+// du/dt = f(u,t), with initial
+//conditions u=u0 at t=t0. The
+//solution is obtained for t = [t0:h:tn]
+//and returned in u
+
+umaxAllowed = 1e+100;
+
+t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0;
+
+for j = 1:n-1
+ k1=h*f(t(j),u(j));
+ k2=h*f(t(j)+h/2,u(j)+k1/2);
+ k3=h*f(t(j)+h/2,u(j)+k2/2);
+ k4=h*f(t(j)+h,u(j)+k3);
+ u(j+1) = u(j) + (1/6)*(k1+2*k2+2*k3+k4);
+ if u(j+1) > umaxAllowed then
+ disp('Euler 1 - WARNING: underflow or overflow');
+ disp('Solution sought in the following range:');
+ disp([t0 h tn]);
+ disp('Solution evaluated in the following range:');
+ disp([t0 h t(j)]);
+ n = j; t = t(1,1:n); u = u(1,1:n);
+ break;
+ end;
+end;
+
+endfunction
diff --git a/50/DEPENDENCIES/Vaitken.sce b/50/DEPENDENCIES/Vaitken.sce new file mode 100755 index 000000000..cf93e1b3a --- /dev/null +++ b/50/DEPENDENCIES/Vaitken.sce @@ -0,0 +1,9 @@ +// 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/50/DEPENDENCIES/Vbisection.sce b/50/DEPENDENCIES/Vbisection.sce new file mode 100755 index 000000000..cb23696bf --- /dev/null +++ b/50/DEPENDENCIES/Vbisection.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/50/DEPENDENCIES/Vbisection5.sce b/50/DEPENDENCIES/Vbisection5.sce new file mode 100755 index 000000000..da2e0f204 --- /dev/null +++ b/50/DEPENDENCIES/Vbisection5.sce @@ -0,0 +1,27 @@ +function x=bisection5(a,b,f)
+ N=5; //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/50/DEPENDENCIES/Vchebyshev.sce b/50/DEPENDENCIES/Vchebyshev.sce new file mode 100755 index 000000000..b883af5fa --- /dev/null +++ b/50/DEPENDENCIES/Vchebyshev.sce @@ -0,0 +1,18 @@ +function x=chebyshev(x,f,fp,fpp)
+ R=100;
+ PE=10^-5;
+ maxval=10^4;
+ if fp(x)==0 then disp("select another initial root x0");
+ break;
+ end
+ for n=1:1:R
+ x=x-f(x)/fp(x)-(1/2)*(f(x)/fp(x))^2 *(fpp(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/50/DEPENDENCIES/Vgausselim.sce b/50/DEPENDENCIES/Vgausselim.sce new file mode 100755 index 000000000..9488608a1 --- /dev/null +++ b/50/DEPENDENCIES/Vgausselim.sce @@ -0,0 +1,44 @@ +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;
+
+endfunction
+
diff --git a/50/DEPENDENCIES/Vgaussseidel.sce b/50/DEPENDENCIES/Vgaussseidel.sce new file mode 100755 index 000000000..ba227af87 --- /dev/null +++ b/50/DEPENDENCIES/Vgaussseidel.sce @@ -0,0 +1,24 @@ +function [X]=gaussseidel(A,n,N,X,b)
+ L=A;
+ U=A;
+ D=A;
+ for i=1:1:n
+ for j=1:1:n
+ if j>i then L(i,j)=0;
+ D(i,j)=0;
+ end
+ if i>j then U(i,j)=0;
+ D(i,j)=0;
+ end
+ if i==j then L(i,j)=0;
+ U(i,j)=0;
+ end
+ end
+
+ end
+ for k=1:1:N
+ X=(D+L)^-1*(-U*X+b);
+ disp(X)
+ end
+
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/Vnewton.sce b/50/DEPENDENCIES/Vnewton.sce new file mode 100755 index 000000000..f0fb337be --- /dev/null +++ b/50/DEPENDENCIES/Vnewton.sce @@ -0,0 +1,17 @@ +
+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/50/DEPENDENCIES/Vnewton4.sce b/50/DEPENDENCIES/Vnewton4.sce new file mode 100755 index 000000000..d35f210b4 --- /dev/null +++ b/50/DEPENDENCIES/Vnewton4.sce @@ -0,0 +1,17 @@ +function x=newton4(x,f,fp)
+ R=4;
+ PE=10^-15;
+ maxval=10^4;
+ for n=1:1:R
+ if fp(x)==0 then disp("select another initial root x0")
+ end
+ 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/50/DEPENDENCIES/Vregulafalsi.sce b/50/DEPENDENCIES/Vregulafalsi.sce new file mode 100755 index 000000000..7b6b0259a --- /dev/null +++ b/50/DEPENDENCIES/Vregulafalsi.sce @@ -0,0 +1,12 @@ +function [x]=regularfalsi(a,b,f)
+ N=100;
+ PE=10^-5;
+ for n=2:1:N
+ x=a-(a-b)*f(a)/(f(a)-f(b));
+ if abs(f(x))<=PE then break;
+ elseif (f(a)*f(x)<0) then b=x;
+ else a=x;
+ end
+ end
+ disp(n," no. of iterations =")
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/Vsecant.sce b/50/DEPENDENCIES/Vsecant.sce new file mode 100755 index 000000000..a44c89157 --- /dev/null +++ b/50/DEPENDENCIES/Vsecant.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/50/DEPENDENCIES/Vsim_eulercauchy.sce b/50/DEPENDENCIES/Vsim_eulercauchy.sce new file mode 100755 index 000000000..129472469 --- /dev/null +++ b/50/DEPENDENCIES/Vsim_eulercauchy.sce @@ -0,0 +1,24 @@ +function [u,v,t] = simeulercauchy(u0,v0,t0,tn,h,f1,f2)
+
+
+// du/dt = f1(t,u,v), dv/dt = f2(t,u,v) with initial
+//conditions u=u0,v=v0 at t=t0. The
+//solution is obtained for t = [t0:h:tn]
+//and returned in u,v
+
+
+umaxAllowed = 1e+100;
+
+t = [t0:h:tn]; u = zeros(t);v= zeros(t); n = length(u); u(1) = u0;v(1)=v0;
+
+for j = 1:n-1
+ k11=h*f1(t(j),u(j),v(j));
+ k21=h*f2(t(j),u(j),v(j));
+ k12=h*f1(t(j)+h,u(j)+k11,v(j)+k21);
+ k22=h*f2(t(j)+h,u(j)+k11,v(j)+k21);
+ u(j+1) = u(j) + (k11+k12)/2;
+ v(j+1) = v(j) + (k21+k22)/2;
+
+end;
+
+endfunction
diff --git a/50/DEPENDENCIES/adamsbashforth3.sci b/50/DEPENDENCIES/adamsbashforth3.sci new file mode 100755 index 000000000..dfc65961f --- /dev/null +++ b/50/DEPENDENCIES/adamsbashforth3.sci @@ -0,0 +1,23 @@ +function [u,t] = adamsbashforth3(u0,t0,tn,h,f)
+
+//adamsbashforth3 3rd order method solving ODE
+// du/dt = f(u,t), with initial
+//conditions u=u0 at t=t0. The
+//solution is obtained for t = [t0:h:tn]
+//and returned in u
+
+umaxAllowed = 1e+100;
+
+t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0;
+for j = 1:n-2
+if j<3 then
+ k1=h*f(t(j),u(j));
+ k2=h*f(t(j)+h,u(j)+k1);
+ u(j+1) = u(j) + (k2+k1)/2;
+end;
+
+if j>=2 then
+ u(j+2) = u(j+1) + (h/12)*(23*f(t(j+1),u(j+1))-16*f(t(j),u(j))+5*f(t(j-1),u(j-1)));
+end;
+end;
+endfunction
diff --git a/50/DEPENDENCIES/aitken.sce b/50/DEPENDENCIES/aitken.sce new file mode 100755 index 000000000..cf93e1b3a --- /dev/null +++ b/50/DEPENDENCIES/aitken.sce @@ -0,0 +1,9 @@ +// 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/50/DEPENDENCIES/aitkeninterpol.sci b/50/DEPENDENCIES/aitkeninterpol.sci new file mode 100755 index 000000000..2b6052d81 --- /dev/null +++ b/50/DEPENDENCIES/aitkeninterpol.sci @@ -0,0 +1,16 @@ +function P1=aitkeninterpol (x0,x1,f0,f1)
+ x=poly(0,"x");
+ P1=(1/(x1-x0))*det([f0 x0-x;f1 x1-x]);
+endfunction
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file diff --git a/50/DEPENDENCIES/back.sce b/50/DEPENDENCIES/back.sce new file mode 100755 index 000000000..d92199af3 --- /dev/null +++ b/50/DEPENDENCIES/back.sce @@ -0,0 +1,16 @@ +function [x] = back(U,Z)
+
+x=zeros(1,n);
+for i = n:-1:1
+ sumk=0
+ for j=i+1:n
+ sumk=sumk+U(i,j)*x(j);
+ end;
+ x(i)=(Z(i)-sumk)/U(i,i);
+end;
+
+
+
+
+endfunction
+
diff --git a/50/DEPENDENCIES/backeuler.sci b/50/DEPENDENCIES/backeuler.sci new file mode 100755 index 000000000..714f130ce --- /dev/null +++ b/50/DEPENDENCIES/backeuler.sci @@ -0,0 +1,22 @@ +function [u] = backeuler(u0,t0,tn,h,f)
+
+//backeuler 1st order method solving ODE
+// du/dt = f(u,t), with initial
+//conditions u=u0 at t=t0. The
+//solution is obtained for t = [t0:h:tn]
+//and returned in u
+
+umaxAllowed = 1e+100;
+
+ t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0;
+
+for j=1:n-1
+ u(j+1)=u(j);
+for i = 0:5
+ u(j+1) = u(j) + h*f(t(j+1),u(j+1));
+ i=i+1;
+end;
+end;
+
+
+endfunction
diff --git a/50/DEPENDENCIES/chebyshev.sce b/50/DEPENDENCIES/chebyshev.sce new file mode 100755 index 000000000..b883af5fa --- /dev/null +++ b/50/DEPENDENCIES/chebyshev.sce @@ -0,0 +1,18 @@ +function x=chebyshev(x,f,fp,fpp)
+ R=100;
+ PE=10^-5;
+ maxval=10^4;
+ if fp(x)==0 then disp("select another initial root x0");
+ break;
+ end
+ for n=1:1:R
+ x=x-f(x)/fp(x)-(1/2)*(f(x)/fp(x))^2 *(fpp(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/50/DEPENDENCIES/cholesky.sce b/50/DEPENDENCIES/cholesky.sce new file mode 100755 index 000000000..b3911d495 --- /dev/null +++ b/50/DEPENDENCIES/cholesky.sce @@ -0,0 +1,16 @@ +function L=cholesky (A,n)
+ L=zeros(n,n);
+ for k=1:1:n
+ S=0;
+ P=0;
+ for j=1:1:k-1
+ S=S+(L(k,j)^2);
+ P=P+L(i,j)*L(k,j)
+ end
+ L(k,k)=sqrt(A(k,k)-S);
+ for i=k+1:1:n
+ L(i,k)=(A(i,k)-P)/L(k,k);
+ end
+ end
+
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/comp_trapezoidal.sci b/50/DEPENDENCIES/comp_trapezoidal.sci new file mode 100755 index 000000000..f76a78571 --- /dev/null +++ b/50/DEPENDENCIES/comp_trapezoidal.sci @@ -0,0 +1,34 @@ +function I=comptrapezoidal(x,h,f)
+ //This function calculates the numerical integration of f(x)dx
+//between limits x(1) and x(n) using composite trapezoidal rule
+//Check that x and y have the same size (which must be an odd number)
+//Also, the values of x must be equally spaced with spacing h
+y=feval(x,f);
+[nrx,ncx]=size(x)
+[nrf,ncf]=size(y)
+if ((nrx<>1)|(nrf<>1)) then
+ error('x or f, or both, not column vector(s)');
+ abort;
+end;
+if ((ncx<>ncf)) then
+ error('x and f are not of the same length');
+ abort;
+end;
+//check that the size of the lists xL and f is odd
+if (modulo(ncx,2)==0) then
+ disp(ncx,"list size =")
+ error('list size must be an odd number');
+ abort
+end;
+n = ncx;
+
+I = f(x(1)) + f(x(n));
+for j = 2:n-1
+ if(modulo(j,2)==0) then
+ I = I + 2*f(x(j));
+ else
+ I = I + 2*f(x(j));
+ end;
+end;
+I = (h/2.0)*I
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/eigenvectors.sci b/50/DEPENDENCIES/eigenvectors.sci new file mode 100755 index 000000000..e36ffaaeb --- /dev/null +++ b/50/DEPENDENCIES/eigenvectors.sci @@ -0,0 +1,37 @@ +function [x,lam] = geigenvectors(A,B)
+
+//Calculates unit eigenvectors of matrix A
+//returning a matrix x whose columns are
+//the eigenvectors. The function also
+//returns the eigenvalues of the matrix.
+
+[nA,mA] = size(A);
+[nB,mB] = size(B);
+
+if (mA<>nA | mB<>nB) then
+ error('geigenvectors - matrix A or B not square');
+ abort;
+end;
+
+if nA<>nB then
+ error('geigenvectors - matrix A and B have different dimensions');
+ abort;
+end;
+
+lam = poly(0,'lam'); //Define variable "lam"
+chPoly = det(A-B*lam); //Characteristic polynomial
+lam = roots(chPoly)'; //Eigenvalues of matrix A
+
+x = []; n = nA;
+
+for k = 1:n
+ BB = A - lam(k)*B; //Characteristic matrix
+ CC = BB(1:n-1,1:n-1); //Coeff. matrix for reduced system
+ bb = -BB(1:n-1,n); //RHS vector for reduced system
+ y = CC\bb; //Solution for reduced system
+ y = [y;1]; //Complete eigenvector
+
+ x = [x y]; //Add eigenvector to matrix
+end;
+
+endfunction
diff --git a/50/DEPENDENCIES/eulermidpoint.sci b/50/DEPENDENCIES/eulermidpoint.sci new file mode 100755 index 000000000..f05453299 --- /dev/null +++ b/50/DEPENDENCIES/eulermidpoint.sci @@ -0,0 +1,26 @@ +function [u] = eulermidpoint(u0,t0,tn,h,f,fp)
+
+//midpoint 1st order method solving ODE
+// du/dt = f(u,t), with initial
+//conditions u=u0 at t=t0. The
+//solution is obtained for t = [t0:h:tn]
+//and returned in u
+
+umaxAllowed = 1e+100;
+
+t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0;
+u(2)=u(1)+h*f(t(1),u(1))+(h^2/2)*fp(t(1),u(1));
+for j = 2:n-1
+ u(j+1) = u(j-1) + 2*h*f(t(j),u(j));
+ if u(j+1) > umaxAllowed then
+ disp('Euler 1 - WARNING: underflow or overflow');
+ disp('Solution sought in the following range:');
+ disp([t0 h tn]);
+ disp('Solution evaluated in the following range:');
+ disp([t0 h t(j)]);
+ n = j; t = t(1,1:n); u = u(1,1:n);
+ break;
+ end;
+end;
+
+endfunction
diff --git a/50/DEPENDENCIES/fact.sci b/50/DEPENDENCIES/fact.sci new file mode 100755 index 000000000..c12441e95 --- /dev/null +++ b/50/DEPENDENCIES/fact.sci @@ -0,0 +1,6 @@ +function x=fact(n)
+ x=1;
+ for i=2:1:n
+ x=x*i;
+ end;
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/fore.sce b/50/DEPENDENCIES/fore.sce new file mode 100755 index 000000000..d8338b7cb --- /dev/null +++ b/50/DEPENDENCIES/fore.sce @@ -0,0 +1,12 @@ +function x=fore(L,b)
+
+for i = 1:1:n
+ sumk=0
+ for j=1:i-1
+ sumk=sumk+L(i,j)*x(j);
+ end;
+ x(i)=(b(i)-sumk)/L(i,i);
+end;
+
+endfunction
+
diff --git a/50/DEPENDENCIES/geigenvectors.sci b/50/DEPENDENCIES/geigenvectors.sci new file mode 100755 index 000000000..e36ffaaeb --- /dev/null +++ b/50/DEPENDENCIES/geigenvectors.sci @@ -0,0 +1,37 @@ +function [x,lam] = geigenvectors(A,B)
+
+//Calculates unit eigenvectors of matrix A
+//returning a matrix x whose columns are
+//the eigenvectors. The function also
+//returns the eigenvalues of the matrix.
+
+[nA,mA] = size(A);
+[nB,mB] = size(B);
+
+if (mA<>nA | mB<>nB) then
+ error('geigenvectors - matrix A or B not square');
+ abort;
+end;
+
+if nA<>nB then
+ error('geigenvectors - matrix A and B have different dimensions');
+ abort;
+end;
+
+lam = poly(0,'lam'); //Define variable "lam"
+chPoly = det(A-B*lam); //Characteristic polynomial
+lam = roots(chPoly)'; //Eigenvalues of matrix A
+
+x = []; n = nA;
+
+for k = 1:n
+ BB = A - lam(k)*B; //Characteristic matrix
+ CC = BB(1:n-1,1:n-1); //Coeff. matrix for reduced system
+ bb = -BB(1:n-1,n); //RHS vector for reduced system
+ y = CC\bb; //Solution for reduced system
+ y = [y;1]; //Complete eigenvector
+
+ x = [x y]; //Add eigenvector to matrix
+end;
+
+endfunction
diff --git a/50/DEPENDENCIES/generaliteration.sce b/50/DEPENDENCIES/generaliteration.sce new file mode 100755 index 000000000..3e89287e1 --- /dev/null +++ b/50/DEPENDENCIES/generaliteration.sce @@ -0,0 +1,21 @@ +
+function x=generaliteration(x,g,gp)
+ R=5;
+ PE=10^-8;
+ maxval=10^4;
+ k=gp(x);
+ if abs(k)>1 then error('function chosen does not converge')
+ abort;
+ end
+ for n=1:1:R
+ x=g(x);
+ disp(x);
+ if abs(g(x))<=PE then break
+ end
+ if (abs(g(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/50/DEPENDENCIES/generaliteration2.sce b/50/DEPENDENCIES/generaliteration2.sce new file mode 100755 index 000000000..29545493a --- /dev/null +++ b/50/DEPENDENCIES/generaliteration2.sce @@ -0,0 +1,22 @@ +
+function x=generaliteration2(x,g,gp)
+ R=2;
+ PE=10^-8;
+ maxval=10^4;
+ A=[0 0];
+ k=gp(x);
+ if abs(k)>1 then error('function chosen does not converge')
+ abort;
+ end
+ for n=1:1:R
+ x=g(x);
+ disp(x);
+ if abs(g(x))<=PE then break
+ end
+ if (abs(g(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/50/DEPENDENCIES/hermiteinterpol.sci b/50/DEPENDENCIES/hermiteinterpol.sci new file mode 100755 index 000000000..d9bc22148 --- /dev/null +++ b/50/DEPENDENCIES/hermiteinterpol.sci @@ -0,0 +1,37 @@ +function P= hermiteinterpol(x,f,fp)
+ X=poly(0,"X");
+ function f0=L0(X)
+
+ f0=(X-x(2))*(X-x(3))/((x(1)-x(2))*(x(1)-x(3)))
+ endfunction;
+ a0=[1-2*(X-x(1))*numdiff(L0,x(1))];
+L0=(X-x(2))*(X-x(3))/((x(1)-x(2))*(x(1)-x(3)));
+ A0=a0*L0*L0;
+ disp(A0)
+ B0=(X-x(1))*L0^2;
+
+ X=poly(0,"X");
+
+ function f1=L1(X)
+
+ f1=(X-x(1))*(X-x(3))/((x(2)-x(1))*(x(2)-x(3)))
+ endfunction;
+a1=[1-2*(X-x(2))*0];
+L1=(X-x(1))*(X-x(3))/((x(2)-x(1))*(x(2)-x(3)));
+A1=a1 *L1*L1;
+disp(A1)
+ B1=(X-x(2))*L1^2;
+ function f2=L2(X)
+
+ f2=(X-x(1))*(X-x(2))/((x(3)-x(1))*(x(3)-x(2)))
+ endfunction;
+a2=[1-2*(X-x(3))*numdiff(L2,x(3))];
+L2=(X-x(1))*(X-x(2))/((x(3)-x(1))*(x(3)-x(2)));
+A2=a2 *L2*L2;
+disp(A2)
+ B2=(X-x(3))*L2^2;
+
+
+
+ P=A0*f(1)+A1*f(2)+A2*f(3)+B0*fp(1)+B1*fp(2)+B2*fp(3);
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/heun.sci b/50/DEPENDENCIES/heun.sci new file mode 100755 index 000000000..f308280b0 --- /dev/null +++ b/50/DEPENDENCIES/heun.sci @@ -0,0 +1,28 @@ +function [u,t] = heun(u0,t0,tn,h,f)
+
+//heun method solving ODE
+// du/dt = f(u,t), with initial
+//conditions u=u0 at t=t0. The
+//solution is obtained for t = [t0:h:tn]
+//and returned in u
+
+umaxAllowed = 1e+100;
+
+t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0;
+
+for j = 1:n-1
+ k1=h*f(t(j),u(j));
+ k2=h*f(t(j)+h,u(j)+k1);
+ u(j+1) = u(j) + (k2+k1)/2;
+ if u(j+1) > umaxAllowed then
+ disp('Euler 1 - WARNING: underflow or overflow');
+ disp('Solution sought in the following range:');
+ disp([t0 h tn]);
+ disp('Solution evaluated in the following range:');
+ disp([t0 h t(j)]);
+ n = j; t = t(1,1:n); u = u(1,1:n);
+ break;
+ end;
+end;
+
+endfunction
diff --git a/50/DEPENDENCIES/iteratedinterpol.sci b/50/DEPENDENCIES/iteratedinterpol.sci new file mode 100755 index 000000000..2997bd6c5 --- /dev/null +++ b/50/DEPENDENCIES/iteratedinterpol.sci @@ -0,0 +1,18 @@ +function [L012,L02,L01]=iteratedinterpol (x,f,n)
+ X=poly(0,"X");
+ L01=(1/(x(2)-x(1)))*det([f(1) x(1)-X;f(2) x(2)-X]);
+ L02=(1/(x(3)-x(1)))*det([f(1) x(1)-X;f(3) x(3)-X]);
+ L012=(1/(x(3)-x(2)))*det([L01 x(2)-X;L02 x(3)-X]);
+
+endfunction
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file diff --git a/50/DEPENDENCIES/jacob28.sce b/50/DEPENDENCIES/jacob28.sce new file mode 100755 index 000000000..0553cd801 --- /dev/null +++ b/50/DEPENDENCIES/jacob28.sce @@ -0,0 +1,9 @@ +function [J] = jacob28(x)
+
+//Evaluates the Jacobian of a 2x2
+//system of non-linear equations
+
+J(1,1) = 2*x(1)+x(2); J(1,2) = x(1)+2*x(2);
+J(2,1) = 3*x(1)^2; J(2,2) =3* x(2)^2;
+endfunction
+
diff --git a/50/DEPENDENCIES/jacobianmat.sci b/50/DEPENDENCIES/jacobianmat.sci new file mode 100755 index 000000000..1a6d8ba75 --- /dev/null +++ b/50/DEPENDENCIES/jacobianmat.sci @@ -0,0 +1,8 @@ +function J= jacobianmat (f1,f2,h,k)
+ J=zeros(2,2);
+J(1,1)=(f1(1+h,1)-f1(1,1))/2*h;
+
+J(1,2)=(f1(1,1+k)-f1(1,1))/2*k;
+J(2,1)=(f2(1+h,1)-f2(1,1))/2*h;
+J(2,2)=(f2(1,1+k)-f2(1,1))/2*k;
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/jacobiiteration.sce b/50/DEPENDENCIES/jacobiiteration.sce new file mode 100755 index 000000000..2cdc4e224 --- /dev/null +++ b/50/DEPENDENCIES/jacobiiteration.sce @@ -0,0 +1,23 @@ +function [X]=jacobiiteration(A,n,N,X,b)
+ L=A;
+ U=A;
+ D=A;
+ for i=1:1:n
+ for j=1:1:n
+ if j>i then L(i,j)=0;
+ D(i,j)=0;
+ end
+ if i>j then U(i,j)=0;
+ D(i,j)=0;
+ end
+ if i==j then L(i,j)=0;
+ U(i,j)=0;
+ end
+ end
+
+ end
+ for k=1:1:N
+ X=-D^-1*(L+U)*X+D^-1*(b);
+ end
+
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/jordan.sce b/50/DEPENDENCIES/jordan.sce new file mode 100755 index 000000000..22ba3c1fa --- /dev/null +++ b/50/DEPENDENCIES/jordan.sce @@ -0,0 +1,20 @@ +function [M]=jorden(A,b)
+ M=[A b];
+ [ra,ca]=size(A);
+ [rb,cb]=size(b);
+ n=ra;
+ for p=1:1:n
+ for k=(p+1):1:n
+ if abs(M(k,p))>abs(M(p,p)) then
+ M({p,k},:)=M({k,p},:);
+ end
+ end
+ M(p,:)=M(p,:)/M(p,p);
+ for i=1:1:p-1
+ M(i,:)=M(i,:)-M(p,:)*(M(i,p)/M(p,p));
+ end
+ for i=p+1:1:n
+ M(i,:)=M(i,:)-M(p,:)*(M(i,p)/M(p,p));
+ end
+ end
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/lagrangefundamentalpoly.sci b/50/DEPENDENCIES/lagrangefundamentalpoly.sci new file mode 100755 index 000000000..b55a42ff6 --- /dev/null +++ b/50/DEPENDENCIES/lagrangefundamentalpoly.sci @@ -0,0 +1,30 @@ +function P2=lagrangefundamentalpoly(x,f,n)
+ [nrx,ncx]=size(x)
+ [nrf,ncf]=size(f)
+if ((nrx<>1)|(nrf<>1)) then
+ error('x or f, or both, not column vector(s)');
+ abort;
+end;
+if ((ncx<>ncf)) then
+ error('x and f are not of the same length');
+ abort;
+end;
+
+ X=poly(0,"X");
+L=zeros(n);
+
+P2=0;
+for i=1:n+1
+ L(i)=1;
+ for j=1:n+1
+ if i~=j then
+ L(i)=L(i)*(X-x(j))/(x(i)-x(j))
+ end;
+end;
+P2=P2+L(i)*f(i);
+end;
+
+
+endfunction
+
+
diff --git a/50/DEPENDENCIES/legrangeinterpol.sci b/50/DEPENDENCIES/legrangeinterpol.sci new file mode 100755 index 000000000..cbaff9113 --- /dev/null +++ b/50/DEPENDENCIES/legrangeinterpol.sci @@ -0,0 +1,18 @@ +function P1=legrangeinterpol (x0,x1,f0,f1)
+ x=poly(0,"x");
+ L0=(x-x1)/(x0-x1);
+ L1=(x-x0)/(x1-x0);
+ P1=L0*f0+L1*f1;
+endfunction
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file diff --git a/50/DEPENDENCIES/linearinterpol.sci b/50/DEPENDENCIES/linearinterpol.sci new file mode 100755 index 000000000..a6c10cbc5 --- /dev/null +++ b/50/DEPENDENCIES/linearinterpol.sci @@ -0,0 +1,3 @@ +function fp=linearinterpol(xL,f)
+ fp=(f(2)-f(1))/(xL(2)-xL(1));
+endfunction;
\ No newline at end of file diff --git a/50/DEPENDENCIES/modified_newton.sce b/50/DEPENDENCIES/modified_newton.sce new file mode 100755 index 000000000..29d60f47b --- /dev/null +++ b/50/DEPENDENCIES/modified_newton.sce @@ -0,0 +1,17 @@ +
+function x=modified_newton(x,f,fp)
+ R=100;
+ PE=10^-8;
+ maxval=10^4;
+
+ for n=1:1:R
+ x=x-m*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/50/DEPENDENCIES/modifiedeuler.sci b/50/DEPENDENCIES/modifiedeuler.sci new file mode 100755 index 000000000..a515eff0e --- /dev/null +++ b/50/DEPENDENCIES/modifiedeuler.sci @@ -0,0 +1,28 @@ +function [u,t] = modifiedeuler(u0,t0,tn,h,f)
+
+//modifiedeuler 1st order method solving ODE
+// du/dt = f(u,t), with initial
+//conditions u=u0 at t=t0. The
+//solution is obtained for t = [t0:h:tn]
+//and returned in u
+
+umaxAllowed = 1e+100;
+
+t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0;
+
+for j = 1:n-1
+ k1=h*f(t(j),u(j));
+ k2=h*f(t(j)+h/2,u(j)+k1/2);
+ u(j+1) = u(j) + k2;
+ if u(j+1) > umaxAllowed then
+ disp('Euler 1 - WARNING: underflow or overflow');
+ disp('Solution sought in the following range:');
+ disp([t0 h tn]);
+ disp('Solution evaluated in the following range:');
+ disp([t0 h t(j)]);
+ n = j; t = t(1,1:n); u = u(1,1:n);
+ break;
+ end;
+end;
+
+endfunction
diff --git a/50/DEPENDENCIES/muller3.sce b/50/DEPENDENCIES/muller3.sce new file mode 100755 index 000000000..75007237c --- /dev/null +++ b/50/DEPENDENCIES/muller3.sce @@ -0,0 +1,32 @@ +function x=muller3(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/50/DEPENDENCIES/muller5.sce b/50/DEPENDENCIES/muller5.sce new file mode 100755 index 000000000..e2cf73b8a --- /dev/null +++ b/50/DEPENDENCIES/muller5.sce @@ -0,0 +1,32 @@ +function x=muller5(x0,x1,x2,f)
+ R=5;
+ 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/50/DEPENDENCIES/multipoint_iteration31.sce b/50/DEPENDENCIES/multipoint_iteration31.sce new file mode 100755 index 000000000..5c016939d --- /dev/null +++ b/50/DEPENDENCIES/multipoint_iteration31.sce @@ -0,0 +1,14 @@ +function x=multipoint_iteration31(x,f,fp,R)
+ R=3;
+ PE=10^-5;
+ maxval=10^4;
+ for n=1:1:R
+ x=x-f(x)/fp(x-(1/2)*(f(x)/fp(x)));
+ if abs(f(x))<=PE then break;
+ end
+ if (abs(f(x))>maxval) then error('Solution diverges');
+ break
+ end
+ end
+ disp(n," no. of iterations =")
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/multipoint_iteration33.sce b/50/DEPENDENCIES/multipoint_iteration33.sce new file mode 100755 index 000000000..bd4bc3308 --- /dev/null +++ b/50/DEPENDENCIES/multipoint_iteration33.sce @@ -0,0 +1,14 @@ +function x=multipoint_iteration33(x,f,fp,R)
+ R=3;
+ PE=10^-5;
+ maxval=10^4;
+ for n=1:1:R
+ x=x-f(x)/fp(x)-f(x-(f(x)/fp(x)))/fp(x);
+ if abs(f(x))<=PE then break;
+ end
+ if (abs(f(x))>maxval) then error('Solution diverges');
+ break
+ end
+ end
+ disp(n," no. of iterations =")
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/newton2.sce b/50/DEPENDENCIES/newton2.sce new file mode 100755 index 000000000..12714759c --- /dev/null +++ b/50/DEPENDENCIES/newton2.sce @@ -0,0 +1,21 @@ +function [x]=newton2(x0,f,fp)
+//newton-raphson algorithm
+N = 100; eps = 1.e-5; // define max. no. iterations and error
+maxval = 10000.0; // define value for divergence
+xx = x0;
+while (N>0)
+ xn = xx-f(xx)/fp(xx);
+ if(abs(f(xn))<=eps) then x=xn;
+ disp(100-N);
+ return(x);
+ end;
+ if (abs(f(xx))>maxval) then disp(100-N);
+ error('Solution diverges');
+ abort;
+ end;
+ N = N - 1;
+ xx = xn;
+end;
+error('No convergence');
+abort;
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/newton63.sce b/50/DEPENDENCIES/newton63.sce new file mode 100755 index 000000000..186634d2e --- /dev/null +++ b/50/DEPENDENCIES/newton63.sce @@ -0,0 +1,17 @@ +
+function x=newton63(x,f,fp,fpp)
+ R=100;
+ PE=10^-15;
+ maxval=10^4;
+
+ for n=1:1:R
+ x=x-(f(x)*fp(x))/(fp(x)^2-f(x)*fpp(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/50/DEPENDENCIES/newtonNLE.sce b/50/DEPENDENCIES/newtonNLE.sce new file mode 100755 index 000000000..cb0e9d5c9 --- /dev/null +++ b/50/DEPENDENCIES/newtonNLE.sce @@ -0,0 +1,31 @@ +function [xn] =newtonNLE(x0,f,J)
+
+//Newton-Raphson method applied to a
+//system of linear equations f(x) = 0,
+//given the jacobian function J, with
+//J is the jacobian matrix of the functions.
+//x = [x1;x2;...;xn], f = [f1;f2;...;fn]
+//x0 is an initial guess of the solution
+
+N = 3; //define max. number of iterations
+PE= 10^-15; //define tolerance
+maxval = 10000.0; //define value for divergence
+xx = x0; //load initial guess
+for n=1:1:N
+ JJ = J(xx);
+ if(abs(det(JJ))<PE) then
+ error('newtonm - Jacobian is singular - try new x0');
+ abort;
+ end;
+ xn = xx - inv(JJ)*f(xx);
+ if(abs(f(xn))<=PE) then break;
+ else xx=xn;
+ end;
+ if (abs(f(xx))>=maxval) then
+ error('Solution diverges');
+ abort;
+ end;
+
+end;
+disp(n," no. of iterations =");
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/newtonrap.sce b/50/DEPENDENCIES/newtonrap.sce new file mode 100755 index 000000000..63fc2abc8 --- /dev/null +++ b/50/DEPENDENCIES/newtonrap.sce @@ -0,0 +1,18 @@ +
+function x=newton(x,f,fp)
+ R=5;
+ PE=10^-15;
+ 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/50/DEPENDENCIES/pivotgausselim.sci b/50/DEPENDENCIES/pivotgausselim.sci new file mode 100755 index 000000000..3fbd3402a --- /dev/null +++ b/50/DEPENDENCIES/pivotgausselim.sci @@ -0,0 +1,27 @@ +function [x]=pivotgausselim(A,b)
+ M=[A b];
+ [ra,ca]=size(A);
+ [rb,cb]=size(b);
+ n=ra;
+ for p=1:1:n
+ for k=(p+1):1:n
+ if abs(M(k,p))>abs(M(p,p)) then
+ M({p,k},:)=M({k,p},:);
+ end
+ end
+ for i=p+1:1:n
+ m(i,p)=M(i,p)/M(p,p);
+ M(i,:)=M(i,:)-M(p,:)*m(i,p);
+
+ end
+ end
+ a=M(1:n,1:n);
+ b=M(:,n+1);
+ for i = n:-1:1
+ sumj=0
+ for j=n:-1:i+1
+ sumj = sumj + a(i,j)*x(j);
+ end;
+ x(i)=(b(i)-sumj)/a(i,i);
+ end
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/quadraticapprox.sci b/50/DEPENDENCIES/quadraticapprox.sci new file mode 100755 index 000000000..ad49267e9 --- /dev/null +++ b/50/DEPENDENCIES/quadraticapprox.sci @@ -0,0 +1,22 @@ +function [P]=quadraticapprox(x,f)
+
+n=length(x);m=length(f);
+if m<>n then
+ error('linreg - Vectors x and f are not of the same length.');
+ abort;
+end;
+s1=0;
+s2=0;
+for i=1:n
+ s1=s1+x(i)*f(i);
+ s2=s2+x(i)^2*f(i);
+end
+c0=det([sum(f) sum(x) sum(x^2);s1 sum(x^2) sum(x^3); s2 sum(x^3) sum(x^4)])/det([n sum(x) sum(x^2);sum(x) sum(x^2) sum(x^3); sum(x^2) sum(x^3) sum(x^4)]);
+
+c1=det([n sum(f) sum(x^2);sum(x) s1 sum(x^3); sum(x^2) s2 sum(x^4)])/det([n sum(x) sum(x^2);sum(x) sum(x^2) sum(x^3); sum(x^2) sum(x^3) sum(x^4)]);
+
+c2=det([n sum(x) sum(f);sum(x) sum(x^2) s1; sum(x^2) sum(x^3) s2])/det([n sum(x) sum(x^2);sum(x) sum(x^2) sum(x^3); sum(x^2) sum(x^3) sum(x^4)]);
+
+X=poly(0,"X");
+P=c2*X^2+c1*X+c0;
+endfunction
diff --git a/50/DEPENDENCIES/regulafalsi.sce b/50/DEPENDENCIES/regulafalsi.sce new file mode 100755 index 000000000..1b93787bf --- /dev/null +++ b/50/DEPENDENCIES/regulafalsi.sce @@ -0,0 +1,12 @@ +function [x]=regulafalsi(a,b,f)
+ N=100;
+ PE=10^-5;
+ for n=2:1:N
+ x=a-(a-b)*f(a)/(f(a)-f(b));
+ if abs(f(x))<=PE then break;
+ elseif (f(a)*f(x)<0) then b=x;
+ else a=x;
+ end
+ end
+ disp(n," no. of iterations =")
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/regulafalsi4.sce b/50/DEPENDENCIES/regulafalsi4.sce new file mode 100755 index 000000000..5028fbdb0 --- /dev/null +++ b/50/DEPENDENCIES/regulafalsi4.sce @@ -0,0 +1,12 @@ +function [x]=regulafalsi4(a,b,f)
+ N=100;
+ PE=10^-5;
+ for n=2:1:N
+ x=a-(a-b)*f(a)/(f(a)-f(b));
+ if abs(f(x))<=PE then break;
+ elseif (f(a)*f(x)<0) then b=x;
+ else a=x;
+ end
+ end
+ disp(n," no. of iterations =")
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/secant4.sce b/50/DEPENDENCIES/secant4.sce new file mode 100755 index 000000000..4cca19886 --- /dev/null +++ b/50/DEPENDENCIES/secant4.sce @@ -0,0 +1,12 @@ +function [x]=secant4(a,b,f)
+ N=4; // 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/50/DEPENDENCIES/secant64.sce b/50/DEPENDENCIES/secant64.sce new file mode 100755 index 000000000..b65a9f3d5 --- /dev/null +++ b/50/DEPENDENCIES/secant64.sce @@ -0,0 +1,12 @@ +function [x]=secant64(a,b,f,fp)
+ N=100; // define max. no. iterations to be performed
+ PE=10^-15 // define tolerance for convergence
+ for n=1:1:N // initiating for loop
+ x=(b*f(a)*fp(b)-a*f(b)*fp(a))/(f(a)*fp(b)-f(b)*fp(a));
+ 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/50/DEPENDENCIES/secant65.sce b/50/DEPENDENCIES/secant65.sce new file mode 100755 index 000000000..5d7d4e52c --- /dev/null +++ b/50/DEPENDENCIES/secant65.sce @@ -0,0 +1,13 @@ +function [x]=secant65(a,b,f)
+ deff('[y]=g(x)','y=-f(x)^2/(f(x-f(x))-f(x))');
+ N=4; // define max. no. iterations to be performed
+ PE=10^-15 // define tolerance for convergence
+ for n=1:1:N // initiating for loop
+ x=a-(b-a)*g(a)/(g(b)-g(a));
+ 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/50/DEPENDENCIES/shooting.sci b/50/DEPENDENCIES/shooting.sci new file mode 100755 index 000000000..a9bf9bee6 --- /dev/null +++ b/50/DEPENDENCIES/shooting.sci @@ -0,0 +1,31 @@ +function [U] = shooting(ub,up,x,f)
+
+//Shooting method for a second order
+//boundary value problem
+//ub = [u0 u1] -> boundary conditions
+//x = a vector showing the range of x
+//f = function defining ODE, i.e.,
+// du/dx = f(x,u), u = [u(1);u(2)].
+//up = vector with range of du/dx at x=x0
+//xuTable = table for interpolating derivatives
+//uderiv = derivative boundary condition
+
+n = length(up);
+m = length(x);
+y1 = zeros(up);
+
+for j = 1:n
+ u0 = [ub(1);up(j)];
+ uu = ode(u0,x(1),x,f);
+ u1(j) = uu(1,m);
+end;
+
+xuTable = [u1';up];
+uderiv = interpln(xuTable,ub(2));
+u0 = [ub(1);uderiv];
+u = ode(u0,x(1),x,f);
+U=u';
+
+endfunction
+
+
diff --git a/50/DEPENDENCIES/simRK4.sci b/50/DEPENDENCIES/simRK4.sci new file mode 100755 index 000000000..d7dc4a913 --- /dev/null +++ b/50/DEPENDENCIES/simRK4.sci @@ -0,0 +1,27 @@ +function [u,v,t] = simRK4(u0,v0,t0,tn,h,f1,f2)
+
+// RK4 method solving simultanious ODE
+// du/dt = f1(t,u,v), dv/dt = f2(t,u,v) with initial
+//conditions u=u0,v=v0 at t=t0. The
+//solution is obtained for t = [t0:h:tn]
+//and returned in u,v
+
+umaxAllowed = 1e+100;
+
+t = [t0:h:tn]; u = zeros(t);v=zeros(t) ;n = length(u); u(1) = u0;v(1)=v0
+
+for j = 1:n-1
+ k11=h*f1(t(j),u(j),v(j));
+ k21=h*f2(t(j),u(j),v(j));
+ k12=h*f1(t(j)+h/2,u(j)+k11/2,v(j)+k21/2);
+ k22=h*f2(t(j)+h/2,u(j)+k11/2,v(j)+k21/2);
+ k13=h*f1(t(j)+h/2,u(j)+k12/2,v(j)+k22/2);
+ k23=h*f2(t(j)+h/2,u(j)+k12/2,v(j)+k22/2);
+ k14=h*f1(t(j)+h,u(j)+k13,v(j)+k23);
+ k24=h*f2(t(j)+h,u(j)+k13,v(j)+k23);
+ u(j+1) = u(j) + (1/6)*(k11+2*k12+2*k13+k14);
+ v(j+1) = v(j) + (1/6)*(k21+2*k22+2*k23+k24);
+
+end;
+
+endfunction
diff --git a/50/DEPENDENCIES/sim_eulercauchy.sce b/50/DEPENDENCIES/sim_eulercauchy.sce new file mode 100755 index 000000000..b16af4509 --- /dev/null +++ b/50/DEPENDENCIES/sim_eulercauchy.sce @@ -0,0 +1,24 @@ +function [u,v,t] = simeulercauchy(u0,v0,t0,tn,h,f1,f2)
+
+
+// du/dt = f1(t,u,v), dv/dt = f2(t,u,v) with initial
+//conditions u=u0,v=v0 at t=t0. The
+//solution is obtained for t = [t0:h:tn]
+//and returned in u,v
+
+
+umaxAllowed = 1e+100;
+
+t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0;
+
+for j = 1:n-1
+ k11=h*f1(t(j),u(j),v(j));
+ k21=h*f2(t(j),u(j),v(j));
+ k12=h*f1(t(j)+h,u(j)+k11,v(j)+k21);
+ k22=h*f2(t(j)+h,u(j)+k11,v(j)+k21);
+ u(j+1) = u(j) + (k11+k12)/2;
+ v(j+1) = v(j) + (k21+k22)/2;
+
+end;
+
+endfunction
diff --git a/50/DEPENDENCIES/simpson.sci b/50/DEPENDENCIES/simpson.sci new file mode 100755 index 000000000..da9992fd8 --- /dev/null +++ b/50/DEPENDENCIES/simpson.sci @@ -0,0 +1,3 @@ +function I=simpson(a,b,f)
+ I=((b-a)/6)*(f(a)+4*f((a+b)/2)+f(b));
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/simpson13.sci b/50/DEPENDENCIES/simpson13.sci new file mode 100755 index 000000000..3e3c6b326 --- /dev/null +++ b/50/DEPENDENCIES/simpson13.sci @@ -0,0 +1,34 @@ +function [I] = simpson13(x,h,f)
+//This function calculates the numerical integration of f(x)dx
+//between limits x(1) and x(n) using Simpson's 1/3 rule
+//Check that x and y have the same size (which must be an odd number)
+//Also, the values of x must be equally spaced with spacing h
+y=feval(x,f);
+[nrx,ncx]=size(x)
+[nrf,ncf]=size(y)
+if ((nrx<>1)|(nrf<>1)) then
+ error('x or f, or both, not column vector(s)');
+ abort;
+end;
+if ((ncx<>ncf)) then
+ error('x and f are not of the same length');
+ abort;
+end;
+//check that the size of the lists xL and f is odd
+if (modulo(ncx,2)==0) then
+ disp(ncx,"list size =")
+ error('list size must be an odd number');
+ abort
+end;
+n = ncx;
+
+I = f(x(1)) + f(x(n));
+for j = 2:n-1
+ if(modulo(j,2)==0) then
+ I = I + 4*f(x(j));
+ else
+ I = I + 2*f(x(j));
+ end;
+end;
+I = (h/3.0)*I
+endfunction
diff --git a/50/DEPENDENCIES/simpson38.sci b/50/DEPENDENCIES/simpson38.sci new file mode 100755 index 000000000..ca057dfa7 --- /dev/null +++ b/50/DEPENDENCIES/simpson38.sci @@ -0,0 +1,37 @@ +function [I] = simpson38(x,f)
+//This function calculates the numerical integration of f(x)dx
+//between limits x(1) and x(n) using Simpson's 3/8 rule
+//Check that x and f have the same size (which must be of the form 3*i+1,
+//where i is an integer number)
+//Also, the values of x must be equally spaced with spacing h
+
+y=feval(x,f);
+[nrx,ncx]=size(x)
+[nrf,ncf]=size(y)
+if ((nrx<>1)|(nrf<>1)) then
+ error('x or f, or both, not column vector(s)');
+ abort;
+end;
+if ((ncx<>ncf)) then
+ error('x and f are not of the same length');
+ abort;
+end;
+//check that the size of the lists xL and f is odd
+if (modulo(ncx-1,3)<>0) then
+ disp(ncx,"list size =")
+ error('list size must be of the form 3*i+1, where i=integer');
+ abort
+end;
+n = ncx;
+xdiff = mtlb_diff(x);
+h = xdiff(1,1);
+I = f(x(1)) + f(x(n));
+for j = 2:n-1
+ if(modulo(j-1,3)==0) then
+ I = I + 2*f(x(j));
+ else
+ I = I + 3*f(x(j));
+ end;
+end;
+I = (3.0/8.0)*h*I
+endfunction
diff --git a/50/DEPENDENCIES/straightlineapprox.sce b/50/DEPENDENCIES/straightlineapprox.sce new file mode 100755 index 000000000..9b8e50e19 --- /dev/null +++ b/50/DEPENDENCIES/straightlineapprox.sce @@ -0,0 +1,16 @@ +function [P]=straightlineapprox(x,f)
+
+n=length(x);m=length(f);
+if m<>n then
+ error('linreg - Vectors x and f are not of the same length.');
+ abort;
+end;
+s=0;
+for i=1:n
+ s=s+x(i)*f(i);
+end
+c0=det([sum(f) sum(x); s sum(x^2)])/det([n sum(x);sum(x) sum(x^2)]);
+c1=det([ n sum(f); sum(x) s])/det([n sum(x);sum(x) sum(x^2)]);
+X=poly(0,"X");
+P=c1*X+c0;
+endfunction
diff --git a/50/DEPENDENCIES/taylor.sci b/50/DEPENDENCIES/taylor.sci new file mode 100755 index 000000000..ea5fdc776 --- /dev/null +++ b/50/DEPENDENCIES/taylor.sci @@ -0,0 +1,3 @@ +function u=taylor(t)
+ u=(t^1*U1)/fact(1)+(t^2*U2)/fact(2)+(t^3*U3)/fact(3)+(t^4*U4)/fact(4)+(t^5*U5)/fact(5)+(t^6*U6)/fact(6)+(t^7*U7)/fact(7)+(t^8*U8)/fact(8)+(t^9*U9)/fact(9)+(t^10*U10)/fact(10)+(t^11*U11)/fact(11)
+endfunction
\ No newline at end of file diff --git a/50/DEPENDENCIES/trapezoidal.sci b/50/DEPENDENCIES/trapezoidal.sci new file mode 100755 index 000000000..0ec3e7e49 --- /dev/null +++ b/50/DEPENDENCIES/trapezoidal.sci @@ -0,0 +1,7 @@ +// solves the definite integral by the trapezoidal rule,
+// given the limits a,b and the function f,
+// returns the integral value I
+
+function I =trapezoidal(a,b,f)
+ I=((b-a)/2)*(f(a)+f(b));
+endfunction
\ No newline at end of file |
