Advertisement
Guest User

Untitled

a guest
Apr 23rd, 2019
135
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Scilab 1.94 KB | None | 0 0
  1.     function y=TrMxAlg(A,B,C,D) //Прогонка
  2.     N = length(B);
  3.     c(1) = C(1)/B(1);
  4.     d(1) = D(1)/B(1);
  5.     for i = 2:N-1
  6.         k = B(i) - A(i)*c(i-1);
  7.         c(i) = C(i)/k;
  8.         d(i) = (D(i) - A(i)*d(i-1))/k;
  9.     end
  10.     y(N) = (D(N) - A(N)*d(N-1))/(B(N) - A(N)*c(N-1));
  11.     for i = N-1:-1:1
  12.         y(i) = d(i) - c(i)*y(i+1);
  13.     end
  14.     endfunction
  15.  
  16.     function y = u(x,t)//точное решение
  17.     y = t*t*exp(x);
  18.     endfunction
  19.  
  20.     function y=f(x,t)
  21.     y= (2*t-t*t)*exp(x); //функция f в нашем решении
  22.     endfunction
  23.  
  24.     function y=solve(Nx,Nt)
  25.     h = (b-a)/Nx;
  26.     tau = T/Nt;
  27.     x = a:h:b;
  28.     t = 0:tau:T;
  29.     A = zeros(Nx+1,1);
  30.     B = zeros(Nx+1,1);
  31.     C = zeros(Nx+1,1);
  32.     D = zeros(Nx+1,1);
  33.     W = zeros(Nt+1,Nx+1);  
  34.     for n = 1:Nt    
  35.         for i = 2:Nx
  36.         B(i) = -h*h-2*tau*sigma;
  37.         A(i) = tau*sigma;
  38.         C(i) = tau*sigma;
  39.         D(i) = (-h*h+2*tau*(1-sigma))*W(n,i)-tau*(1-sigma)*(W(n,i+1)+W(n,i-1))-tau*h*h*f(x(i),t(n));
  40.         end
  41.     C(1) = 4 + B(2)/C(2);
  42.     B(1) = -3 + A(2)/C(2);
  43.     D(1) = 2*t(n+1)*t(n+1)*h + D(2)/C(2);
  44.     B(Nx+1) = 3  - C(Nx)/A(Nx);
  45.     A(Nx+1) = -4 - B(Nx)/A(Nx);
  46.     D(Nx+1) = 2*t(n+1)*t(n+1)*h*exp(b) - D(Nx)/A(Nx);  
  47.     W(n+1,:) = TrMxAlg(A,B,C,D)';
  48.     end
  49.     y = W;
  50. endfunction
  51. //W-Приближенное рещение
  52. //U-Точное
  53.     a = 0;b = 2;
  54.     T = 1;
  55.     Nx = 200; Nt = 400;
  56.     sigma = 0.5;
  57.     W = solve(Nx,Nt);
  58.     h = (b-a)/Nx;
  59.     tau = T/Nt;
  60.     x = a:h:b;
  61.     t = 0:tau:T;
  62. for i=1:Nx+1
  63.     for n=1:Nt+1
  64.         U(n,i) = u(x(i),t(n));
  65.         end
  66. end
  67.  
  68. Epsilon=max(abs( W - U )); //
  69.  
  70.  
  71. Pucture1=scf(1);
  72.  
  73. xtitle ('x', 't', 'Точное решение');
  74.  
  75.  
  76. plot3d1(t,x,U);
  77.  
  78. Picture2=scf(2);
  79.  
  80. xtitle ('x', 't', 'Приближенное');
  81.  
  82.  
  83. plot3d1(t,x,W);
  84.  
  85. Picture3=scf(3);
  86.  
  87. xtitle ('x', 't', 'Погрешность');
  88.  
  89. plot3d1(t,x,abs(W - U));
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement