　　這篇文章給出（1）Jacobi與SOR迭代法的實作與性能比較及（2）均勻間距與Chebyshev插值的實作、性能分析及二者生成的插值誤差比較，給出完整的實作代碼，沒有進行性能優化，僅供參考，

（1）Jacobi與SOR迭代法的實作與性能比較

一、舉例計算

給出線性方程組：

　　其中n=100或者n=1000（任選一種，在本報告測驗中，選取了n=100），使用Jacobi迭代法和SOR迭代法（=1,1.25,1.5）解此方程，計算結果精確到小數點后8位，結果輸出小數點后至少12位，報告所需要的步數和誤差（用無窮范數表示），比較計算結果，

二、Jacobi迭代法的實作

2.1、實作環境

MATLAB2018b

2.2、實作原理

　　整個Jacobi實作分下面幾步實作：

　　首先需要給定這個系統的輸入：矩陣A、向量b及初始向量x0，根據題目A與b是確定的，選定n=100，又矩陣A的是三對角矩陣，實作方式用到了MATLAB的內置函式diag來快速實作，具體實作如下：

　　A = diag(2*ones(100,1))+diag(ones(99,1),1)+diag(ones(99,1),-1);

而向量b直接取值即可，初始向量x0給定為：x0= 0.01*ones(100,1);

　　在確定好輸入之后，其二需要給出Jacobi的運算方式，而其關鍵是計算向量的各個元素值，依據Jacobi的的計算公式：

來進行運算，撰寫的函式在規定的迭代次數之前進行運算；

　　如果在規定的迭代次數之前所計算的誤差精度大于規定的誤差值時，即或者，繼續執行下面的步驟；

　　在上一步結束之后?需要更新的值，即的值賦給，更新迭代的次數，之后再轉到第二步繼續回圈運行，如此反復運行，直到結束，

2.3、資料測驗

　　在n=100時，在MATLAB輸入如下命令：function [y] = myJocobi(u_x0,u_A,u_b,S_error,N_times);而為了探討經過Jacobi迭代運算之后誤差的計算程序用二范數和無窮范數區別，特地做了實驗來進行說明，當采用二范數時，計算結果如如圖1所示：

圖1 Jacobi迭代結果（誤差精度計算采用二范數)

　　當誤差精度計算采用無窮范數時，計算結果如如圖2所示：

圖2 Jacobi迭代結果（誤差精度計算采用無窮范數)

2.4、資料分析

　　當n=100時，在實作精度等要求的前提下，誤差精度計算采用二范數和無窮范數的Jacobi迭代的次數均達到3萬多次，分別為34541及30536次，可以預測當n值取得很大的時候，Jacobi迭代效率會很低，

三、SOR迭代法的實作

3.1、實作環境

MATLAB2018b

3.2、實作原理

　　整個SOR迭代演算法實作分下面幾步實作：

　　首先需要給定這個系統的輸入：矩陣A、向量b及初始向量x0，根據題目A與b是確定的，選定n=100，又矩陣A的是三對角矩陣，而相對于Jacobi迭代，SOR迭代多了一個超松弛因子w（w=1,1.25,1.5），三對角矩陣A實作方式用到了MATLAB的內置函式diag來快速實作，具體實作如下：

A = diag(2*ones(100,1))+diag(ones(99,1),1)+diag(ones(99,1),-1);

而向量b直接依題取值即可，初始向量x0給定為：x0= 0.01*ones(100,1);

　　在確定好輸入之后，其二需要給出SOR迭代的運算函式，而其關鍵是計算向量x的各個元素值x_i，依據SOR的的計算公式：

來進行運算，撰寫的函式在規定的迭代次數之前進行運算；

　　如果在規定的迭代次數之前所計算的誤差精度大于規定的誤差值時，即,繼續執行下面步驟；

　　在上一步結束之后?需要更新的值，即的值賦給，更新迭代的次數，之后再轉到第二步繼續回圈運行，如此反復運行，直到結束，

3.3、資料測驗

　　在n=100時，在MATLAB輸入如下命令：

