#### /chapter12.m

https://github.com/rwinklerwilkes/Matlab
Objective C | 368 lines | 367 code | 1 blank | 0 comment | 30 complexity | 5469064d5a6c4da8f68b03e313aff046 MD5 | raw file
```  1%Poisson Equation Finite-Difference
2%Variable key
3%a - Lowest x-value
4%b - Highest x-value
5%c - Lowest y-value
6%d - Highest y-value
7%m - Number of y partitions, >=3
8%n - Number of x partitions, >=3
9%f - Function they satisfy
10%gB - Boundary condition for y=c
11%gT - Boundary condition for y=d
12%gL - Boundary condition for x=a
13%gR - Boundary condition for x=b
14function p = pefd(a,b,c,d,m,n,f,gB,gT,gL,gR,TOL,maxIter)
15	%Step 1
16	h = (b-a)/n;
17	k = (d-c)/m;
18	x = zeros(1,n-1);
19	y = zeros(1,m-1);
20	%Step 4
21	w = zeros(n,m);
22	%Step 2
23	i=1;
24	while(i<=(n-1))
25		x(i) = a+i*h;
26		i++;
27	endwhile
28	%Step 3
29	i = 1;
30	while(i<=(m-1))
31		y(i) = c+i*k;
32		i++;
33	endwhile
34	%Step 5
35	lambda = (h**2)/(k**2);
36	mu = 2*(1+lambda);
37	l = 1;
38	%Step 6
39	while(l<=maxIter)
40		%Step 7
41		z = ((-h**2)*f(x(1),y(m-1)) +gL(y(m-1)) + lambda*gT(x(1)) + lambda*w(1,m-2) + w(2,m-1))/mu;
42		NORM = abs(z-w(1,m-1));
43		w(1,m-1) = z;
44		%Step 8
45		i = 2;
46		while(i<=(n-2))
47			z = ((-h**2)*f(x(i),y(m-1)) + lambda*gT(x(i)) + w(i-1,m-1) + w(i+1,m-1) + lambda*w(i,m-2))/mu;
48			if(abs(w(i,m-1) - z) > NORM)
49				NORM = abs(w(i,m-1)-z);
50			endif
51			w(i,m-1) = z;
52			i++;
53		endwhile
54		%Step 9
55		z = ((-h**2)*f(x(n-1),y(m-1)) + gR(y(m-1)) + lambda*gT(x(n-1)) + w(n-2,m-1) + lambda*w(n-1,m-2))/mu;
56		if(abs(w(n-1,m-1)-z) > NORM)
57			NORM = abs(w(n-1,m-1)-z);
58		endif
59		w(n-1,m-1) = z;
60		%Step 10
61		j = m-2;
62		while(j>=2)
63			%Step 11
64			z = ((-h**2)*f(x(1),y(j)) + gL(y(j)) + lambda*w(1,j+1)+lambda*w(1,j-1)+w(2,j))/mu;
65			if(abs(w(1,j)-z) > NORM)
66				NORM = abs(w(1,j)-z);
67			endif
68			w(1,j) = z;
69			%Step 12
70			i = 2;
71			while(i<=n-2)
72				z = ((-h**2)*f(x(i),y(j)) + w(i-1,j) + lambda*w(i,j+1)+w(i+1,j)+lambda*w(i,j-1))/mu;
73				if(abs(w(i,j)-z) > NORM)
74					NORM = abs(w(i,j)-z);
75				endif
76				w(i,j) = z;
77				i++;
78			endwhile
79			%Step 13
80			z= (-(h**2)*f(x(n-1),y(j)) + gR(y(j)) + w(n-2,j) + lambda*w(n-1,j+1) + lambda*w(n-1,j-1))/mu;
81			if(abs(w(n-1,j)-z) > NORM)
82				NORM = abs(w(n-1,j)-z);
83			endif
84			w(n-1,j) = z;
85			j--;
86		endwhile
87		%Step 14
88		z= (-(h**2)*f(x(1),y(1)) + gL(y(1))+lambda*gB(x(1)) + lambda*w(1,2) + w(2,1))/mu;
89		if(abs(w(1,1)-z) > NORM)
90			NORM = abs(w(1,1)-z);
91		endif
92		w(1,1) = z;
93		%Step 15
94		i = 2;
95		while(i<=(n-2))
96			z= (-(h**2)*f(x(i),y(1)) + lambda*gB(x(i))+w(i-1,1) + lambda*w(i,2) + w(i+1,1))/mu;
97			if(abs(w(i,1)-z) > NORM)
98				NORM = abs(w(i,1)-z);
99			endif
100			w(i,1) = z;
101			i++;
102		endwhile
103		%Step 16
104		z = (-(h**2)*f(x(n-1),y(1)) + gR(y(1)) + lambda*gB(x(n-1)) + w(n-2,1) + lambda*w(n-1,2))/mu;
105		if(abs(w(n-1,1) - z) > NORM)
106			NORM = abs(w(n-1,1) - z);
107		endif
108		w(n-1,1) = z;
109		%Step 17
110		if(NORM<=TOL)
111			i = 1;
112			while(i<=(n-1))
113				j = 1;
114				while(j<=(m-1))
115					printf("x(%d): %g, y(%d): %g, w(%d,%d): %g\n", i,x(i),j,y(j),i,j,w(i,j));
116					j++;
117				endwhile
118				i++;
119			endwhile
120			p = 0;
121			return;
122		endif
123		l++;
124	endwhile
125	printf("Max iterations exceeded");
126	p = -1;
127	return;
128endfunction
129
130%g - Dirichlet boundary condition
131%g1,g2 - Neumann boundary conditions
132%verts - Vertices of the triangulation
133%nodes - each node corresponds to the row of verts which contains the vertex for that node
134%triVerts - each row has 3 values - the three node numbers of the triangle
135%K - Interior triangles
136%N - Triangles with at least one Neumann edge
137%M - All other triangles (total triangles - K+N+M)
138%n - Nodes either in D or on the Neumann edge
139%m - Dirichlet nodes
140function f = fe(p,q,r,f,g,g1,g2,verts,nodes,triVerts,K,N,M,n,m)
141	gamma = zeros(1,m);
142	l = n+1;
143	while(l<=m)
144		gamma(l) = g(verts(nodes(l),1),verts(nodes(l),2));
145		l++;
146	endwhile
147	alpha = zeros(n,n);
148	beta = zeros(1,n);
149	delta = zeros(1,M);
150	a = zeros(M,3);
151	b = zeros(M,3);
152	c = zeros(M,3);
153	z = zeros(M,3,3);
154	J = zeros(M,3,3);
155	I = zeros(M,3);
156	H = zeros(M,3);
157	i = 1;
158	while(i<=M)
159		x1 = verts(nodes(triVerts(i,1)),1);
160		x2 = verts(nodes(triVerts(i,2)),1);
161		x3 = verts(nodes(triVerts(i,3)),1);
162		y1 = verts(nodes(triVerts(i,1)),2);
163		y2 = verts(nodes(triVerts(i,2)),2);
164		y3 = verts(nodes(triVerts(i,3)),2);
165		delta(i) = det([1,x1,y1;1,x2,y2;1,x3,y3]);
166		a(i,1) = (x2*y3-y2*x3)/delta(i);
167		a(i,2) = (x3*y1-y3*x1)/delta(i);
168		a(i,3) = (x1*y2-y1*x2)/delta(i);
169		b(i,1) = (y2-y3)/delta(i);
170		b(i,2) = (y3-y1)/delta(i);
171		b(i,3) = (y1-y2)/delta(i);
172		c(i,1) = (x3-x2)/delta(i);
173		c(i,2) = (x1-x3)/delta(i);
174		c(i,3) = (x2-x1)/delta(i);
175		i++;
176	endwhile
177	i = 1;
178	while(i<=M)
179		j = 1;
180		while(j<=3)
181			k = 1;
182			while(k<=j)
183				nj = @(x,y)(a(i,j) + b(i,j)*x + c(i,j)*y);
184				nk = @(x,y)(a(i,k) + b(i,k)*x + c(i,k)*y);
185				%Calculate the integrals over the triangles
186				zIntegral = @(x,y)(r(x,y)*nj(x,y)*nk(x,y));
187				%Approximate the integral using quadrature rule for triangles
188				%node 1
189				x1 = verts(nodes(triVerts(i,1)),1);
190				y1 = verts(nodes(triVerts(i,1)),2);
191				%node 2
192				x2 = verts(nodes(triVerts(i,2)),1);
193				y2 = verts(nodes(triVerts(i,2)),2);
194				%node 3
195				x3 = verts(nodes(triVerts(i,3)),1);
196				y3 = verts(nodes(triVerts(i,3)),2);
197				area = (1/6)*abs(det([x1,x2,x3;y1,y2,y3;1,1,1]));
198				%mid1 = midpoint between node 1 and node 2
199				midx1 = (x1+x2)/2;
200				midy1 = (y1+y2)/2;
201				%mid1 = midpoint between node 2 and node 3
202				midx2 = (x2+x3)/2;
203				midy2 = (y2+y3)/2;
204				%mid1 = midpoint between node 1 and node 3
205				midx3 = (x1+x3)/2;
206				midy3 = (y1+y3)/2;
207				firstTerm = b(i,j)*b(i,k)*area*(p(x1,y1)+p(x2,y2)+p(x3,y3));
208				secondTerm = c(i,j)*c(i,k)*area*(q(x1,y1)+q(x2,y2)+q(x3,y3));
209				thirdTerm = -1*area*(zIntegral(x1,y1)+zIntegral(x2,y2)+zIntegral(x3,y3));
210				z(i,j,k) = firstTerm+secondTerm+thirdTerm;
211				hInt = @(x,y)(f(x,y)*nj(x,y));
212				H(i,j) = -1*area*(hInt(x1,y1)+hInt(x2,y2)+hInt(x3,y3));
213				k++;
214			endwhile
215			j++;
216		endwhile
217		i++;
218	endwhile
219	%Approximate the line integrals for the Neumann nodes
220	%%%%%%%%%%%%%%%%%%%%%%%%%%%%
221	%Unfinished, probably doesn't work
222	i = K+1;
223	while(i<N)
224		j = 1;
225		while(j<=3)
226			k = 1;
227			while(k<=j)
228				node1X = verts(nodes(i),1);
229				node1Y = verts(nodes(i),2);
230				node2X = verts(nodes(i+1),1);
231				node2Y = verts(nodes(i+1),2);
232				length = sqrt((node2X-node1X)**2+(node2Y-node1Y)**2);
233				nj = @(x,y)(a(i,j) + b(i,j)*x + c(i,j)*y);
234				nk = @(x,y)(a(i,k) + b(i,k)*x + c(i,k)*y);
235				firstInt = @(x,y)(g1(x,y)*nj(x,y)*nk(x,y));
236				secondInt = @(x,y)(g2(x,y)*nj(x,y));
237				J(i,j,k) = (length/4)*(firstInt(node1X,node1Y)+firstInt(node2X,node2Y));
238				I(i,j) = (length/4)*(secondInt(node1X,node1Y)+secondInt(node2X,node2Y));
239				k++;
240			endwhile
241			j++;
242		endwhile
243		i++;
244	endwhile
245	%%%%%%%%%%%%%%%%%%%%%%%%%%%%
246	%Step 7
247	i = 1;
248	while(i<=M)
249		k = 1;
250		while(k<=3)
251			%Step 8
252			triNodeIKX = verts(nodes(triVerts(i,k)),1);
253			triNodeIKY = verts(nodes(triVerts(i,k)),2);
254			l = 1;
255			while(l<=M)
256				if((verts(nodes(l),1)==triNodeIKX)&&(verts(nodes(l),2)==triNodeIKY))
257					break;
258				endif
259				l++;
260			endwhile
261			%Step 9
262			if(k>1)
263				j = 1;
264				while(j<=k-1)
265					%Step 10
266					triNodeIKX = verts(nodes(triVerts(i,j)),1);
267					triNodeIKY = verts(nodes(triVerts(i,j)),2);
268					t = 1;
269					while(t<=M)
270						if((verts(nodes(t),1)==triNodeIKX)&&(verts(nodes(t),2)==triNodeIKY))
271							break;
272						endif
273						t++;
274					endwhile
275					%Step 11
276					if(l<=n)
277						if(t<=n)
278							alpha(l,t) = alpha(l,t)+z(i,k,j);
279							alpha(t,l) = alpha(t,l)+z(i,k,j);
280						else
281							beta(l) = beta(l) - gamma(t)*z(i,k,j);
282						endif
283					else
284						if(t<=n)
285							beta(t) = beta(t) - gamma(l)*z(i,k,j);
286						endif
287					endif
288					j++;
289				endwhile
290			endif
291			if(l<=n)
292				alpha(l,l) = alpha(l,l) + z(i,k,k);
293				beta(l) = beta(l) + H(i,k);
294			endif
295			k++;
296		endwhile
297		i++;
298	endwhile
299	%Step 13
300	i = K+1;
301	while(i<=N)
302		%Step 14
303		k = 1;
304		while(k<=3)
305			%Step 15
306			triNodeIKX = verts(nodes(triVerts(i,k)),1);
307			triNodeIKY = verts(nodes(triVerts(i,k)),2);
308			l = 1;
309			while(l<=M)
310				if((verts(nodes(l),1)==triNodeIKX)&&(verts(nodes(l),2)==triNodeIKY))
311					break;
312				endif
313				l++;
314			endwhile
315			%Step 16
316			if(k>1)
317				j = 1;
318				while(j<=k-1)
319					%Step 17
320					triNodeIKX = verts(nodes(triVerts(i,j)),1);
321					triNodeIKY = verts(nodes(triVerts(i,j)),2);
322					t = 1;
323					while(t<=M)
324						if((verts(nodes(t),1)==triNodeIKX)&&(verts(nodes(t),2)==triNodeIKY))
325							break;
326						endif
327						t++;
328					endwhile
329					%Step 18
330					if(l<=n)
331						if(t<=n)
332							alpha(l,t) = alpha(l,t) + J(i,k,j);
333							alpha(t,l) = alpha(t,l) + J(i,k,j);
334						else
335							beta(l) = beta(l) - gamma(t)*J(i,k,j);
336						endif
337					else
338						if(t<=n)
339							beta(t) = beta(t)-gamma(l)*J(i,k,j);
340						endif
341					endif
342					j++;
343				endwhile
344			endif
345			%Step 19
346			if(l<=n)
347				alpha(l,l) = alpha(l,l) + J(i,k,k);
348				beta(l) = beta(l) + I(i,k);
349			endif
350			k++;
351		endwhile
352		i++;
353	endwhile
354	disp(alpha);
355	disp(beta');
356	solution = alpha \ beta';
357	disp(solution);
358	i = 1;
359	while(i<=M)
360		j = 1;
361		while(j<=3)
362			printf("Triangle %i, a(%i): %4g, b(%i): %4g, c(%i): %4g\n",i,j,a(i,j),j,b(i,j),j,c(i,j));
363			j++;
364		endwhile
365		i++;
366	endwhile
367	return;
368endfunction```