Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- clear();
- a = 0;
- b = 2;
- T = 1;
- Nx = 200;
- Nt = 400;
- sigma = 0.5;
- function y = u_exact(x,t)
- y = t*t*exp(x);
- endfunction
- function y=f(x,t)
- y= (2*t-t*t)*exp(x);
- endfunction
- function y=err(Nx,Nt,u_ex)
- y = abs( solve(Nx,Nt) - u_ex );
- endfunction
- function y=progonka(A,B,C,D)
- // A,B,C,D - массивы поддиагональных эл-в, диагональных, наддиагональных и правой части соответственно
- N = length(B);
- c(1) = C(1)/B(1);
- d(1) = D(1)/B(1);
- for i = 2:N-1
- k = B(i) - A(i)*c(i-1);
- c(i) = C(i)/k;
- d(i) = (D(i) - A(i)*d(i-1))/k;
- end
- y(N) = (D(N) - A(N)*d(N-1))/(B(N) - A(N)*c(N-1));
- for i = N-1:-1:1
- y(i) = d(i) - c(i)*y(i+1);
- end
- endfunction
- //u_0 = 0 в нашем случае
- function y=solve(Nx,Nt)
- h = (b-a)/Nx;
- tau = T/Nt;
- x = a:h:b;
- t = 0:tau:T;
- U = zeros(Nt+1,Nx+1);
- A = zeros(Nx+1,1);
- B = zeros(Nx+1,1);
- C = zeros(Nx+1,1);
- D = zeros(Nx+1,1);
- //цикл по t
- for n = 1:Nt
- //Заполняем диагонали
- for i = 2:Nx
- B(i) = -h*h-2*tau*sigma;
- A(i) = tau*sigma;
- C(i) = tau*sigma;
- 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));
- end
- // Граничные условия
- C(1) = 4 + B(2)/C(2);
- B(1) = -3 + A(2)/C(2);
- D(1) = 2*t(n+1)*t(n+1)*h + D(2)/C(2);
- B(Nx+1) = 3 - C(Nx)/A(Nx);
- A(Nx+1) = -4 - B(Nx)/A(Nx);
- D(Nx+1) = 2*t(n+1)*t(n+1)*h*exp(b) - D(Nx)/A(Nx);
- U(n+1,:) = progonka(A,B,C,D)';
- end
- y = U;
- endfunction
- U = solve(Nx,Nt);
- h = (b-a)/Nx;
- tau = T/Nt;
- x = a:h:b;
- t = 0:tau:T;
- for i=1:Nx+1
- for n=1:Nt+1
- U_ex(n,i) = u_exact(x(i),t(n));
- end
- end
- f1=scf(1); //creates figure with id==4 and make it the current one
- f1.color_map = jetcolormap(64);
- xtitle ('x', 't', 'exact_solution');
- plot3d1(t,x,U_ex);
- f2=scf(2); //creates figure with id==4 and make it the current one
- f2.color_map = jetcolormap(64);
- xtitle ('x', 't', 'U');
- plot3d1(t,x,U);
- f3=scf(3); //creates figure with id==4 and make it the current one
- f3.color_map = jetcolormap(64);
- xtitle ('x', 't', 'error');
- mis=max(err(Nx,Nt,U_ex));
- plot3d1(t,x,err(Nx,Nt,U_ex));
- /*
- // pogresnosti
- NxL = [50, 100,200,400,800];
- NtL = [25,100,400,1600,6400];
- errL = zeros(5);
- for j = 1:5
- Nx = NxL(j);
- Nt = NtL(j);
- h = (b-a)/Nx;
- tau = T/Nt;
- x = a:h:b;
- t = 0:tau:T;
- U_ex = zeros(Nt+1,Nx+1);
- for i=1:Nx+1
- for n=1:Nt+1
- U_ex(n,i) = u_exact(x(i),t(n));
- end
- end
- errL(j) = max(err(Nx,Nt,U_ex));
- end
- */
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement