Part F

README.md

@@ -104,3 +104,53 @@ end
 se = E*t*h^2/(2*(1-v^2))*sum(sum(0.25.*dwdx.^4+.25.*dwdy.^4+0.5.*(dwdx.*dwdy).^2));
 pw_se = pw-se;
 ```
+
+## Part F
+```
+function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)
+if nargin<3,error('at least 3 input arguments required'),end
+test = func(xl,varargin{:})*func(xu,varargin{:});
+if nargin<4||isempty(es), es=0.0001;end
+if nargin<5||isempty(maxit), maxit=50;end
+iter = 0; xr = xl; ea = 100;
+while (1)
+  xr_old = xr;
+  xr = (xl + xu)/2;
+  iter = iter + 1;
+  if xr ~= 0,ea = abs((xr - xr_old)/xr) * 100;end
+  test = func(xl,varargin{:})*func(xr,varargin{:});
+  if test < 0
+    xu = xr;
+  elseif test > 0
+    xl = xr;
+  else
+    ea = 0;
+  end
+  if ea <= es || iter >= maxit,break,end
+end
+root = xr; fx = func(xr, varargin{:});
+```
+```
+function [T,ea] = tension_sol(P,n)
+y =@(T) SE_diff(T,P,n);
+[T,fx,ea,iter]=bisect(y,.01,1);
+```
+```
+n=[3,20:5:40];
+P=0.001;
+T = zeros(1,length(n));
+ea = zeros(1,length(n));
+
+ for i = 1:length(n)
+[T(i), ea(i)] = tension_sol(P,n(i));
+ end
+ ```
+
+ |number of nodes|Tension (uN/um)| rel. error|
+|---|---|---|
+|3 |0.049 |n/a|
+|20|0.0599|20.5%|
+|25|0.0601|0.27%|
+|30|0.0602|0.16%|
+|35|0.0602|0.09%|
+|40|0.0603|0.05%|

SE_diff.asv

@@ -1,27 +1,58 @@
-function [pw_se,w]=SE_diff(T,P,n)
-
-E = 1e6;
-v = .31;
-t = .0003;
-h = 10/(n+1);
-w = membrane_solution(T,P,n);
-z = zeros(n+2);
-z(2:end-1,2:end-1) = reshape(w,[n n]);
-nt = n + 1;
-wn = zeros(nt);
-for i = 1:nt
-    for j = 1:nt
-        wn(i,j) = mean([z(i,j),z(i+1,j),z(i,j+1),z(i+1,j+1)]);
+function [pw_se,w] = SE_diff(T,P,n)
+
+    % Call function from part c to obtain vector w
+    w = membrane_solution(T,P,n);
+    m = zeros(n+2);
+    c = vec2mat(w,n);
+    m(2:n+1,2:n+1) = c;
+
+
+    h = 10/(n+1);
+
+    w_bar = zeros(n+1,n+1);
+    for i = 1:n+1
+        for j = 1:n+1
+            w_bar(i,j) = (m(i,j)+m(i,j+1)+m(i+1,j)+m(i+1,j+1))/4;
+        end
+    end
+
+    left_side = zeros(n+1,n+1);
+    for  i = 1:(n+1)^2
+        left_side(i) = P * h^2*w_bar(i);
     end
-end
-pw = sum(sum(wn.*h^2.*P));
-dwdx = zeros(nt);
-dwdy = zeros(nt);
-for i = 1:nt
-    for j = 1:nt
-        dwdx(i,j) = mean([z(i+1,j)-z(i,j),z(i+1,j+1)-z(i,j+1)]);
-        dwdy(i,j) = mean([z(i,j+1)-z(i,j),z(i+1,j+1)-z(i+1,j)]);
+    L = sum(sum(left_side));
+
+    dx = zeros(n+2,n+1);
+    for i = 2:n+2
+        dx(:,i-1)= (m(:,i)-m(:,i-1))/(h);
     end
-end
-se = E*t*h^2/(2*(1-v^2))*sum(sum(0.25.*dwdx.^4+.25.*dwdy.^4+0.5.*(dwdx.*dwdy).^2));
-pw_se = pw-se;
+
+    dw_dx = zeros(n+1,n+1);
+    for i = 1:n+1
+        for j = 1:n+1
+        dw_dx(i,j) = (dx(i,j)+dx(i+1,j))/2;
+        end 
+    end 
+
+    dy = zeros(n+1,n+2);
+    for i = 2:n+2
+        dy(i-1,:)= (m(i-1,:)-m(i,:))/(h);
+    end
+
+    dw_dy = zeros(n+1,n+1);
+    for i = 1:n+1
+        for j = 1:n+1
+        dw_dy(i,j) = (dy(i,j)+dy(i,j+1))/2;
+        end 
+    end 
+
+    right_side = zeros(n+1,n+1);
+
+    for i = 1:(n+1)^2
+        right_side(i) = (0.25*(dw_dx(i))^4+0.25*(dw_dy(i))^4+0.5*(dw_dx(i)*dw_dy(i))^2);
+    end 
+    R = sum(sum(right_side));
+    R = (1000*10^3*0.3*10^-3*h^2)/(2*(1-0.31^2))*R;
+
+    pw_se = R - L;
+end

bisect.asv

@@ -0,0 +1,23 @@
+function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)
+if nargin<3,error('at least 3 input arguments required'),end
+test = func(xl,varargin{:})*func(xu,varargin{:});
+if test>0,error('no sign change'),end
+if nargin<4||isempty(es), es=0.0001;end
+if nargin<5||isempty(maxit), maxit=50;end
+iter = 0; xr = xl; ea = 100;
+while (1)
+  xrold = xr;
+  xr = (xl + xu)/2;
+  iter = iter + 1;
+  if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end
+  test = func(xl,varargin{:})*func(xr,varargin{:});
+  if test < 0
+    xu = xr;
+  elseif test > 0
+    xl = xr;
+  else
+    ea = 0;
+  end
+  if ea <= es || iter >= maxit,break,end
+end
+root = xr; fx = func(xr, varargin{:});

bisect.m

@@ -0,0 +1,22 @@
+function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)
+if nargin<3,error('at least 3 input arguments required'),end
+test = func(xl,varargin{:})*func(xu,varargin{:});
+if nargin<4||isempty(es), es=0.0001;end
+if nargin<5||isempty(maxit), maxit=50;end
+iter = 0; xr = xl; ea = 100;
+while (1)
+  xr_old = xr;
+  xr = (xl + xu)/2;
+  iter = iter + 1;
+  if xr ~= 0,ea = abs((xr - xr_old)/xr) * 100;end
+  test = func(xl,varargin{:})*func(xr,varargin{:});
+  if test < 0
+    xu = xr;
+  elseif test > 0
+    xl = xr;
+  else
+    ea = 0;
+  end
+  if ea <= es || iter >= maxit,break,end
+end
+root = xr; fx = func(xr, varargin{:});

partFsetup.m

@@ -0,0 +1,8 @@
+n=[3,20:5:40];
+P=0.001;
+T = zeros(1,length(n));
+ea = zeros(1,length(n));
+
+ for i = 1:length(n)
+[T(i), ea(i)] = tension_sol(P,n(i));
+ end

rel_error.m

@@ -0,0 +1,5 @@
+function re = rel_error (T)
+re = zeros(1,length(T)-1);
+for i = 2:length(T)
+    re(i-1)= abs(T(i)-T(i-1))/T(i-1);
+end

tension_sol.m

@@ -0,0 +1,3 @@
+function [T,ea] = tension_sol(P,n)
+y =@(T) SE_diff(T,P,n);
+[T,fx,ea,iter]=bisect(y,.01,1);