• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Untitled

a guest Apr 24th, 2019 108 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. clear()
2.
3. function y=TrMxAlg(A,B,C,D) //Прогонка
4.     N = length(B);
5.     c(1) = C(1)/B(1);
6.     d(1) = D(1)/B(1);
7.     for i = 2:N-1
8.         k = B(i) - A(i)*c(i-1);
9.         c(i) = C(i)/k;
10.         d(i) = (D(i) - A(i)*d(i-1))/k;
11.     end
12.     y(N) = (D(N) - A(N)*d(N-1))/(B(N) - A(N)*c(N-1));
13.     for i = N-1:-1:1
14.         y(i) = d(i) - c(i)*y(i+1);
15.     end
16.     endfunction
17.
18.     function y = u(x,t)//точное решение
19.     y = t*t*exp(x);
20.     endfunction
21.
22.     function y=f(x,t)
23.     y = (2*t-t*t)*exp(x); //функция f в нашем решении
24.     endfunction
25.     function y=g(x,t)
26.     y = (2*t-t*t)*exp(x)+tau/2*(2-2*t)*exp(x)+(h*h)/12*(2*t-t*t)*exp(x);
27.     endfunction
28.
29.
30.
31.
32.     function y=solve(Nx,Nt)
33.     h = (b-a)/Nx;
34.     tau = T/Nt;
35.     x = a:h:b;
36.     t = 0:tau:T;
37.     A = zeros(Nx+1,1);
38.     B = zeros(Nx+1,1);
39.     C = zeros(Nx+1,1);
40.     D = zeros(Nx+1,1);
41.     W = zeros(Nt+1,Nx+1);
42.     for n = 1:Nt
43.         for i = 2:Nx
44.         B(i) = -h*h-2*tau*sigma;
45.         A(i) = tau*sigma;
46.         C(i) = tau*sigma;
47.         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*g(x(i),t(n));
48.     //  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),n*tau+tau/2)+h*h/12*(f(x(i+1),n*tau+tau/2)-2*f(x(i),n*tau+tau/2)+f(x(i-1),n*tau+tau/2))/(h*h));
49.         end
50.         S1 = [-25 , 48,-36, 16, - 3, 12*t(n+1)*t(n+1)*h];
51.         S2 = [3, - 16, 36, -48, 25 ,12*t(n+1)*t(n+1)*h*exp(b)] ;
52.    for j = 1:3
53.         R1 = [A(2),B(2),C(2),D(1+4-j)];
54.         S1(6) = S1(6) - S1(6-j)*R1(4)/R1(3);
55.         S1(4-j:5-j) = S1(4-j:5-j) - S1(6-j)*R1(1:2)/R1(3);
56.
57.
58.         R2 = [A(2),B(2),C(2),D(Nx+1-(4-j))];
59.         S2(6) = S2(6) - S2(j)*R2(4)/R2(1);
60.         S2(j+1:j+2) = S2(j+1:j+2) - S2(j)*R2(2:3)/R2(1);
61.     end
62.
63. // Граничные условия
64. B(1) = S1(1);
65. C(1) = S1(2);
66. D(1) = S1(6);
67.
68. A(Nx+1) = S2(4);
69. B(Nx+1) = S2(5);
70. D(Nx+1) = S2(6);
71.   /*  C(1) = 48 - 36*A(2)/C(2)+B(2)*(18/(C(2))+36*B(2)/(C(2)*C(2)));
72.     B(1) = -25+A(2)*(18/C(2)+36*B(2)/(C(2)*C(2)))
73.     D(1) = 12*t(n+1)*t(n+1)*h + 3*D(6)/A(2)-36*D(3)/C(2)+D(2)*(18/C(2)+36*B(2)/(C(2)*C(2)));
74.     B(Nx+1) = 3  - C(Nx)/A(Nx);
75.     A(Nx+1) = -4 - B(Nx)/A(Nx);
76.     D(Nx+1) = 2*t(n+1)*t(n+1)*h*exp(b) - D(Nx)/A(Nx);*/
77.     W(n+1,:) = TrMxAlg(A,B,C,D)';
78.     end
79.     y = W;
80. endfunction
81. //W-Приближенное рещение
82. //U-Точное
83.     a = 0;b = 2;
84.     T = 1;
85.     Nx = 100; Nt = 100;
86.      h = (b-a)/Nx;
87.
88.
89.     tau = T/Nt;
90.     x = a:h:b;
91.     t = 0:tau:T;
92.     sigma = 0.5-h*h/(12*tau);
93.
94.     W = solve(Nx,Nt);
95. for i=1:Nx+1
96.     for n=1:Nt+1
97.         U(n,i) = u(x(i),t(n));
98.         end
99. end
100.
101. Epsilon=max(abs( W - U ));
102.
103.
104. Pucture1=scf(1);
105.
106. xtitle ('x', 't', 'Точное решение');
107.
108.
109. plot3d1(t,x,U);
110.
111. Picture2=scf(2);
112.
113. xtitle ('x', 't', 'Приближенное');
114.
115.
116. plot3d1(t,x,W);
117.
118. Picture3=scf(3);
119.
120. xtitle ('x', 't', 'Погрешность');
121.
122. plot3d1(t,x,abs(W - U));
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.

Top