　　[x_result,final_error] = mySOR(x0,A,b,1e-8,50000,1);

　　而為了探討經過SOR迭代運算程序中超松弛因子w值不同對迭代的影響，現進行了w=1,1.25,1.5三次測驗，測驗結果如下：

n=100，超松弛因子=1時，SOR迭代結果如圖3所示：

圖3 超松弛因子=1時，SOR迭代結果

n=100，超松弛因子=1.25時，SOR迭代結果如圖4所示：

圖4 超松弛因子=1.25時，SOR迭代結果

n=100，超松弛因子=1.5時，SOR迭代結果如圖5所示：

圖5 超松弛因子=1.5時，SOR迭代結果

3.4、資料分析

　　當n=100時，在實作精度等要求的前提下，誤差精度計算采用無窮范數的SOR迭代的次數在超松弛因子w分別取1、1.25、1.5時分別為12115、7584及4411次，可以看到，SOR迭代的整體速度較快，而在選取不同的超松弛因子值時，帶來的最終結果卻相差較大，因此可以預測當n值取得很大的時候，SOR迭代效率會比較高，且在選取合適的超松弛因子w時，效率提高會有數量級的提高，

四、Jacobi與SOR迭代法的運算結果比較

　　在前提條件相同的基礎上，n=100時，Jacobi的迭代次數為34541及30536次，而SOR迭代的次數為12115、7584及4411次，誤差精度計算采用無窮范數時，最大的次數比為30536/4411，比例約為7點多，可以看到Jacobi相比與SOR迭代在效率上差很多，而且隨著n的增大，兩者之間的效率會相差越來越明顯，

（2）均勻間距與Chebyshev插值的實作、性能分析及二者生成的插值誤差比較

一、舉例計算

　　令f(x)=exp(|x|)，通過同時在區間[-1,1]上畫出兩種型別的n階多項式，其中n(N)=10,20，比較均勻間距的插值和Chebyshev插值，對于均勻間距的插值，左右插值的基點為-1和1，以0.01的步長采樣，對每種插值型別生成插值誤差，并畫出來進行比較，在這個問題里看是否能看到Runge現象，

二、均勻間距插值的實作

2.1、實作環境

MATLAB2018b

2.2、實作原理

　　在實驗中，采用的是Lagrange來實作均勻間距插值和Chebyshev插值，實驗中首先需要實作的就是Lagrange演算法插值，對于n次的Lagrange插值，構造的基函式為：

即得到滿足插值條件：

得到插值函式為：

　　為了實作n=10,20的均勻間距插值，即對區間[-1,1]進行劃相應的等分，核心實作為Lagrange函式，而Lagrange函式介面輸入引數有三個，其一為[-1,1]區間中0.01步長得到的u_x，其二為經過均勻間距的變數值x_hk，其三為將劃分的變數值經過函式y=exp(|x|)得到的輸出值y_hk(hk =11對應n=10或者12對應n=20)，在Lagrange函式中，通過使用上文已經構造好的基函式及插值函式公式來進行計算，最后得到經過均勻間距插值的資料結果，

2.3、資料測驗

　　在資料進行測驗時，在MATLAB分別輸入如下命令：

　　Ln1 = myLagrange(x11,y11,u_x);

　　Ln2 = myLagrange(x12,y12,u_x);

　　其中，x11與y11均為均勻間隔插值N(n)=10時，x12與y12均為均勻間隔插值N(n)=20時對資料的處理，最后得到結果如圖6和圖7所示，其中圖6為將原函式與均勻間隔插值N(n)=10及均勻間隔插值N(n)=20影像分開比較，圖7為三者一同進行比較：

圖6 均勻間隔插值與原函式進行比較-1

圖7 均勻間隔插值與原函式進行比較-2

2.4、資料分析

　　從測驗結果可以看到，使用均勻間隔插值的時候時可以逼近原函式的，但是存在當n太小的時候，插值函式整體上與原函式有一定的區別，但是在n增大的情況下，插值函式中間段部分與原函式逼近，但是兩端出現了震蕩，即出現了Runge現象，

