initial commit / add all books

18 files changed, 392 insertions, 0 deletions

diff --git a/964/CH25/EX25.1/25_1.sce b/964/CH25/EX25.1/25_1.sce
new file mode 100755
index 000000000..3b6647efd
--- /dev/null
+++ b/964/CH25/EX25.1/25_1.sce
@@ -0,0 +1,20 @@
+//clc()

+//dy/dx = -2*x^3 + 12*x^2 - 20*x + 8.5

+//therefore, y = -0.5*x^4 + 4*x^3 - 10*x^2 + 8.5 + c

+x1 = 0;

+y1 = 1;

+h = 0.5;

+c =-(-0.5*x1^4 + 4*x1^3 - 10*x1^2 + 8.5*x1 - y1);

+x = 0:0.5:4;

+disp(x,"x = ")

+y = -0.5*x^4 + 4*x^3 - 10*x^2 + 8.5*x + c;

+disp(y,"true values of y = ")

+fxy = -2*x^3 + 12*x^2 - 20*x + 8.5;

+y2(1) = y(1);

+e(1) = (y(1) - y2(1)) * 100 / y(1);

+for i = 2:9

+    y2(i) = y2(i-1) + fxy(i-1)*h;

+    e(i) = (y(i) - y2(i))*100/y(i);

+end

+disp(y2,"y by euler method =")

+disp(e,"error =")

diff --git a/964/CH25/EX25.10/25_10.sce b/964/CH25/EX25.10/25_10.sce
new file mode 100755
index 000000000..95feeae4c
--- /dev/null
+++ b/964/CH25/EX25.10/25_10.sce
@@ -0,0 +1,30 @@
+//clc()

+//dy1/dx = -0.5*y1

+//dy2/dx = 4 - 0.3*y2 - 0.1*y1

+x1 = 0;

+h =0.5;

+y1(1) = 4;

+y2(1) = 6;

+x = 0:h:2;

+for i = 1:5

+    k11(i) = -0.5*y1(i);

+    k12(i) = 4 - 0.3*y2(i) - 0.1*y1(i);

+    y11 = y1(i) + k11(i) * h/2;

+    y21 = y2(i) + k12(i) * h/2;

+    k21(i) = -0.5*y11;

+    k22(i) = 4 - 0.3*y21 - 0.1*y11;

+    y11 = y1(i) + k21(i) * h/2;

+    y21 = y2(i) + k22(i) * h/2;

+    k31(i) = -0.5*y11;

+    k32(i) = 4 - 0.3*y21 - 0.1*y11;

+    y11 = y1(i) + k31(i) * h;

+    y21 = y2(i) + k32(i) * h;

+    k41(i) = -0.5*y11;

+    k42(i) = 4 - 0.3*y21 - 0.1*y11;

+    y1(i+1) = y1(i) + ( k11(i) + 2*k21(i) + 2*k31(i) + k41(i) )*h / 6;

+    y2(i+1) = y2(i) + ( k12(i) + 2*k22(i) + 2*k32(i) + k42(i) ) *h/ 6;

+end

+disp("using fourth order RK method")

+disp(x,"x =")

+disp(y1(1:5),"y1 = ")

+disp(y2(1:5),"y2 = ")

diff --git a/964/CH25/EX25.11/25_11.sce b/964/CH25/EX25.11/25_11.sce
new file mode 100755
index 000000000..ec5a7c7c2
--- /dev/null
+++ b/964/CH25/EX25.11/25_11.sce
@@ -0,0 +1,57 @@
+clc;

+clear;

+function yp=f(x,y)

+    yp=[y(2);-16.1*y(1)];

+endfunction

+x=0:0.1:4

+y0=[0.1 0];

+sol=ode(y0,0,x,f);

+count=1;

+disp(sol)

+for i=1:2:81

+    a(count)=sol(i);

+    b(count)=sol(i+1);

+    count=count+1;

+end

+plot(x,a)

+plot(x,b,".-")

+h1=legend(["y1,y3","y2,y4"])

+xtitle("y vs x","x","y")

+function yp=g(x,y)

