Heat Equation程式-(本程式使用向后差分法).pdfVIP

Heat Equation程式 (本程式使用向後差分法) 授課教師: 李義明 教授 助 教: 余紹銘、郭益廷 日 期: 2007/12/7 科學計算(一) 1 國立交通大學科管所 2007年12月7 日 Heat Equation程式使用說明 李義明老師 function [u,x,t] = my_heat(t_0,t_f,x_0,x_f,M,N,type,T1,T2) % define the mesh in space dx = (x_f-x_0)/N; x = x_0:dx:x_f; x = x; % define the mesh in time dt = (t_f-t_0)/M; t = t_0:dt:t_f; % define the ratio r r = dt/dx^2; fprintf(r = %g\n,r); if r 0.5 warning(r 0.5) dt = 0.499*dx^2; t = t_0:dt:t_f; fprintf(modified r = %g\n,dt/dx^2); end 科學計算(一) 2 國立交通大學科管所 2007年12月7 日 Heat Equation程式使用說明 李義明老師 M = length(t); N = length(x); u = zeros(N,M); if type == 0 % Dirichlet B.C u(1,:) = T1; u(N,:) = T2; for i = 2:M % Time runs A = zeros(N); A(1,1) = 1; u(1,i) = T1; A(N,N) = 1; u(N,i) = T2; for j = 2:N-1 A(j,j) = 2*r + 1; A(j,j+1) = -r; A(j,j-1) = -r; end u(:,i) = A\u(:,i-1); end 科學計算(一) 3 國立交通大學科管所 2007年12月7 日 Heat Equation程式使用說明 李義明老師 elseif type == 1 % Neumann at the left bounary, Dirichlet at the right boundary; u(N,:) = T2; for i = 2:M % Time runs A = zeros(N); A(1,1) = 2*r+1; A(1,2) = -2*r; A(N,N) = 1; u(N,i) = T2; for j = 2:N-1 A(j,j) = 2*r + 1; A(j,j+1) = -r; A(j,j-1) = -r; end u(:,i) = A\u(:,i-1); end 科學計算(一) 4 國立交通大學科管所 2007年12月7 日 Heat Equation程式使用說明