三、Chebyshev插值的實作

3.1、實作環境

MATLAB2018b

3.2、實作原理

　　對于Chebyshev插值，其插值區間必須為[-1,1]，而題目是剛好吻合的，所以能夠直接使用Chebyshev插值，而對于Chebyshev插值，選取的插值節點以如下公式進行計算：

而在MATLAB中，向量的索引值從1開始，所以需要對公式進行變形，即：

　　為了實作n=10,20的Chebyshev插值，即通過上面所列寫的公式取節點值，同樣Lagrange函式介面輸入引數有三個，其一為[-1,1]區間中0.01步長得到的，其二為通過上面所列寫的公式所取節點值x_hk，其三為將所取得節點值經過函式y=exp(|x|)得到的輸出值y_hk(hk =21對應n=10或者22對應n=20)，在Lagrange函式中，通過使用上文已經構造好的基函式及插值函式公式來進行計算，最后得到經過Chebyshev插值的資料結果，

3.3、資料測驗

　　在資料進行測驗時，在MATLAB分別輸入如下命令：

　　Ln3 = myLagrange(x21,y21,u_x);

　　Ln4 = myLagrange(x22,y22,u_x);

　　其中，x21與y21均為Chebyshev插值N(n)=10時，x22與y22均為Chebyshev插值N(n)=20時對資料的處理，最后得到結果如圖8和圖9所示，其中圖8為將原函式與Chebyshev插值N(n)=10及Chebyshev插值N(n)=20影像分開比較，圖9為三者一同進行比較：

圖8 Chebyshev插值與原函式進行比較-1

圖9 Chebyshev插值與原函式進行比較-2

3.4、資料分析

　　從測驗結果可以看到，使用Chebyshev插值的時候時可以更快地逼近原函式的，且當n比較小的時候，插值函式整體上與原函式區別不大，隨著n的增大，插值函式與原函式越來越逼近，在誤差允許情況下，能夠很好的替換原函式，且避免了均勻間隔插值兩端出現震蕩的情況（Runge現象），

四、均勻間距插值與Chebyshev插值誤差比較

　　為了進行均勻間隔插值與Chebyshev插值誤差比較，需要計算誤差余項，通過對誤差余項的定義得到余項的運算式：

　　通過對誤差余項的化簡可以得到：

而在本次實驗中f(x)=exp(|x|)，通過對化簡，得到余項公式：

　　均勻間隔插值與Chebyshev插值通過計算化簡后的余項即進行比較，其中在N(n)=10時，輸入Rn_1 = myRn_solve(u_x,x11);

Rn_3 = myRn_solve(u_x,x21);

得到當前N(n)值時,間隔與Chebyshev插值誤差，結果如圖10所示：

圖10 n=10時，插值誤差（余項）的對比

　　在N(n)=20時，輸入Rn_2 = myRn_solve(u_x,x12);

　　　　　　　　　 Rn_4 = myRn_solve(u_x,x22);

得到當前N(n)值時,均勻間隔與Chebyshev插值誤差，結果如圖11所示：

圖11 n=20時，插值誤差（余項）的對比

　　通過測驗結果可以看出，在同一基礎條件下，均勻間隔插值與Chebyshev插值的誤差主要是在兩端點處相差較大，而且隨著n的增大，端點出的相差程度會更明顯，從這也反映了均勻間隔插值帶來的Runge現象影響，

附錄

（1）Jacobi與SOR迭代法的實作代碼

math_test.m

 1 %Jacobi迭代計算
 2 %function [y] = myJocobi(u_x0,u_A,u_b,S_error,N_times) 
 3 % [x_result,final_error] = myJacobi(x0,A,b,1e-8,50000);
 4 % 
 5 % if abs(final_error) < 1e-8
 6 %     fprintf('final_error*10^8=%f\n',abs(final_error)*1e8);
 7 % end
 8 
 9 
10 
11 
12 
13 %SOR迭代計算
14 %function [xx,y] = mySOR(u_x0,u_A,u_b,S_error,N_times,u_w)
15 [x_result,final_error] = mySOR(x0,A,b,1e-8,50000,1);
16 
17 if abs(final_error) < 1e-8
18     fprintf('final_error*10^8=%f\n',abs(final_error)*1e8);
19 end

myJacobi.m

 1 function [xx,y] = myJacobi(u_x0,u_A,u_b,S_error,N_times)
 2 %S_error 為規定的需要達到誤差值
 3 %N_times 為規定的最大迭代次數
 4     u_count = 1; 
 5     u_error = 1; %初始化默認誤差為1
 6     
 7     u_x = zeros(100,1); %初始化
 8     u_ax_toal = 0;
 9     
10     while(abs(u_error) > S_error)
11         
12        
13         
14         for i=1:100
15             
16             for j=1:100
17                 if j ~= i  %如果 cnt不等于i
18                     u_ax = u_A(i,j)*u_x0(j);
19                     u_ax_toal = u_ax_toal + u_ax;
20                 end
21             end
22             
23             u_x(i) = (u_b(i) - u_ax_toal)/u_A(i,i);
24             u_ax_toal = 0;
25         end
26         
27         
28         %u_error = norm(u_x-u_x0,2);   %求二范數
29         u_error = norm(u_x-u_x0,inf);   %求無窮范數
30         
31         u_x0 = u_x;
32         
33         if u_count>=N_times
34             break;
35         end
36         
37         fprintf('u_error  %d: %16.14f\n',u_count,u_error);
38         
39         u_count = u_count + 1;
40     end
41     
42     y = u_error;
43     xx = u_x0;
44 end

mySOR.m

 1 function [xx,y] = mySOR(u_x0,u_A,u_b,S_error,N_times,u_w)
 2    
 3 %S_error 為規定的需要達到誤差值
 4 %N_times 為規定的最大迭代次數
 5     u_count = 1; 
 6     u_error = 1; %初始化默認誤差為1
 7     
 8     u_x = zeros(100,1); %初始化
 9     u_ax_toal_last = 0;
10     u_ax_toal_now = 0;
11      
12     
13     while(abs(u_error) > S_error)
14         
15         for i=1:100
16             if i==1
17                 for j=2:100
18 
19                         u_ax = u_A(i,j)*u_x0(j);
20                         u_ax_toal_last = u_ax_toal_last + u_ax;
21                 end
22                 
23                 u_x(i) = (1-u_w)*u_x0(i) + u_w*((u_b(i) - u_ax_toal_last)/u_A(i,i));
24                 
25                 u_ax_toal_last = 0;
26             end
27             
28             if i>=2 && i<=99
29                 for j=1:i-1
30 
31                         u_ax = u_A(i,j)*u_x(j);
32                         u_ax_toal_now = u_ax_toal_now + u_ax;
33                 end
34 
35                 for j=i+1:100
36 
37                         u_ax = u_A(i,j)*u_x0(j);
38                         u_ax_toal_last = u_ax_toal_last + u_ax;
39                 end
40                 
41                 u_x(i) = (1-u_w)*u_x0(i) + u_w*((u_b(i) - u_ax_toal_last - u_ax_toal_now)/u_A(i,i));
42                 
43                 u_ax_toal_last = 0;
44                 u_ax_toal_now = 0;
45             end
46             
47             
48             if i==100
49                 for j=1:99
50 
51                         u_ax = u_A(i,j)*u_x(j);
52                         u_ax_toal_now = u_ax_toal_now + u_ax;
53                 end
54                 
55                 u_x(100) = (1-u_w)*u_x0(100) + u_w*((u_b(100) - u_ax_toal_now)/u_A(100,100));
56                 
57                 u_ax_toal_now = 0;
58             end
59             
60             
61         end
62         
63         
64         %u_error = norm(u_x-u_x0,2);   %求二范數
65         u_error = norm(u_x-u_x0,inf);   %求無窮范數
66         
67         u_x0 = u_x;
68         
69         if u_count>=N_times
70             break;
71         end
72         
73         fprintf('u_error  %d: %16.14f\n',u_count,u_error);
74         
75         u_count = u_count + 1;
76     end
77     
78     y = u_error;
79     xx = u_x0; 
80 end