+    yp=[y(2);-16.1*sin(y(1))];

+endfunction

+sol=ode(y0,0,x,g);

+count=1;

+disp(sol)

+for i=1:2:81

+    a(count)=sol(i);

+    b(count)=sol(i+1);

+    count=count+1;

+end

+plot(x,a)

+plot(x,b,".-")

+clf();

+y0=[%pi/4 0];

+sol=ode(y0,0,x,f);

+count=1;

+disp(sol)

+for i=1:2:81

+    a(count)=sol(i);

+    b(count)=sol(i+1);

+    count=count+1;

+end

+plot(x,a)

+plot(x,b,".-")

+

+xtitle("y vs x","x","y")

+sol=ode(y0,0,x,g);

+count=1;

+disp(sol)

+for i=1:2:81

+    a(count)=sol(i);

+    b(count)=sol(i+1);

+    count=count+1;

+end

+plot(x,a,'o-')

+plot(x,b,"x-")

+h1=legend(["y1","y2","y3","y4"])

diff --git a/964/CH25/EX25.11/25_11_1.jpeg b/964/CH25/EX25.11/25_11_1.jpeg
new file mode 100755
index 000000000..b8f1e82d5
--- /dev/null
+++ b/964/CH25/EX25.11/25_11_1.jpeg
diff --git a/964/CH25/EX25.11/25_11_2.jpeg b/964/CH25/EX25.11/25_11_2.jpeg
new file mode 100755
index 000000000..952c87781
--- /dev/null
+++ b/964/CH25/EX25.11/25_11_2.jpeg
diff --git a/964/CH25/EX25.12/25_12.sce b/964/CH25/EX25.12/25_12.sce
new file mode 100755
index 000000000..ae67f03a0
--- /dev/null
+++ b/964/CH25/EX25.12/25_12.sce
@@ -0,0 +1,39 @@
+//clc()

+disp("f(x,y) = 4*exp(0.8*x) - 0.5*y")

+//f'(x,y) = 4*exp(0.8*x) - 0.5*y

+x = 0:2:2;

+y(1) = 2;

+h =2;

+t = 14.84392;

+k1 = 4*exp(0.8*x(1)) - 0.5*y(1);

+x1 = x(1) + h/2;

+y1 = y(1) + k1*h/2;

+k2 = 4*exp(0.8*x1) - 0.5*y1;

+x1 = x(1) + h/2;

+y1 = y(1) + k2*h/2;

+k3 = 4*exp(0.8*x1) - 0.5*y1;

+x1 = x(1) + h;

+y1 = y(1) + k3*h;

+k4 = 4*exp(0.8*x1) - 0.5*y1;

+y(2) = y(1) + (k1 + 2*k2 + 2*k3 + k4)*h/6;

+e = (t - y(2))/(t);

+disp(y(1:2),"y by h = 2 is")

+disp(e,"error = ")

+h = 1;

+x = 0:h:2;

+for i=1:3

+    k1(i) = 4*exp(0.8*x(i)) - 0.5*y(i);

+    x1 = x(i) + h/2;

+    y1 = y(i) + k1(i)*h/2;

+    k2(i) = 4*exp(0.8*x1) - 0.5*y1;

+    x1 = x(i) + h/2;

+    y1 = y(i) + k2(i)*h/2;

+    k3(i) = 4*exp(0.8*x1) - 0.5*y1;

+    x1 = x(i) + h;

+    y1 = y(i) + k3(i)*h;

+    k4(i) = 4*exp(0.8*x1) - 0.5*y1;

+    y(i+1) = y(i) + (k1(i) + 2*k2(i) + 2*k3(i) + k4(i))*h/6;    

+end

+e = (t - (y(3)))/t;

+disp(y(1:3),"y by h = 1 is")

+disp(e,"error = ")

diff --git a/964/CH25/EX25.13/25_13.sce b/964/CH25/EX25.13/25_13.sce
new file mode 100755
index 000000000..a0c5ba8ae
--- /dev/null
+++ b/964/CH25/EX25.13/25_13.sce
@@ -0,0 +1,27 @@
+//clc()

