From b1f5c3f8d6671b4331cef1dcebdf63b7a43a3a2b Mon Sep 17 00:00:00 2001 From: priyanka Date: Wed, 24 Jun 2015 15:03:17 +0530 Subject: initial commit / add all books --- 845/CH9/EX9.1/Ex9_1.sce | 2 ++ 845/CH9/EX9.2/Ex9_2.sce | 47 ++++++++++++++++++++++++++ 845/CH9/EX9.3/Ex9_3.sce | 42 +++++++++++++++++++++++ 845/CH9/EX9.4/Ex9_4.sce | 46 ++++++++++++++++++++++++++ 845/CH9/EX9.5/Ex9_5.sce | 45 +++++++++++++++++++++++++ 845/CH9/EX9.6/Ex9_6.sce | 2 ++ 845/CH9/EX9.7/Ex9_7.sce | 88 +++++++++++++++++++++++++++++++++++++++++++++++++ 7 files changed, 272 insertions(+) create mode 100755 845/CH9/EX9.1/Ex9_1.sce create mode 100755 845/CH9/EX9.2/Ex9_2.sce create mode 100755 845/CH9/EX9.3/Ex9_3.sce create mode 100755 845/CH9/EX9.4/Ex9_4.sce create mode 100755 845/CH9/EX9.5/Ex9_5.sce create mode 100755 845/CH9/EX9.6/Ex9_6.sce create mode 100755 845/CH9/EX9.7/Ex9_7.sce (limited to '845/CH9') diff --git a/845/CH9/EX9.1/Ex9_1.sce b/845/CH9/EX9.1/Ex9_1.sce new file mode 100755 index 000000000..939a282c2 --- /dev/null +++ b/845/CH9/EX9.1/Ex9_1.sce @@ -0,0 +1,2 @@ +// Example 9.1 +// This is an analytical problem and need not be coded. diff --git a/845/CH9/EX9.2/Ex9_2.sce b/845/CH9/EX9.2/Ex9_2.sce new file mode 100755 index 000000000..642893b92 --- /dev/null +++ b/845/CH9/EX9.2/Ex9_2.sce @@ -0,0 +1,47 @@ +//Example 9.2 + +clc +clear + +delx = 0.1; +delt = 0.002; +xf = 1; +tf = 0.006; +x = 0:delx:xf; +t = 0:delt:tf; +m = length(x); +n = length(t); +lamda = delt/delx^2; + +y = zeros(n,m); +y(1:n,1) = 0; +y(1:n,m) = 0; +y(1,1:m) = sin(%pi*x); +for k = 2:n + M1 = zeros(m-2); + M2 = zeros(m-2,1); + for i = 1:m-2 + M1(i,i) = 1+2*lamda; + if i==1 + M1(i,i+1) = -lamda; + M2(i) = y(k-1,i+1) + lamda*y(k,i); + elseif i==m-2 + M1(i,i-1) = -lamda; + M2(i) = y(k-1,i+1) + lamda*y(k,i+2); + else + M1(i,i+1) = -lamda; + M1(i,i-1) = -lamda; + M2(i) = y(k-1,i+1); + end + end + y(k,2:m-1) = (M1\M2)'; +end +y = round(y*10^4)/10^4; +mprintf("%4s %7s %9s %8s %9s %9s %9s %9s %9s %9s %9s %9s %9s",'n','t','x = 0.0','x = 0.1','x = 0.2','x = 0.3','x = 0.4','x = 0.5','x = 0.6','x = 0.7','x = 0.8','x = 0.9','x = 1.0'); +disp([(0:n-1)' t' y]) + +disp("At t = 0.006:") +disp(y(n,1:m),"Computed Solution:") +Texact = exp(-%pi^2*tf)*sin(%pi*x); +Texact = round(Texact*10^4)/10^4; +disp(Texact,"Analytical Solution:") diff --git a/845/CH9/EX9.3/Ex9_3.sce b/845/CH9/EX9.3/Ex9_3.sce new file mode 100755 index 000000000..ae53f775f --- /dev/null +++ b/845/CH9/EX9.3/Ex9_3.sce @@ -0,0 +1,42 @@ +//Example 9.3 + +clc +clear + +delx = 0.1; +delt = 0.001; +xf = 0.5; +tf = 0.003; +x = 0:delx:xf; +t = 0:delt:tf; +m = length(x); +n = length(t); +r = delt/delx^2; + +T = zeros(m,n); +T(1:m,1) = 0; +delTxi = 0; +delTxf = 1; + +for j = 1:n + M1 = zeros(m,m); + M2 = zeros(m,1); + for i = 1:m + if i == 1 then + M1(i,i) = 1; + M1(i,i+1) = -1; + M2(i) = -delx * delTxi; + elseif i == m then + M1(i,i) = 1; + M1(i,i-1) = -1; + M2(i) = delx * delTxf; + else + M1(i,i) = 1; + M2(i) = r*T(i+1,j) + (1-2*r) * T(i,j) + r*T(i-1,j); + end + end + T(1:m,j+1) = (M1\M2); +end +T = T(:,2:n+1); +mprintf("%4s %7s %9s %5s %7s %9s %9s %9s",'n','t','x = 0.0','x=0.1','x = 0.2','x = 0.3','x = 0.4','x = 0.5'); +disp([(0:n-1)' t' T']) diff --git a/845/CH9/EX9.4/Ex9_4.sce b/845/CH9/EX9.4/Ex9_4.sce new file mode 100755 index 000000000..a50cc5902 --- /dev/null +++ b/845/CH9/EX9.4/Ex9_4.sce @@ -0,0 +1,46 @@ +//Example 9.4 + +clc +clear + +delx = 0.25; +delt = 1/32; +xf = 1; +tf = delt; +x = 0:delx:xf; +t = 0:delt:tf; +m = length(x); +n = length(t); +r = delt/delx^2; + + +T = zeros(m,n); +T(1:m,1) = 1; +T(1,1:n) = 0; +T(m,1:n) = 0; + +for j = 1:n-1 + M1 = zeros(m-2,m-2); + M2 = zeros(m-2,1); + for i = 2:m-1 + if i == 2 then + M1(i-1,i-1) = -2*(1+r); + M1(i-1,i) = r; + M2(i-1) = -(r*T(i+1,j) + 2*(1-r)*T(i,j) + r*T(i-1,j)); + elseif i == m-1 + M1(i-1,i-2) = r; + M1(i-1,i-1) = -2*(1+r); + M2(i-1) = -(r*T(i+1,j) + 2*(1-r)*T(i,j) + r*T(i-1,j)); + else + M1(i-1,i-2) = r; + M1(i-1,i-1) = -2*(1+r); + M1(i-1,i) = r; + M2(i-1) = -(r*T(i+1,j) + 2*(1-r)*T(i,j) + r*T(i-1,j)); + end + end + T(2:m-1,j+1) = (M1\M2); +end +T1 = M1\M2; +for i = 1:length(T1) + disp(strcat(["T",string(i)," = ",string(T1(i))])); +end diff --git a/845/CH9/EX9.5/Ex9_5.sce b/845/CH9/EX9.5/Ex9_5.sce new file mode 100755 index 000000000..19884c709 --- /dev/null +++ b/845/CH9/EX9.5/Ex9_5.sce @@ -0,0 +1,45 @@ +//Example 9.5 + +clc +clear + +delx = 1; +delt = 1.893; +alpha = 0.132; +xf = 4; +tf = delt; +x = 0:delx:xf; +t = 0:delt:tf; +m = length(x); +n = length(t); +r = alpha*delt/delx^2; +r = round(r*10^2)/10^2; + +T = zeros(m,n); +T(1:m,1) = 1000; + +for j = 1:n-1 + M1 = zeros(m,m); + M2 = zeros(m,1); + for i = 1:m + if i == 1 then + M1(i,i) = -(2+2.15*r); + M1(i,i+1) = 2*r; + M2(i) = -(2*r*T(2,j) + (2-2.15*r)*T(1,j) + 21*r); + elseif i == m then + M1(i,i) = -(2+1.75*r); + M1(i,i-1) = 2*r; + M2(i) = -(2*r*T(m-1,j) + (2-1.75*r)*T(m,j) - 35*r); + else + M1(i,i-1) = r; + M1(i,i) = -2*(1+r); + M1(i,i+1) = r; + M2(i) = -(r*T(i+1,j) + 2*(1-r)*T(i,j) + r*T(i-1,j)); + end + end + T(1:m,j+1) = (M1\M2); +end +disp(M1,"Coefficient Matrix:") +disp(M2,"Constant Matrix:") +T = round(T*10^4)/10^4; +disp(T',"Table:") diff --git a/845/CH9/EX9.6/Ex9_6.sce b/845/CH9/EX9.6/Ex9_6.sce new file mode 100755 index 000000000..8da94fe79 --- /dev/null +++ b/845/CH9/EX9.6/Ex9_6.sce @@ -0,0 +1,2 @@ +// Example 9.6 +// This is an analytical problem and need not be coded. diff --git a/845/CH9/EX9.7/Ex9_7.sce b/845/CH9/EX9.7/Ex9_7.sce new file mode 100755 index 000000000..e5be14367 --- /dev/null +++ b/845/CH9/EX9.7/Ex9_7.sce @@ -0,0 +1,88 @@ +//Example 9.7 + +clc +clear + +h = 2; +delt = 4; +tf = 8; +xf = 8; +yf = 6; +x = 0:h:xf; +y = 0:h:yf; +t = 0:delt:tf; +m = length(x); +n = length(y); +p = length(t); +r = delt/h^2; +r = round(r*10^2)/10^2; + +T = 50*ones(n,m); +T0 = T; +T(1,1:m) = 110:-10:70; +T(n,1:m) = 0:10:40; +T(2:n-1,1) = [65; 25]; +T(2:n-1,m) = [60; 50]; + +u = (m-2)*(n-2); +index = [repmat(2:m-1,1,n-2); gsort(repmat(2:n-1,1,m-2))]; + +M1 = zeros(u,u); +M2 = zeros(u,1); +for j = 2:m-1 + for i = 2:n-1 + ind = find(index(1,:)== j& index(2,:) == i); + if j == 2 then + M1(ind,ind) = 1+2*r; + M1(ind,ind+1) = -r; + M2(ind) = r*T(i,j-1) + r*T0(i-1,j) + (1-2*r)*T0(i,j) + r*T0(i+1,j); + elseif j == m-1 then + M1(ind,ind-1) = -r; + M1(ind,ind) = 1+2*r; + M2(ind) = r*T(i,j+1) + r*T0(i-1,j) + (1-2*r)*T0(i,j) + r*T0(i+1,j); + else + M1(ind,ind-1) = -r; + M1(ind,ind) = 1+2*r; + M1(ind,ind+1) = -r; + M2(ind) = r*T0(i-1,j) + (1-2*r)*T0(i,j) + r*T0(i+1,j); + end + end +end +value = M1\M2; +value = round(value*10^4)/10^4; +for i = 1:length(index(1,:)) + t = index(:,i); + T(t(2),t(1)) = value(i); +end +disp(T,"At t = 4:") +T0 = T; + +index = gsort(index,'lc','i'); +M1 = zeros(u,u); +M2 = zeros(u,1); +for j = 2:m-1 + for i = 2:n-1 + ind = find(index(1,:)== j& index(2,:) == i); + if i == 2 then + M1(ind,ind) = 1+2*r; + M1(ind,ind+1) = -r; + M2(ind) = r*T(i-1,j) + r*T0(i,j-1) + (1-2*r)*T0(i,j) + r*T0(i,j+1); + elseif i == n-1 then + M1(ind,ind-1) = -r; + M1(ind,ind) = 1+2*r; + M2(ind) = r*T(i+1,j) + r*T0(i,j-1) + (1-2*r)*T0(i,j) + r*T0(i,j+1); + else + M1(ind,ind-1) = -r; + M1(ind,ind) = 1+2*r; + M1(ind,ind+1) = -r; + M2(ind) = + r*T0(i,j-1) + (1-2*r)*T0(i,j) + r*T0(i,j+1); + end + end +end +value = M1\M2; +value = round(value*10^4)/10^4; +for i = 1:length(index(1,:)) + t = index(:,i); + T(t(2),t(1)) = value(i); +end +disp(T,"At t = 8:") -- cgit