SHARE
TWEET

Untitled

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