three_matrix.m

1 A = diag(2*ones(100,1))+diag(ones(99,1),1)+diag(ones(99,1),-1);
2 
3 x0 = 0.01*ones(100,1);
4 
5 
6 b = zeros(100,0);
7 b(1) = 1;
8 b(100) = -1;
9 b = b';

（2）均勻間距與Chebyshev插值的實作

math_test_plot_Ln.m

 1 u_x = -1:0.01:1;
 2 u_y = exp(1).^abs(u_x);  % y = e^|x|
 3 
 4 
 5 
 6 % %均勻間距插值，使用等間距節點作為插值節點
 7 % % title('均勻間距插值')
 8 % 
 9 % %subplot(3,1,1);
10 % plot(u_x,u_y,'-.b')
11 % hold on
12 % %title('原函式')
13 % title('均勻間距插值')
14 % 
15 % 
16 % %subplot(3,1,2);
17 % x11 = -1:0.2:1; %n=10
18 % y11 = exp(1).^abs(x11);
19 % Ln1 = myLagrange(x11,y11,u_x);
20 % % Rn_1 = exp(1)*myRn_solve(u_x,x11);
21 % Rn_1 = myRn_solve(u_x,x11);
22 % plot(u_x,Ln1,'-k')
23 % % title('均勻間距插值 N = 10')
24 % hold on
25 % 
26 % %subplot(3,1,3);
27 % x12 = -1:0.1:1; %n=20
28 % y12 = exp(1).^abs(x12);
29 % Ln2 = myLagrange(x12,y12,u_x);
30 % Rn_2 = myRn_solve(u_x,x12);
31 % plot(u_x,Ln2,'-r')
32 % % title('均勻間距插值 N = 20')
33 % hold on
34 % 
35 % legend('原函式','均勻間距插值 N = 10','均勻間距插值 N = 20')
36 
37 
38 
39 
40 % %切比雪夫插值，使用切比雪夫節點作為插值節點
41 % subplot(3,1,1);
42 % 
43 plot(u_x,u_y,'-.b')
44 % title(' 原函式 ')
45 hold on
46 title('Chebyshev插值 N=10');
47 
48 % subplot(3,1,2);
49 x21 = myChebyshev(10);%n=10
50 y21 = exp(1).^abs(x21);
51 Ln3 = myLagrange(x21,y21,u_x);
52 Rn_3 = myRn_solve(u_x,x21);
53 plot(u_x,Ln3,'-r');
54 % title('Chebyshev插值 N=10');
55 hold on
56 
57 % subplot(3,1,3);
58 x22 = myChebyshev(20);%n=20
59 y22 = exp(1).^abs(x22);
60 Ln4 = myLagrange(x22,y22,u_x);
61 Rn_4 = myRn_solve(u_x,x22);
62 plot(u_x,Ln4,'-b');
63 % title('Chebyshev插值 N=20');
64 
65 legend('原函式','Chebyshev插值 N=10','Chebyshev插值 N=20')

math_test_plot_Rn.m

 1 u_x = -1:0.01:1;
 2 u_y = exp(1).^abs(u_x);  % y = e^|x|
 3 
 4 
 5 
 6 % 均勻間距插值
 7 
 8 x11 = -1:0.2:1; %n=10
 9 y11 = exp(1).^abs(x11);
