• API
• FAQ
• Tools
• Archive
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.

Top