• API
• FAQ
• Tools
• Archive
SHARE
TWEET # Untitled a guest Apr 18th, 2019 97 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. for k=20:1:100
2.     % anylitical solution
3.     x=(0:4/k:4); y=(0:4/k:4);
4. for i=1:length(x); j=1:length(y);
5.     u_a(i,j)=exp(y(j))*cos(x(i))-exp(x(i))*cos(y(j));
6. end
7. x0 = 0; xf = 4; % minimal and maximal value of x-coordinate
8. Mx = k; % number of intervals in x direction
9. y0 = 0; yf = 4; % minimal and maximal value of y-coordinate
10. My = k; % number of intervals in y direction
11.
12. bx0 = inline('exp(y) - cos(y)','y'); % Boundary condition for left edge
13. bxf = inline('exp(y)*cos(4) - exp(4)*cos(y)','y'); % Boundary condition for right edge
14. by0 = inline('cos(x) - exp(x)','x'); % Boundary condition for bottom edge
15. byf = inline('exp(4)*cos(x) - exp(x)*cos(4)','x'); % Boundary condition for top edge
16.
17. D = [x0 xf y0 yf];
18.
19. MaxIter = 50000; % Maximal number of iteration
20. tol = 1e-6; % tollerance of Jacobi method
21.
22. % Solving Laplace equation via function sollap
23. [u_n,x,y] = sollap(bx0,bxf,by0,byf,D,Mx,My,tol,MaxIter);
24.
25. %RMS Error
26. U=-u_n-u_a;
27. U_m=mean(U);
28. E(k)=rms(U_m);
29. h(k)=4/k;
30. end
31.
32. plot(h,E,'*')
33. xlabel('Grid Spacing (h)')
34. ylabel('Error')
35. grid on
36.
37.
38. function [u_n,x,y] = sollap(bx0,bxf,by0,byf,D,Mx,My,tol,MaxIter)
39. x0 = D(1); xf = D(2); y0 = D(3); yf = D(4); % end points for each coordinate
40. dx = (xf - x0)/Mx; x = x0 + [0:Mx]*dx; % step size and grid for x-coordinate
41. dy = (yf - y0)/My; y = y0 + [0:My]'*dy; % step size and grid for y-coordinate
42. Mx1 = Mx + 1; My1 = My + 1; % number of nodes for x and y directions
43. %boundary conditions
44. for m = 1:My1, u0(m,[1 Mx1])=[bx0(y(m)) bxf(y(m))]; end % left and right sides
45. for n = 1:Mx1, u0([1 My1],n) = [by0(x(n)); byf(x(n))]; end % bottom and top sides
46. u_n=u0;
47. sum_of_bv = sum(sum([u0(2:My,[1 Mx1]) u0([1 My1],2:Mx)']));
48. u0(2:My,2:Mx) = sum_of_bv/(2*(Mx + My - 2));
49. dx2 = dx*dx; dy2 = dy*dy; dxy2 = 2*(dx2 + dy2);
50. rx = dx2/dxy2; ry = dy2/dxy2;
51.     for itr = 1:MaxIter
52.     for j = 2:Mx
53.     for i = 2:My
54. u_n(i,j) = ry*(u0(i,j + 1)+u0(i,j - 1)) + rx*(u0(i + 1,j)+u0(i - 1,j)) ; % Iterative equation
55.     end
56.     end
57. if itr > 1 & max(max(abs(u_n - u0))) < tol, break; end % stopping condition
58.
59. u0 = u_n;
60.     end
61. end
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!

Top