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

Top