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. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top