+disp("f(x,y) = 4*exp(0.8*x) - 0.5*y")

+//f'(x,y) = 4*exp(0.8*x) - 0.5*y

+h =2;

+x = 0:h:2;

+y(1) = 2;

+t = 14.84392;

+k1 = 4*exp(0.8*x(1)) - 0.5*y(1);

+x1 = x(1) + h/5;

+y1 = y(1) + k1*h/5;

+k2 = 4*exp(0.8*x1) - 0.5*y1;

+x1 = x(1) + 3*h/10;

+y1 = y(1) + 3*k1*h/40 + 9*k2*h/40;

+k3 = 4*exp(0.8*x1) - 0.5*y1;

+x1 = x(1) + 3*h/5;

+y1 = y(1) + 3*k1*h/10 - 9*k2*h/10 + 6*k3*h/5;

+k4 = 4*exp(0.8*x1) - 0.5*y1;

+x1 = x(1) + h;

+y1 = y(1) -11*k1*h/54 + 5*k2*h/2 - 70*k3*h/27 + 35*k4*h/27;

+k5 = 4*exp(0.8*x1) - 0.5*y1;

+x1 = x(1) + 7*h/8;

+y1 = y(1) + 1631*k1*h/55296 + 175*k2*h/512 + 575*k3*h/13824 + 44275*k4*h/110592 +253*k5*h/4096;

+k6 = 4*exp(0.8*x1) - 0.5*y1;

+y1 = y(1) + (37*k1/378 + 250*k3/621 + 125*k4/594 + 512*k6/1771)*h;

+y2 = y(1) + (2825*k1/27648 + 18575*k3/48384 + 13525*k4/55296 + 277*k5/14336 + k6/4)*h;

+disp(y1,"y ( fourth order prediction ) = ")

+disp(y2,"y ( fifth order prediction ) = ")

diff --git a/964/CH25/EX25.14/25_14.jpeg b/964/CH25/EX25.14/25_14.jpeg
new file mode 100755
index 000000000..73b0a44b9
--- /dev/null
+++ b/964/CH25/EX25.14/25_14.jpeg
diff --git a/964/CH25/EX25.14/25_14.sce b/964/CH25/EX25.14/25_14.sce
new file mode 100755
index 000000000..66b236d37
--- /dev/null
+++ b/964/CH25/EX25.14/25_14.sce
@@ -0,0 +1,10 @@
+clc;

+clear;

+function yp=f(x,y)

+    yp=10*exp(-(x-2)^2/(2*(0.075^2)))-0.6*y

+endfunction

+x=0:0.1:4

+y0=0.5;

+sol=ode(y0,0,x,f);

+plot(x,sol)

+xtitle("y vs x","x","y")

diff --git a/964/CH25/EX25.2/25_2.sce b/964/CH25/EX25.2/25_2.sce
new file mode 100755
index 000000000..02feb97b4
--- /dev/null
+++ b/964/CH25/EX25.2/25_2.sce
@@ -0,0 +1,11 @@
+//clc()

+//f(x,y) = dy/dx = -2*x^3 + 12*x^2 - 20*x + 8.5

+//f'(x,y) = -6*x^2 + 24*x - 20

+//f"(x,y) = -12*x + 24

+//f"'(x,y) = -12

+x = 0;

+Et2 = (-6*x^2 + 24*x - 20) * 0.5^2 / 2;

+Et3 = (-12*x + 24) * (0.5)^3 / 6;

+Et4 = (-12) *(0.5 ^ 4) / 24;

+Et = Et2 + Et3 + Et4;

+disp(Et,"Total truncation error =")

diff --git a/964/CH25/EX25.3/25_3.sce b/964/CH25/EX25.3/25_3.sce
new file mode 100755
index 000000000..7cf519d87
--- /dev/null
+++ b/964/CH25/EX25.3/25_3.sce
@@ -0,0 +1,20 @@
+//clc()

+//dy/dx = -2*x^3 + 12*x^2 - 20*x + 8.5

+//therefore, y = -0.5*x^4 + 4*x^3 - 10*x^2 + 8.5 + c

