急求一维抛物型方程的差分方法,要有数值解,精确解和误差分析,急求!
1条回答 默认 最新
- 技术专家团-Joel 2021-09-22 18:54关注
你好!比如一维热传导方程就是抛物线型的
% Number of Control Volumes N = 5; % Domain length L = 0.02; %[m] % Grid Spacing h = L/(N); % Diffusivity (Thermal conductivity) k = 0.5; %[W/m-K] % Uniform Heat Generation, A Source Term q = 1000e3; % [W/m^(3)] % The center of the first cell is at dx/2 & the last cell is at L-dx/2. x = h/2 : h : L-(h/2); % Cross sectional area of the 1D domain A = 1; %[m^2] % Left Boundary Condition T_a = 100; %[\circ c] % Right Boundary Condition T_b = 200; %[\circ c] % Initialising Temperature Nh = 2*N + 1; % the number of node with temperature T=zeros(Nh,1); % Left Boundary Condition T(1)= T_a; %[\circ c] % Right Boundary Condition T(Nh) = T_b; %[\circ c] % Interior of Domain index_x = 2:Nh-1; %% Formation of Coefficient Matrices, Numerical Solution % Tolerence eps = 1.e-6; % Iteration It = 0; T_old = T; Error = 2*eps; % Calculation while (Error > eps) It = It +1; T(index_x)= 0.5*(T(index_x-1) + T(index_x+1)+ (q*h*h/4)/(k)); % here divided by 4, cause h/2*h/2 Error = max(abs(T(:)-T_old(:))); if any(isnan(T(:))) || any(isinf(T(:))) fprintf(' Iteration Diverge \n'); return; end T_old = T; end T_Numeric = T(2:2:end); %% Analytical Results X_G = x; T_Analytical=((T_b-T_a)/L + (q/(2*k))*( L - X_G )).*(X_G) + T_a;%解析解 %% Error Analysis M_Error(N)=0; for i=1:N % This x is only use to compare results on Versteeg book k=find(abs(X_G - x(i))<=1e-3); % No need to find k (You can but no Need) % Absolute Error M_Error(i) = abs(T_Numeric(i)-T_Analytical(i)); % Relative Error Relative_Error(i) = abs(T_Numeric(i)-T_Analytical(i))/abs(T_Numeric(i)); % Percentage Error Percentage_Error(i) = Relative_Error(i).*100; end %% Post Processing % Data Table Comparison of Different Results Nodes = 1 : N; fprintf('Node\tGrid_{Loc.}\t\tT_{Num.}\t\tT_{An.}\t\tError\n\n'); fprintf('===\t\t=======\t\t=========\t\t=========\t\t=======\n'); for i = 1: N; fprintf('%3.0f\t\t%3.5f\t\t%3.5f\t\t%3.5f\t\t%3.5f\n',[Nodes(i); X_G(i); T_Numeric(i); T_Analytical(i); Percentage_Error(i)]); end % Graphical Representation plot(100*x,T_Numeric,'-o','MarkerSize',5,'MarkerFaceColor','b'); % hold on % plot([0, 100.*L], [T_a, T_b], 'd', 'LineWidth',3,'MarkerFaceColor','b'); hold on plot(100*X_G,T_Analytical,'r--','LineWidth',2); % hold on % plot([0, 100.*L], [T_a, T_b], 'd', 'LineWidth',3,'MarkerFaceColor','r'); % xlim([0 2]); % ylim([50 300]); %Boundary Points Inclusion xlabel('Distance (m)'); ylabel('Temperature [\circ C]'); legend('T_{Numeric}','T_{Analytical}'); title ('Numerical & Analytical Results Comparison'); grid on; %通过不断迭代误差控制在1e-6内 % Number of Control Volumes N = 5; % Domain length L = 0.02; %[m] % Grid Spacing h = L/(N); % Diffusivity (Thermal conductivity) k = 0.5; %[W/m-K] % Uniform Heat Generation, A Source Term q = 1000e3; % [W/m^(3)] % The center of the first cell is at dx/2 & the last cell is at L-dx/2. x = h/2 : h : L-(h/2); % Cross sectional area of the 1D domain A = 1; %[m^2] % Left Boundary Condition T_a = 100; %[\circ c] % Right Boundary Condition T_b = 200; %[\circ c] % Initialising Temperature Nh = 2*N + 1; % the number of node with temperature T=zeros(Nh,1); % Left Boundary Condition T(1)= T_a; %[\circ c] % Right Boundary Condition T(Nh) = T_b; %[\circ c] % Interior of Domain index_x = 2:Nh-1; %% Formation of Coefficient Matrices, Numerical Solution % Tolerence eps = 1.e-6; % Iteration It = 0; T_old = T; Error = 2*eps; % Calculation while (Error > eps) It = It +1; T(index_x)= 0.5*(T(index_x-1) + T(index_x+1)+ (q*h*h/4)/(k)); % here divided by 4, cause h/2*h/2 Error = max(abs(T(:)-T_old(:))); if any(isnan(T(:))) || any(isinf(T(:))) fprintf(' Iteration Diverge \n'); return; end T_old = T; end T_Numeric = T(2:2:end); %% Analytical Results X_G = x; T_Analytical=((T_b-T_a)/L + (q/(2*k))*( L - X_G )).*(X_G) + T_a;%解析解 %% 误差分析 M_Error(N)=0; for i=1:N % This x is only use to compare results on Versteeg book k=find(abs(X_G - x(i))<=1e-3); % No need to find k (You can but no Need) % Absolute Error M_Error(i) = abs(T_Numeric(i)-T_Analytical(i)); % Relative Error Relative_Error(i) = abs(T_Numeric(i)-T_Analytical(i))/abs(T_Numeric(i)); % Percentage Error Percentage_Error(i) = Relative_Error(i).*100; end %% Post Processing % Data Table Comparison of Different Results Nodes = 1 : N; fprintf('Node\tGrid_{Loc.}\t\tT_{Num.}\t\tT_{An.}\t\tError\n\n'); fprintf('===\t\t=======\t\t=========\t\t=========\t\t=======\n'); for i = 1: N; fprintf('%3.0f\t\t%3.5f\t\t%3.5f\t\t%3.5f\t\t%3.5f\n',[Nodes(i); X_G(i); T_Numeric(i); T_Analytical(i); Percentage_Error(i)]); end % Graphical Representation plot(100*x,T_Numeric,'-o','MarkerSize',5,'MarkerFaceColor','b'); % hold on % plot([0, 100.*L], [T_a, T_b], 'd', 'LineWidth',3,'MarkerFaceColor','b'); hold on plot(100*X_G,T_Analytical,'r--','LineWidth',2); % hold on % plot([0, 100.*L], [T_a, T_b], 'd', 'LineWidth',3,'MarkerFaceColor','r'); % xlim([0 2]); % ylim([50 300]); %Boundary Points Inclusion xlabel('Distance (m)'); ylabel('Temperature [\circ C]'); legend('T_{Numeric}','T_{Analytical}'); title ('Numerical & Analytical Results Comparison'); grid on;
误差分析
Node Grid_{Loc.} T_{Num.} T_{An.} Error === ======= ========= ========= ======= 1 0.00200 146.00000 146.00000 0.00000 2 0.00600 213.99999 214.00000 0.00000 3 0.01000 249.99999 250.00000 0.00000 4 0.01400 253.99999 254.00000 0.00000 5 0.01800 226.00000 226.00000 0.00000
本回答被题主选为最佳回答 , 对您是否有帮助呢?解决 无用评论 打赏 举报
悬赏问题
- ¥15 fpga自动售货机数码管(相关搜索:数字时钟)
- ¥20 Python安装cvxpy库出问题
- ¥15 用前端向数据库插入数据,通过debug发现数据能走到后端,但是放行之后就会提示错误
- ¥15 python天天向上类似问题,但没有清零
- ¥30 3天&7天&&15天&销量如何统计同一行
- ¥30 帮我写一段可以读取LD2450数据并计算距离的Arduino代码
- ¥15 C#调用python代码(python带有库)
- ¥15 活动选择题。最多可以参加几个项目?
- ¥15 飞机曲面部件如机翼,壁板等具体的孔位模型
- ¥15 vs2019中数据导出问题