
<!DOCTYPE html
  PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head>
      <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
   <!--
This HTML was auto-generated from MATLAB code.
To make changes, update the MATLAB code and republish this document.
      --><title>coolForBisect</title><meta name="generator" content="MATLAB 8.1"><link rel="schema.DC" href="http://purl.org/dc/elements/1.1/"><meta name="DC.date" content="2014-04-15"><meta name="DC.source" content="coolForBisect.m"><style type="text/css">
html,body,div,span,applet,object,iframe,h1,h2,h3,h4,h5,h6,p,blockquote,pre,a,abbr,acronym,address,big,cite,code,del,dfn,em,font,img,ins,kbd,q,s,samp,small,strike,strong,sub,sup,tt,var,b,u,i,center,dl,dt,dd,ol,ul,li,fieldset,form,label,legend,table,caption,tbody,tfoot,thead,tr,th,td{margin:0;padding:0;border:0;outline:0;font-size:100%;vertical-align:baseline;background:transparent}body{line-height:1}ol,ul{list-style:none}blockquote,q{quotes:none}blockquote:before,blockquote:after,q:before,q:after{content:'';content:none}:focus{outine:0}ins{text-decoration:none}del{text-decoration:line-through}table{border-collapse:collapse;border-spacing:0}

html { min-height:100%; margin-bottom:1px; }
html body { height:100%; margin:0px; font-family:Arial, Helvetica, sans-serif; font-size:10px; color:#000; line-height:140%; background:#fff none; overflow-y:scroll; }
html body td { vertical-align:top; text-align:left; }

h1 { padding:0px; margin:0px 0px 25px; font-family:Arial, Helvetica, sans-serif; font-size:1.5em; color:#d55000; line-height:100%; font-weight:normal; }
h2 { padding:0px; margin:0px 0px 8px; font-family:Arial, Helvetica, sans-serif; font-size:1.2em; color:#000; font-weight:bold; line-height:140%; border-bottom:1px solid #d6d4d4; display:block; }
h3 { padding:0px; margin:0px 0px 5px; font-family:Arial, Helvetica, sans-serif; font-size:1.1em; color:#000; font-weight:bold; line-height:140%; }

a { color:#005fce; text-decoration:none; }
a:hover { color:#005fce; text-decoration:underline; }
a:visited { color:#004aa0; text-decoration:none; }

p { padding:0px; margin:0px 0px 20px; }
img { padding:0px; margin:0px 0px 20px; border:none; }
p img, pre img, tt img, li img { margin-bottom:0px; } 

ul { padding:0px; margin:0px 0px 20px 23px; list-style:square; }
ul li { padding:0px; margin:0px 0px 7px 0px; }
ul li ul { padding:5px 0px 0px; margin:0px 0px 7px 23px; }
ul li ol li { list-style:decimal; }
ol { padding:0px; margin:0px 0px 20px 0px; list-style:decimal; }
ol li { padding:0px; margin:0px 0px 7px 23px; list-style-type:decimal; }
ol li ol { padding:5px 0px 0px; margin:0px 0px 7px 0px; }
ol li ol li { list-style-type:lower-alpha; }
ol li ul { padding-top:7px; }
ol li ul li { list-style:square; }

.content { font-size:1.2em; line-height:140%; padding: 20px; }

pre, tt, code { font-size:12px; }
pre { margin:0px 0px 20px; }
pre.error { color:red; }
pre.codeinput { padding:10px; border:1px solid #d3d3d3; background:#f7f7f7; }
pre.codeoutput { padding:10px 11px; margin:0px 0px 20px; color:#4c4c4c; }

@media print { pre.codeinput, pre.codeoutput { word-wrap:break-word; width:100%; } }

span.keyword { color:#0000FF }
span.comment { color:#228B22 }
span.string { color:#A020F0 }
span.untermstring { color:#B20000 }
span.syscmd { color:#B28C00 }

.footer { width:auto; padding:10px 0px; margin:25px 0px 0px; border-top:1px dotted #878787; font-size:0.8em; line-height:140%; font-style:italic; color:#878787; text-align:left; float:none; }
.footer p { margin:0px; }
.footer a { color:#878787; }
.footer a:hover { color:#878787; text-decoration:underline; }
.footer a:visited { color:#878787; }

table th { padding:7px 5px; text-align:left; vertical-align:middle; border: 1px solid #d6d4d4; font-weight:bold; }
table td { padding:7px 5px; text-align:left; vertical-align:top; border:1px solid #d6d4d4; }





  </style></head><body><div class="content"><pre class="codeinput"><span class="comment">% Program 8.5 Finite difference solver for 2D Poisson equation</span>
<span class="comment">%   with Dirichlet boundary conditions on a rectangle</span>
<span class="comment">% Input: rectangle domain [xl,xr]x[yb,yt] with MxN space steps</span>


<span class="comment">% Output: matrix w holding solution values</span>
<span class="comment">% Example usage: w=cool(0,4,0,4,10,10,4,1.68,.005,3,3)</span>
<span class="keyword">function</span> w=coolForBisect(xl,xr,yb,yt,M,N,P,K,H,B,C)

<span class="comment">% B (length) and C (width) are the dimensions of the notch we will cut out</span>
<span class="comment">% of the fin</span>

ub = 20; <span class="comment">%degrees Celcius</span>
delta=.1; <span class="comment">%cm thickness of fin</span>

L=2; <span class="comment">%length along which power is input</span>

m=M+1;n=N+1; mn=m*n;
h=(xr-xl)/M;  h2=h^2;      k=(yt-yb)/N;  k2=k^2;
x=xl+(0:M)*h; <span class="comment">% set mesh values</span>
y=yb+(0:N)*k;
A=zeros(mn,mn);      b=zeros(mn,1); <span class="comment">% initially set A and b to zeros</span>
size(A) <span class="comment">%check the size of both A and b to make sure they match</span>
size(b)

<span class="comment">% Interior Points</span>
<span class="keyword">for</span> i=2:m-1
  <span class="keyword">for</span> j=2:n-1
     <span class="keyword">if</span> (i &gt; m-C &amp;&amp; j &lt;= B) <span class="comment">% If we are in the notch, then they cannot be</span>
                            <span class="comment">% interior points. =&gt; We set them to zero.</span>
         A(i+(j-1)*m,i-1+(j-1)*m)=0;       A(i+(j-1)*m,i+1+(j-1)*m)=0;
         A(i+(j-1)*m,i+(j-1)*m)=0;
         A(i+(j-1)*m,i+(j-2)*m)=0;         A(i+(j-1)*m,i+j*m)=0;
         disp(<span class="string">'In the notch'</span>)
     <span class="keyword">else</span>                   <span class="comment">% Otherwise, set up the interior points as</span>
                            <span class="comment">% before.</span>
         A(i+(j-1)*m,i-1+(j-1)*m)=1/h2;       A(i+(j-1)*m,i+1+(j-1)*m)=1/h2;
         A(i+(j-1)*m,i+(j-1)*m)=-2/h2-2/k2-(2*H)/(K*delta);
         A(i+(j-1)*m,i+(j-2)*m)=1/k2;         A(i+(j-1)*m,i+j*m)=1/k2;
         disp(<span class="string">'Regular interior points'</span>)
     <span class="keyword">end</span>
  <span class="keyword">end</span>
       <span class="comment">% disp(i)</span>
<span class="keyword">end</span>

disp(<span class="string">'Part 2'</span>)
<span class="keyword">for</span> i=1:m <span class="comment">% bottom and top boundary points</span>
  <span class="keyword">if</span> i &gt; (m - C) <span class="comment">% If our i gets to the edge of the notch</span>
                 <span class="comment">% move the bottom boundary up to the edge of the notch</span>
      j=1+B;     A(i+(1-1)*m,i+(j-1)*m) = -3/(2*k) + H/K;      <span class="comment">%bottom</span>
      j=2+B;     A(i+(1-1)*m,i+(j-1)*m) = 4/(2*k);
      j=3+B;     A(i+(1-1)*m,i+(j-1)*m) = -1/(2*k);
      disp(<span class="string">'A'</span>)
  <span class="keyword">else</span>
      <span class="comment">% Otherwise, leave it alone.</span>
      j=1;     A(i+(1-1)*m,i+(j-1)*m) = -3/(2*k) + H/K;      <span class="comment">%bottom</span>
      j=2;     A(i+(1-1)*m,i+(j-1)*m) = 4/(2*k);
      j=3;     A(i+(1-1)*m,i+(j-1)*m) = -1/(2*k);
      disp(<span class="string">'Else'</span>)
  <span class="keyword">end</span>
  j=n-2;   A(i+(n-1)*m,i+(j-1)*m) =  -1/(2*k);      <span class="comment">%top</span>
  j=n-1;   A(i+(n-1)*m,i+(j-1)*m) = 4/(2*k);
  j=n;     A(i+(n-1)*m,i+(j-1)*m) =  -3/(2*k) + H/K;


<span class="keyword">end</span>
<span class="keyword">for</span> j=2:n-1 <span class="comment">% left and right boundary points</span>
  i=1;     A(1+(j-1)*m,i+(j-1)*m)=-3/(2*h);      b(i+(j-1)*m)= -P/(L*delta*K); <span class="comment">%left</span>
  i=2;     A(1+(j-1)*m,i+(j-1)*m)=4/(2*h);      <span class="comment">%b(i+(j-1)*m)= -P/(L*delta*K);</span>
  i=3;     A(1+(j-1)*m,i+(j-1)*m)=-1/(2*h);      <span class="comment">%b(i+(j-1)*m)= -P/(L*delta*K);</span>
  <span class="keyword">if</span> j &lt;= B <span class="comment">% If j is within the notch, then move the boundary to the left,</span>
            <span class="comment">% to the left of the notch</span>
      i=m-(C+2);   A(m+(j-1)*m,i+(j-1)*m)= -1/(2*h);     <span class="comment">%right</span>
      i=m-(C+1);   A(m+(j-1)*m,i+(j-1)*m)=4/(2*h);
      i=m-C;     A(m+(j-1)*m,i+(j-1)*m)=  -3/(2*h) + H/K;
      disp(<span class="string">'j &lt; n-(B+1)'</span>)
  <span class="keyword">else</span>
      i=m-2;   A(m+(j-1)*m,i+(j-1)*m)= -1/(2*h);     <span class="comment">%right</span>
      i=m-1;   A(m+(j-1)*m,i+(j-1)*m)=4/(2*h);
      i=m;     A(m+(j-1)*m,i+(j-1)*m)=  -3/(2*h) + H/K;
      disp(<span class="string">'else'</span>)
  <span class="keyword">end</span>
<span class="keyword">end</span>
A;
b;

v=A\b <span class="comment">% solve for solution in v labeling</span>
v = v+ ub;
w=reshape(v(1:mn),m,n); <span class="comment">%translate from v to w</span>
mesh(x,y,w')
</pre><pre class="codeoutput">Error using coolForBisect (line 18)
Not enough input arguments.
</pre><p class="footer"><br><a href="http://www.mathworks.com/products/matlab/">Published with MATLAB&reg; R2013a</a><br></p></div><!--
##### SOURCE BEGIN #####
% Program 8.5 Finite difference solver for 2D Poisson equation 
%   with Dirichlet boundary conditions on a rectangle 
% Input: rectangle domain [xl,xr]x[yb,yt] with MxN space steps 


% Output: matrix w holding solution values 
% Example usage: w=cool(0,4,0,4,10,10,4,1.68,.005,3,3) 
function w=coolForBisect(xl,xr,yb,yt,M,N,P,K,H,B,C)

% B (length) and C (width) are the dimensions of the notch we will cut out 
% of the fin

ub = 20; %degrees Celcius
delta=.1; %cm thickness of fin

L=2; %length along which power is input

m=M+1;n=N+1; mn=m*n; 
h=(xr-xl)/M;  h2=h^2;      k=(yt-yb)/N;  k2=k^2; 
x=xl+(0:M)*h; % set mesh values 
y=yb+(0:N)*k; 
A=zeros(mn,mn);      b=zeros(mn,1); % initially set A and b to zeros
size(A) %check the size of both A and b to make sure they match
size(b)

% Interior Points
for i=2:m-1  
  for j=2:n-1 
     if (i > m-C && j <= B) % If we are in the notch, then they cannot be
                            % interior points. => We set them to zero.
         A(i+(j-1)*m,i-1+(j-1)*m)=0;       A(i+(j-1)*m,i+1+(j-1)*m)=0; 
         A(i+(j-1)*m,i+(j-1)*m)=0;
         A(i+(j-1)*m,i+(j-2)*m)=0;         A(i+(j-1)*m,i+j*m)=0; 
         disp('In the notch')
     else                   % Otherwise, set up the interior points as
                            % before.
         A(i+(j-1)*m,i-1+(j-1)*m)=1/h2;       A(i+(j-1)*m,i+1+(j-1)*m)=1/h2; 
         A(i+(j-1)*m,i+(j-1)*m)=-2/h2-2/k2-(2*H)/(K*delta); 
         A(i+(j-1)*m,i+(j-2)*m)=1/k2;         A(i+(j-1)*m,i+j*m)=1/k2; 
         disp('Regular interior points')
     end
  end 
       % disp(i)
end 

disp('Part 2')
for i=1:m % bottom and top boundary points 
  if i > (m - C) % If our i gets to the edge of the notch
                 % move the bottom boundary up to the edge of the notch
      j=1+B;     A(i+(1-1)*m,i+(j-1)*m) = -3/(2*k) + H/K;      %bottom
      j=2+B;     A(i+(1-1)*m,i+(j-1)*m) = 4/(2*k);
      j=3+B;     A(i+(1-1)*m,i+(j-1)*m) = -1/(2*k);
      disp('A')
  else
      % Otherwise, leave it alone. 
      j=1;     A(i+(1-1)*m,i+(j-1)*m) = -3/(2*k) + H/K;      %bottom
      j=2;     A(i+(1-1)*m,i+(j-1)*m) = 4/(2*k);
      j=3;     A(i+(1-1)*m,i+(j-1)*m) = -1/(2*k);
      disp('Else')
  end
  j=n-2;   A(i+(n-1)*m,i+(j-1)*m) =  -1/(2*k);      %top
  j=n-1;   A(i+(n-1)*m,i+(j-1)*m) = 4/(2*k);
  j=n;     A(i+(n-1)*m,i+(j-1)*m) =  -3/(2*k) + H/K;
  
  
end 
for j=2:n-1 % left and right boundary points 
  i=1;     A(1+(j-1)*m,i+(j-1)*m)=-3/(2*h);      b(i+(j-1)*m)= -P/(L*delta*K); %left
  i=2;     A(1+(j-1)*m,i+(j-1)*m)=4/(2*h);      %b(i+(j-1)*m)= -P/(L*delta*K);
  i=3;     A(1+(j-1)*m,i+(j-1)*m)=-1/(2*h);      %b(i+(j-1)*m)= -P/(L*delta*K);
  if j <= B % If j is within the notch, then move the boundary to the left,
            % to the left of the notch
      i=m-(C+2);   A(m+(j-1)*m,i+(j-1)*m)= -1/(2*h);     %right
      i=m-(C+1);   A(m+(j-1)*m,i+(j-1)*m)=4/(2*h);       
      i=m-C;     A(m+(j-1)*m,i+(j-1)*m)=  -3/(2*h) + H/K;  
      disp('j < n-(B+1)')
  else
      i=m-2;   A(m+(j-1)*m,i+(j-1)*m)= -1/(2*h);     %right
      i=m-1;   A(m+(j-1)*m,i+(j-1)*m)=4/(2*h);       
      i=m;     A(m+(j-1)*m,i+(j-1)*m)=  -3/(2*h) + H/K;
      disp('else')
  end
end 
A;
b;

v=A\b % solve for solution in v labeling 
v = v+ ub;
w=reshape(v(1:mn),m,n); %translate from v to w 
mesh(x,y,w')
##### SOURCE END #####
--></body></html>