+x1 = 0;

+y1 = 1;

+h = 0.25;

+c =-(-0.5*x1^4 + 4*x1^3 - 10*x1^2 + 8.5*x1 - y1);

+x = 0:h:4;

+disp(x,"x = ")

+y = -0.5*x^4 + 4*x^3 - 10*x^2 + 8.5*x + c;

+disp(y,"true values of y = ")

+fxy = -2*x^3 + 12*x^2 - 20*x + 8.5;

+y2(1) = y(1);

+e(1) = (y(1) - y2(1)) * 100 / y(1);

+for i = 2:17

+    y2(i) = y2(i-1) + fxy(i-1)*h;

+    e(i) = (y(i) - y2(i))*100/y(i);

+end

+disp(y2,"y by euler method =")

+disp(e,"error =")

diff --git a/964/CH25/EX25.4/25_4.jpeg b/964/CH25/EX25.4/25_4.jpeg
new file mode 100755
index 000000000..88a957dfb
--- /dev/null
+++ b/964/CH25/EX25.4/25_4.jpeg
diff --git a/964/CH25/EX25.4/25_4.sce b/964/CH25/EX25.4/25_4.sce
new file mode 100755
index 000000000..e2ef29404
--- /dev/null
+++ b/964/CH25/EX25.4/25_4.sce
@@ -0,0 +1,24 @@
+clc;

+clear;

+m=68.1;

+g=9.8;

+c=12.5;

+a=8.3;

+b=2.2;

+vmax=46;

+function yp=f(t,v)

+    yp=g-c*v/m;

+endfunction

+v0=0;

+t=0:15;

+x=ode(v0,0,t,f);

+disp(x)

+plot(t,x,'.-')

+

+function yp=f1(t,v)

+    yp=g-(c/m)*(v+a*(v/vmax)^b)

+endfunction

+x1=ode(v0,0,t,f1);

+plot(t,x1)

+xtitle("velocity vs time","t (s)","v (m/s)")

+h1=legend(["Linear";"Nonlinear"])

diff --git a/964/CH25/EX25.5/25_5.sce b/964/CH25/EX25.5/25_5.sce
new file mode 100755
index 000000000..111e8d0ca
--- /dev/null
+++ b/964/CH25/EX25.5/25_5.sce
@@ -0,0 +1,21 @@
+//clc()

+//y' = 4*exp(0.8*x) - 0.5*y

+//y = 4*(exp(0.8*x) - exp(-0.5*x))/1.3 + 2*exp(-0.5*x)

+x = 0:1:4;

+disp(x)

+x1 = 0;

+y1 = 2;

+y2(1) = y1;

+for i = 1:5

+    y(i) = 4*(exp(0.8*x(i)) - exp(-0.5*x(i)))/1.3 + 2*exp(-0.5*x(i));

+    dy(i) = 4*exp(0.8*x(i)) - 0.5*y2(i);

+    y2(i + 1) = y2(i) + dy(i);

+    if i>1 then

+        m(i) = (dy(i) + dy(i-1))/2;

+        y2(i) = y2(i-1) + m(i);

+        dy(i) = 4*exp(0.8*x(i)) - 0.5*y2(i);

+    end

+    e(i) = (y(i) - y2(i)) * 100 / y(i);

+end

+disp(y2(1:5),"By heuns method(1 iteration)")

+disp("%",e(1:5),"error = ")

diff --git a/964/CH25/EX25.6/25_6.sce b/964/CH25/EX25.6/25_6.sce
new file mode 100755
index 000000000..9bc14bc87
--- /dev/null
+++ b/964/CH25/EX25.6/25_6.sce
@@ -0,0 +1,28 @@
+//clc()

+//f'(x,y) = -2*x^3 + 12*x^2 -20*x + 8.5

+//f(x,y) = -x^4 / 2 + 4*x^3 - 10*x^2 + 8.5*x + 1

+h = 0.5;

+x = 0:h:4;

+y1 = -x^4 / 2 + 4*x^3 - 10*x^2 + 8.5*x + 1;

