SHARE
TWEET

Untitled

a guest Apr 23rd, 2019 98 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. clear();
  2. a = 0;
  3. b = 2;
  4. T = 1;
  5. Nx = 200;
  6. Nt = 400;
  7. sigma = 0.5;
  8.  
  9. function y = u_exact(x,t)
  10.     y = t*t*exp(x);
  11. endfunction
  12.  
  13. function y=f(x,t)
  14.     y= (2*t-t*t)*exp(x);
  15. endfunction
  16.  
  17. function y=err(Nx,Nt,u_ex)
  18.     y = abs( solve(Nx,Nt) - u_ex );
  19. endfunction
  20.  
  21. function y=progonka(A,B,C,D)
  22.     // A,B,C,D - массивы поддиагональных эл-в, диагональных, наддиагональных и правой части соответственно
  23.     N = length(B);
  24.     c(1) = C(1)/B(1);
  25.     d(1) = D(1)/B(1);
  26.     for i = 2:N-1
  27.         k = B(i) - A(i)*c(i-1);
  28.         c(i) = C(i)/k;
  29.         d(i) = (D(i) - A(i)*d(i-1))/k;
  30.     end
  31.     y(N) = (D(N) - A(N)*d(N-1))/(B(N) - A(N)*c(N-1));
  32.     for i = N-1:-1:1
  33.         y(i) = d(i) - c(i)*y(i+1);
  34.     end
  35. endfunction
  36.  
  37. //u_0 = 0 в нашем случае
  38. function y=solve(Nx,Nt)
  39.     h = (b-a)/Nx;
  40.     tau = T/Nt;
  41.     x = a:h:b;
  42.     t = 0:tau:T;
  43.    
  44.     U = zeros(Nt+1,Nx+1);
  45.     A = zeros(Nx+1,1);
  46.     B = zeros(Nx+1,1);
  47.     C = zeros(Nx+1,1);
  48.     D = zeros(Nx+1,1);
  49.    
  50.     //цикл по t
  51.     for n = 1:Nt
  52.        
  53.        
  54.     //Заполняем диагонали
  55.     for i = 2:Nx
  56.         B(i) = -h*h-2*tau*sigma;
  57.         A(i) = tau*sigma;
  58.         C(i) = tau*sigma;
  59.         D(i) = (-h*h+2*tau*(1-sigma))*U(n,i)-tau*(1-sigma)*(U(n,i+1)+U(n,i-1))-tau*h*h*f(x(i),t(n));
  60.     end
  61.    
  62. // Граничные условия
  63. C(1) = 4 + B(2)/C(2);
  64. B(1) = -3 + A(2)/C(2);
  65. D(1) = 2*t(n+1)*t(n+1)*h + D(2)/C(2);
  66.  
  67. B(Nx+1) = 3  - C(Nx)/A(Nx);
  68. A(Nx+1) = -4 - B(Nx)/A(Nx);
  69. D(Nx+1) = 2*t(n+1)*t(n+1)*h*exp(b) - D(Nx)/A(Nx);  
  70.  
  71. U(n+1,:) = progonka(A,B,C,D)';
  72.  
  73. end
  74. y = U;
  75. endfunction
  76.  
  77.  
  78. U = solve(Nx,Nt);
  79.  
  80.     h = (b-a)/Nx;
  81.     tau = T/Nt;
  82.     x = a:h:b;
  83.     t = 0:tau:T;
  84.    
  85.    
  86.  
  87.  
  88. for i=1:Nx+1
  89.     for n=1:Nt+1
  90.         U_ex(n,i) = u_exact(x(i),t(n));
  91.         end
  92. end
  93.  
  94.  
  95. f1=scf(1);  //creates figure with id==4 and make it the current one
  96. f1.color_map = jetcolormap(64);
  97. xtitle ('x', 't', 'exact_solution');
  98. plot3d1(t,x,U_ex);
  99.  
  100. f2=scf(2);  //creates figure with id==4 and make it the current one
  101. f2.color_map = jetcolormap(64);
  102. xtitle ('x', 't', 'U');
  103. plot3d1(t,x,U);
  104.  
  105. f3=scf(3);  //creates figure with id==4 and make it the current one
  106. f3.color_map = jetcolormap(64);
  107. xtitle ('x', 't', 'error');
  108. mis=max(err(Nx,Nt,U_ex));
  109. plot3d1(t,x,err(Nx,Nt,U_ex));
  110. /*
  111. // pogresnosti
  112. NxL = [50, 100,200,400,800];
  113. NtL = [25,100,400,1600,6400];
  114. errL = zeros(5);
  115.  
  116. for j = 1:5
  117.     Nx = NxL(j);
  118.     Nt = NtL(j);
  119.  
  120.     h = (b-a)/Nx;
  121.     tau = T/Nt;
  122.     x = a:h:b;
  123.     t = 0:tau:T;
  124.    
  125.     U_ex = zeros(Nt+1,Nx+1);
  126.  
  127.  
  128. for i=1:Nx+1
  129.     for n=1:Nt+1
  130.         U_ex(n,i) = u_exact(x(i),t(n));
  131.         end
  132. end
  133.     errL(j) = max(err(Nx,Nt,U_ex));
  134. end
  135.  
  136. */
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