summaryrefslogtreecommitdiff
path: root/845/CH9
diff options
context:
space:
mode:
Diffstat (limited to '845/CH9')
-rwxr-xr-x845/CH9/EX9.1/Ex9_1.sce2
-rwxr-xr-x845/CH9/EX9.2/Ex9_2.sce47
-rwxr-xr-x845/CH9/EX9.3/Ex9_3.sce42
-rwxr-xr-x845/CH9/EX9.4/Ex9_4.sce46
-rwxr-xr-x845/CH9/EX9.5/Ex9_5.sce45
-rwxr-xr-x845/CH9/EX9.6/Ex9_6.sce2
-rwxr-xr-x845/CH9/EX9.7/Ex9_7.sce88
7 files changed, 272 insertions, 0 deletions
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:")
generated by cgit (git 2.34.1) at 2025-03-13 01:51:58 +0000