+y(1) = 1;

+disp(x,"x =")

+disp(y1,"true value of y =")

+for i = 1:8

+    k1(i) = -2*x(i)^3 + 12*x(i)^2 -20*x(i) + 8.5;

+    x1(i) = x(i) + h/2;

+    k2(i) = -2*x1(i)^3 + 12*x1(i)^2 -20*x1(i) + 8.5;

+    y(i+1) = y(i) + k2(i)*h;

+    e(i) = (y1(i) - y(i))*100/y1(i);

+end

+disp(y(1:9),"y by midpoint method")

+disp(e,"error = ")

+for i = 1:8

+    k1(i) = -2*x(i)^3 + 12*x(i)^2 -20*x(i) + 8.5;

+    x(i) = x(i) + 3*h/4;

+    k2(i) = -2*x(i)^3 + 12*x(i)^2 -20*x(i) + 8.5;

+    y(i+1) = y(i) + (k1(i)/3 + 2*k2(i)/3)*h;

+    e(i) = (y1(i) - y(i))*100/y1(i);

+end

+disp(y(1:9),"y by second order Ralston RK")

+disp(e,"error = ")

+

diff --git a/964/CH25/EX25.7/25_7.sce b/964/CH25/EX25.7/25_7.sce
new file mode 100755
index 000000000..ce01ffb75
--- /dev/null
+++ b/964/CH25/EX25.7/25_7.sce
@@ -0,0 +1,33 @@
+//clc()

+//f'(x,y) = -2*x^3 + 12*x^2 -20*x + 8.5

+//f(x,y) = -x^4 / 2 + 4*x^3 - 10*x^2 + 8.5*x + 1

+h = 0.5;

+x = 0:h:4;

+y1 = -x^4 / 2 + 4*x^3 - 10*x^2 + 8.5*x + 1;

+y(1) = 1;

+for i=1:8

+     k1(i) = -2*x(i)^3 + 12*x(i)^2 -20*x(i) + 8.5;

+    x1(i) = x(i) + h/2;

+    k2(i) = -2*x1(i)^3 + 12*x1(i)^2 -20*x1(i) + 8.5;

+    k3(i) = -2*x1(i)^3 + 12*x1(i)^2 -20*x1(i) + 8.5;

+    x2(i) = x(i) + h;

+    k4(i) = -2*x2(i)^3 + 12*x2(i)^2 -20*x2(i) + 8.5;

+    y(i+1) = y(i) + (k1(i)+2*k2(i)+2*k3(i)+k4(i))*h/6;

+    e(i) = (y1(i) - y(i))*100/y1(i);

+end

+disp("f(x,y) = -2*x^3 + 12*x^2 -20*x + 8.5")

+disp(y(1:9),"y by fourth order Ralston RK")

+disp("f(x,y) = 4*exp(0.8*x) - 0.5*y")

+x = 0:h:0.5;

+y(1) = 2;

+k1 = 4*(exp(0.8*x(1)))-0.5*y(1);

+x1 = x(1) + 0.5*h;

+y1 = y(1) + 0.5*k1*h;

+k2 = 4*exp(0.8*x1) - 0.5*y1;

+y2 = y(1) + 0.5*k2*h;

+k3 = 4*exp(0.8*x1) - 0.5*y2;

+x1 = x(1) + h;

+y1 = y(1) + k3*h;

+k4 = 4*exp(0.8*x1) - 0.5*y1;

+y(2) = y(1) + (k1+2*k2+2*k3+k4)*h/6;

+disp(y(1:2),"y = ")

diff --git a/964/CH25/EX25.8/25_8.sce b/964/CH25/EX25.8/25_8.sce
new file mode 100755
index 000000000..9ec20d157
--- /dev/null
+++ b/964/CH25/EX25.8/25_8.sce
@@ -0,0 +1,57 @@
+//clc()

+disp("f(x,y) = 4*exp(0.8*x) - 0.5*y")

+h = 1;

+x = 0:h:4;

+y(1) = 2;

+for i = 1:5

+    k1(i) = 4*(exp(0.8*x(i)))-0.5*y(i);