10 Ln1 = myLagrange(x11,y11,u_x);
11 % Rn_1 = exp(1)*myRn_solve(u_x,x11);
12 Rn_1 = myRn_solve(u_x,x11);
13 
14 
15 
16 x12 = -1:0.1:1; %n=20
17 y12 = exp(1).^abs(x12);
18 Ln2 = myLagrange(x12,y12,u_x);
19 Rn_2 = myRn_solve(u_x,x12);
20 
21 
22 
23 
24 
25 % 切比雪夫插值，使用切比雪夫節點作為插值節點
26 
27 x21 = myChebyshev(10);%n=10
28 y21 = exp(1).^abs(x21);
29 Ln3 = myLagrange(x21,y21,u_x);
30 Rn_3 = myRn_solve(u_x,x21);
31 
32 
33 
34 x22 = myChebyshev(20);%n=20
35 y22 = exp(1).^abs(x22);
36 Ln4 = myLagrange(x22,y22,u_x);
37 Rn_4 = myRn_solve(u_x,x22);
38 
39 
40 
41 %繪制插值余項的最大值
42 plot(u_x,Rn_1,'-k',u_x,Rn_3,'-r');
43 title('N=10時 插值誤差(余項)的對比');
44 legend('均勻間距插值','Chebyshev插值');
45 
46 % %繪制插值余項的最大值
47 % plot(u_x,Rn_2,'-k',u_x,Rn_4,'-r');
48 % title('N=20時 插值誤差(余項)的對比');
49 % legend('均勻間距插值','Chebyshev插值');

myChebyshev.m

1 function x = myChebyshev(n)
2    
3 x = zeros(1,n+1);
4 for k = 1:n+1
5     x(k) = cos(((2*k-1)*pi)/(2*(n+1)));  %此處從1開始，所以是2*k-1,否則為2*k+1
6 end

myLagrange.m

 1 function y = myLagrange(x0,y0,x)
 2 %f11 = myLagrange(x11,y11,u_x);
 3 % x11 = -1:0.2:1; %n=10
 4 % y11 = exp(1).^abs(x11);
 5 
 6 
 7 u_N=1:length(x0); 
 8 
 9 
10 
11 y=zeros(size(x));
12 
13 for i=u_N
14     u_ij=find(u_N~=i);   
15     y1=1;
16     for j=1:length(u_ij)
17         y1=y1.*(x-x0(u_ij(j)));  
18     end
19     
20     y2=1;
21     for j=1:length(u_ij)
22         y2=y2.*(x0(i)-x0(u_ij(j)));  
23     end
24     
25     y=y+y0(i)*y1/y2;
26 
27 end

myRn_solve.m

 1 function u_Rn = myRn_solve(x,x0)
 2 
 3 % x11 = -1:0.2:1; %n=10
 4 % y11 = exp(1).^abs(x11);
 5 % f11 = myLagrange(x11,y11,u_x);
 6 % Rn_1 = exp(1)*myRn_solve(u_x,x11);
 7 
 8 u_Wn = zeros(1,length(x));
 9 u_Rn = zeros(1,length(x));
10 for i = 1:length(x)
11     u_Wn(i) = prod(x(i)*ones(1,length(x0))-x0);
12 end
13 
14 u_Rn = u_Wn*exp(1)/factorial(length(x0));
15 end

轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/256221.html

標籤：其他

上一篇：不動點及牛頓法迭代工程實作

下一篇：我憑借這份999頁Java面試pdf！拿下了美團、螞蟻金服、騰訊、位元組跳動offer

Jacobi與SOR迭代法的實作與性能比較及均勻間距與Chebyshev插值的實作、性能分析及二者生成的插值誤差比較

（1）Jacobi與SOR迭代法的實作與性能比較

一、舉例計算

二、Jacobi迭代法的實作

2.1、實作環境

2.2、實作原理

2.3、資料測驗

2.4、資料分析

三、SOR迭代法的實作

3.1、實作環境

3.2、實作原理

3.3、資料測驗

3.4、資料分析

四、Jacobi與SOR迭代法的運算結果比較

（2）均勻間距與Chebyshev插值的實作、性能分析及二者生成的插值誤差比較

一、舉例計算

二、均勻間距插值的實作

2.1、實作環境

2.2、實作原理

2.3、資料測驗

2.4、資料分析

三、Chebyshev插值的實作

3.1、實作環境

3.2、實作原理

3.3、資料測驗

3.4、資料分析

四、均勻間距插值與Chebyshev插值誤差比較

附錄

（1）Jacobi與SOR迭代法的實作代碼

（2）均勻間距與Chebyshev插值的實作