diff --git a/README.md b/README.md index 1ea9650..243bbf4 100644 --- a/README.md +++ b/README.md @@ -68,3 +68,90 @@ Solve for the vector w using the function from part C with a P value of 0.001, T ``` ![nxn Interior Nodes, n = 10](./figures/figure2.png) + +### Part E +``` matlab +function [pw_se,w] = SE_diff(T,P,n) + + % Call function from part c to obtain vector w + w = membrane_solution(T,P,n); + % Create All Nodes Matrix (Interior & Exterior) + m = zeros(n+2); + c = vec2mat(w,n); % Vector w rearranged to n x n matrix. + m(2:n+1,2:n+1) = c; % insert matrix c into middle of m matrix + % m = symmetric matrix + + h = 10/(n+1); % h = delta_x = delta_y (micron) + + % Compute w bar + 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 of the equation + left_side = zeros(n+1,n+1); + for i = 1:(n+1)^2 + left_side(i) = P * h^2*w_bar(i); + end + L = sum(sum(left_side)); + + % Compute dw_dx + dx = zeros(n+2,n+1); + for i = 2:n+2 + dx(:,i-1)= (m(:,i)-m(:,i-1))/(h); + end + + dw_dx = zeros(n+1,n+1); + for i = 1:n+1 % rows + for j = 1:n+1 % columns + dw_dx(i,j) = (dx(i,j)+dx(i+1,j))/2; + end + end + + % Compute dw_dy + 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 % rows + for j = 1:n+1 % columns + 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 +``` + +- At 10 nodes, pw_se = 53.4902 + +### Part F + +``` matlab +[pw_se,w] = SE_diff(T,P,n); +[root,fx,ea,iter] = bisect_final_project(@(T) SE_diff(T,0.001,25),.001,1,.1) +``` + +The same code was run five times, changing the value of n from 20 to 40 for intervals of 5. + +| number of nodes | Tension (uN/um) | rel. error| +| --- | --- | --- | +| 3 | 0.0489 | n/a | +| 20 | 0.0599 | 18.4% | +| 25 | 0.0601 | 0.30% | +| 30 | 0.0602 | 0.17% | +| 35 | 0.0603 | 0.16% | +| 40 | 0.0603 | 0.00% | diff --git a/bisect_final_project.m~ b/bisect_final_project.m~ new file mode 100644 index 0000000..0fae3c9 --- /dev/null +++ b/bisect_final_project.m~ @@ -0,0 +1,37 @@ +function [root,fx,ea,iter]=bisect_final_project(func,xl,xu,es,maxit,varargin) +% bisect: root location zeroes +% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...): +% uses bisection method to find the root of func +% input: +% func = name of function +% xl, xu = lower and upper guesses +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by func +% output: +% root = real root +% fx = function value at root +% ea = approximate relative error (%) +% iter = number of iterations +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{:}); \ No newline at end of file diff --git a/bisect_final_project_script.m b/bisect_final_project_script.m new file mode 100644 index 0000000..a3e3cf0 --- /dev/null +++ b/bisect_final_project_script.m @@ -0,0 +1,2 @@ +[pw_se,w] = SE_diff(T,P,n); +[root,fx,ea,iter] = bisect_final_project(@(T) SE_diff(T,0.001,40),.001,1,.1) \ No newline at end of file diff --git a/partg.m b/partg.m new file mode 100644 index 0000000..1ef04b4 --- /dev/null +++ b/partg.m @@ -0,0 +1,16 @@ + +P = linspace(0.001,0.01,10); + +for i = 1:length(P) + func = @(T) SE_diff(T,P(i),10); + [root(i),fx,ea,iter] = bisect_final_project(func,0.001,1,.1); +end + +for i = 1:length(root) + w = membrane_solution(root(i),P(i),10); + w1(:,i) = w; % displays w values as column +end + +w_max = max(w1); +v_cub = interp1(w_max,P,'pchip'); +plot(w_max,P) \ No newline at end of file