+    x1 = x(i) + h;

+    y1 = y(i) + k1(i)*h;

+    k2(i) = 4*(exp(0.8*x1))-0.5*y1;

+    y(i+1) = y(i) + (k1(i)/2 + k2(i)/2)*h;

+end

+disp(y(1:5),"y(second order RK method) = ")

+for i = 1:5

+    k1(i) = 4*(exp(0.8*x(i)))-0.5*y(i);

+    x1 = x(i) + 0.5*h;

+    y1 = y(i) + 0.5*h*k1(i);

+    k2(i) = 4*(exp(0.8*x1))-0.5*y1;

+    x1 = x(i) + h;

+    y1 = y(i) -k1(i)*h + 2*k2(i)*h;

+    k3(i) = 4*(exp(0.8*x1))-0.5*y1;

+    y(i+1) = y(i) + (k1(i) + 4*k2(i) + k3(i))*h/6;

+end

+disp(y(1:5),"y(third order RK method) = ")

+for i = 1:5

+k1(i) = 4*(exp(0.8*x(i)))-0.5*y(i);

+x1 = x(i) + 0.5*h;

+y1 = y(i) + 0.5*k1(i)*h;

+k2(i) = 4*exp(0.8*x1) - 0.5*y1;

+y2 = y(i) + 0.5*k2(i)*h;

+k3(i) = 4*exp(0.8*x1) - 0.5*y2;

+x1 = x(i) + h;

+y1 = y(i) + k3(i)*h;

+k4(i) = 4*exp(0.8*x1) - 0.5*y1;

+y(i+1) = y(i) + (k1(i)+2*k2(i)+2*k3(i)+k4(i))*h/6;

+end

+disp(y(1:5),"y(fourth order RK method) = ")

+for i = 1:5

+    

+k1(i) = 4*(exp(0.8*x(i)))-0.5*y(i);

+x1 = x(i) + 0.25*h;

+y1 = y(i) + 0.25*k1(i)*h;

+k2(i) = 4*exp(0.8*x1) - 0.5*y1;

+y2 = y(i) + 0.125*k2(i)*h + 0.125*k1(i)*h;

+k3(i) = 4*exp(0.8*x1) - 0.5*y2;

+x1 = x(i) + 0.5*h;

+y1 = y(i) -0.5*k2(i)*h + k3(i)*h;

+k4(i) = 4*exp(0.8*x1) - 0.5*y1;

+x1 = x(i) + 0.75*h;

+y1 = y(i) + 3*k1(i)*h/16 + 9*k4(i)*h/16;

+k5(i) = 4*exp(0.8*x1) - 0.5*y1;

+x1 = x(i) + h;

+y1 = y(i) - 3*k1(i)*h/7 + 2*k2(i)*h/7 + 12*k3(i)*h/7 - 12*k4(i)*h/7 + 8*k5(i)*h/7;

+k6(i) = 4*exp(0.8*x1) - 0.5*y1;

+y(i+1) = y(i) + (7*k1(i)+32*k3(i)+12*k4(i) + 32*k5(i) + 7*k6(i))*h/90;

+end

+disp(y(1:5),"y(fifth order RK method)")

diff --git a/964/CH25/EX25.9/25_9.sce b/964/CH25/EX25.9/25_9.sce
new file mode 100755
index 000000000..ba1c4f88b
--- /dev/null
+++ b/964/CH25/EX25.9/25_9.sce
@@ -0,0 +1,15 @@
+//clc()

+//dy1/dx = -0.5*y1

+//dy2/dx = 4 - 0.3*y2 - 0.1*y1

+x1 = 0;

+h =0.5;

+y1(1) = 4;

+y2(1) = 6;

+x = 0:h:2;

+for i = 2:5

+    y1(i) = y1(i-1)  -0.5*y1(i-1)*h;

+    y2(i) = y2(i-1) + (4 - 0.3*y2(i-1) - 0.1*y1(i-1))*h;

+end

+disp(x,"x = ")

+disp(y1,"y1 = ")

+disp(y2,"y2 = ")