SHARE
TWEET

Untitled

a guest Apr 23rd, 2019 75 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. d = 70;
  2. M = d-1;
  3. N = d-1;
  4.  
  5. %solve Laplace equation
  6. % nabla^2 u(x,y) = 0
  7.  
  8. % Specification of rectangular area
  9. x0 = 0; xf = 4; % minimal and maximal value of x-coordinate
  10. Mx = d; % number of intervals in x direction
  11. y0 = 0; yf = 4; % minimal and maximal value of y-coordinate
  12. My = d; % number of intervals in y direction
  13.  
  14. bx0 = inline('exp(y) - cos(y)','y'); % Boundary condition for left edge
  15. bxf = inline('exp(y)*cos(4) - exp(4)*cos(y)','y'); % Boundary condition for right edge
  16. by0 = inline('cos(x) - exp(x)','x'); % Boundary condition for bottom edge
  17. byf = inline('exp(4)*cos(x) - exp(x)*cos(4)','x'); % Boundary condition for top edge
  18.  
  19. % a vector for defining rectangular area
  20. D = [x0 xf y0 yf];
  21.  
  22. %Parameters for Jacobi Method
  23. MaxIter = 1000; % Maximal number of iteration
  24. tol = 1e-3; % tolerance of Jacobi method
  25.  
  26. % Solving Laplace equation via function sollap
  27. [Ua,x,y] = sollap(bx0,bxf,by0,byf,D,Mx,My,tol,MaxIter);
  28.  
  29. Un = zeros(d+1);
  30.  
  31. for i = 1:length(x)
  32.     for j = 1:length(y)
  33.         Un(i,j) = -((exp(y(j)) * cos(x(i))) - (exp(x(i)) * cos(y(j))));
  34.     end
  35. end
  36.  
  37. Ud = Un - Ua;
  38.  
  39. Ud2 = Ud.^2;
  40.  
  41. totals = zeros(1, length(Ud2));
  42.  
  43. for i = 1:length(Ud2)
  44.     totals(i) = sum(Ud2(i,:));
  45. end
  46.  
  47. USum = sum(totals);
  48.  
  49. error = sqrt((1/(M*N)) * USum);
  50.  
  51.  
  52.  
  53.  
  54. In new script thing....
  55. d=[5 10 15 20 25 30 35 40 45 50 55 60 65 70]
  56. Error= [0.2459297 0.057170249 0.019832251 0.008179454 0.016384573 0.028631307 0.041942391 0.05724528 0.073729247 0.091975748 0.111647881 0.176497771 0.294126947 0.440691868]
  57.  
  58.  
  59. In new script....
  60. zz=zeros(1,70);
  61.  
  62.  
  63. for i=1:1:70
  64.     zz(i) = ((3.424e-07*i^4) -(5.145e-05*i^3) + (0.002808*i^2) -(0.06226*i) + 0.4775)
  65. end
  66.  
  67. figure(1)
  68.  
  69. plot(zz,':b')
  70.  
  71. hold on
  72. box on
  73. plot(d,Error, 'or')
  74. xlabel('Grid Size, d')
  75. ylabel('Error')
  76. title('Error between the Analytical and Numerical solution')
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
 
Top