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

　　這篇文章給出（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

%Jacobi迭代計算
%function [y] = myJocobi(u_x0,u_A,u_b,S_error,N_times) 
% [x_result,final_error] = myJacobi(x0,A,b,1e-8,50000);
% 
% if abs(final_error) < 1e-8
%     fprintf('final_error*10^8=%f\n',abs(final_error)*1e8);
% end





%SOR迭代計算
%function [xx,y] = mySOR(u_x0,u_A,u_b,S_error,N_times,u_w)
[x_result,final_error] = mySOR(x0,A,b,1e-8,50000,1);

if abs(final_error) < 1e-8
    fprintf('final_error*10^8=%f\n',abs(final_error)*1e8);
end

 

myJacobi.m

function [xx,y] = myJacobi(u_x0,u_A,u_b,S_error,N_times)
%S_error 為規定的需要達到誤差值
%N_times 為規定的最大迭代次數
    u_count = 1; 
    u_error = 1; %初始化默認誤差為1
    
    u_x = zeros(100,1); %初始化
    u_ax_toal = 0;
    
    while(abs(u_error) > S_error)
        
       
        
        for i=1:100
            
            for j=1:100
                if j ~= i  %如果 cnt不等於i
                    u_ax = u_A(i,j)*u_x0(j);
                    u_ax_toal = u_ax_toal + u_ax;
                end
            end
            
            u_x(i) = (u_b(i) - u_ax_toal)/u_A(i,i);
            u_ax_toal = 0;
        end
        
        
        %u_error = norm(u_x-u_x0,2);   %求二范數
        u_error = norm(u_x-u_x0,inf);   %求無窮范數
        
        u_x0 = u_x;
        
        if u_count>=N_times
            break;
        end
        
        fprintf('u_error  %d: %16.14f\n',u_count,u_error);
        
        u_count = u_count + 1;
    end
    
    y = u_error;
    xx = u_x0;
end

 

mySOR.m

function [xx,y] = mySOR(u_x0,u_A,u_b,S_error,N_times,u_w)
   
%S_error 為規定的需要達到誤差值
%N_times 為規定的最大迭代次數
    u_count = 1; 
    u_error = 1; %初始化默認誤差為1
    
    u_x = zeros(100,1); %初始化
    u_ax_toal_last = 0;
    u_ax_toal_now = 0;
     
    
    while(abs(u_error) > S_error)
        
        for i=1:100
            if i==1
                for j=2:100

                        u_ax = u_A(i,j)*u_x0(j);
                        u_ax_toal_last = u_ax_toal_last + u_ax;
                end
                
                u_x(i) = (1-u_w)*u_x0(i) + u_w*((u_b(i) - u_ax_toal_last)/u_A(i,i));
                
                u_ax_toal_last = 0;
            end
            
            if i>=2 && i<=99
                for j=1:i-1

                        u_ax = u_A(i,j)*u_x(j);
                        u_ax_toal_now = u_ax_toal_now + u_ax;
                end

                for j=i+1:100

                        u_ax = u_A(i,j)*u_x0(j);
                        u_ax_toal_last = u_ax_toal_last + u_ax;
                end
                
                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));
                
                u_ax_toal_last = 0;
                u_ax_toal_now = 0;
            end
            
            
            if i==100
                for j=1:99

                        u_ax = u_A(i,j)*u_x(j);
                        u_ax_toal_now = u_ax_toal_now + u_ax;
                end
                
                u_x(100) = (1-u_w)*u_x0(100) + u_w*((u_b(100) - u_ax_toal_now)/u_A(100,100));
                
                u_ax_toal_now = 0;
            end
            
            
        end
        
        
        %u_error = norm(u_x-u_x0,2);   %求二范數
        u_error = norm(u_x-u_x0,inf);   %求無窮范數
        
        u_x0 = u_x;
        
        if u_count>=N_times
            break;
        end
        
        fprintf('u_error  %d: %16.14f\n',u_count,u_error);
        
        u_count = u_count + 1;
    end
    
    y = u_error;
    xx = u_x0; 
end

 

three_matrix.m

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

x0 = 0.01*ones(100,1);


b = zeros(100,0);
b(1) = 1;
b(100) = -1;
b = b';

 

 

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

math_test_plot_Ln.m

u_x = -1:0.01:1;
u_y = exp(1).^abs(u_x);  % y = e^|x|



% %均勻間距插值，使用等間距節點作為插值節點
% % title('均勻間距插值')
% 
% %subplot(3,1,1);
% plot(u_x,u_y,'-.b')
% hold on
% %title('原函數')
% title('均勻間距插值')
% 
% 
% %subplot(3,1,2);
% x11 = -1:0.2:1; %n=10
% y11 = exp(1).^abs(x11);
% Ln1 = myLagrange(x11,y11,u_x);
% % Rn_1 = exp(1)*myRn_solve(u_x,x11);
% Rn_1 = myRn_solve(u_x,x11);
% plot(u_x,Ln1,'-k')
% % title('均勻間距插值 N = 10')
% hold on
% 
% %subplot(3,1,3);
% x12 = -1:0.1:1; %n=20
% y12 = exp(1).^abs(x12);
% Ln2 = myLagrange(x12,y12,u_x);
% Rn_2 = myRn_solve(u_x,x12);
% plot(u_x,Ln2,'-r')
% % title('均勻間距插值 N = 20')
% hold on
% 
% legend('原函數','均勻間距插值 N = 10','均勻間距插值 N = 20')




% %切比雪夫插值，使用切比雪夫節點作為插值節點
% subplot(3,1,1);
% 
plot(u_x,u_y,'-.b')
% title(' 原函數 ')
hold on
title('Chebyshev插值 N=10');

% subplot(3,1,2);
x21 = myChebyshev(10);%n=10
y21 = exp(1).^abs(x21);
Ln3 = myLagrange(x21,y21,u_x);
Rn_3 = myRn_solve(u_x,x21);
plot(u_x,Ln3,'-r');
% title('Chebyshev插值 N=10');
hold on

% subplot(3,1,3);
x22 = myChebyshev(20);%n=20
y22 = exp(1).^abs(x22);
Ln4 = myLagrange(x22,y22,u_x);
Rn_4 = myRn_solve(u_x,x22);
plot(u_x,Ln4,'-b');
% title('Chebyshev插值 N=20');

legend('原函數','Chebyshev插值 N=10','Chebyshev插值 N=20')

 

math_test_plot_Rn.m

u_x = -1:0.01:1;
u_y = exp(1).^abs(u_x);  % y = e^|x|



% 均勻間距插值

x11 = -1:0.2:1; %n=10
y11 = exp(1).^abs(x11);
Ln1 = myLagrange(x11,y11,u_x);
% Rn_1 = exp(1)*myRn_solve(u_x,x11);
Rn_1 = myRn_solve(u_x,x11);



x12 = -1:0.1:1; %n=20
y12 = exp(1).^abs(x12);
Ln2 = myLagrange(x12,y12,u_x);
Rn_2 = myRn_solve(u_x,x12);





% 切比雪夫插值，使用切比雪夫節點作為插值節點

x21 = myChebyshev(10);%n=10
y21 = exp(1).^abs(x21);
Ln3 = myLagrange(x21,y21,u_x);
Rn_3 = myRn_solve(u_x,x21);



x22 = myChebyshev(20);%n=20
y22 = exp(1).^abs(x22);
Ln4 = myLagrange(x22,y22,u_x);
Rn_4 = myRn_solve(u_x,x22);



%繪制插值余項的最大值
plot(u_x,Rn_1,'-k',u_x,Rn_3,'-r');
title('N=10時 插值誤差(余項)的對比');
legend('均勻間距插值','Chebyshev插值');

% %繪制插值余項的最大值
% plot(u_x,Rn_2,'-k',u_x,Rn_4,'-r');
% title('N=20時 插值誤差(余項)的對比');
% legend('均勻間距插值','Chebyshev插值');

 

myChebyshev.m

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

 

myLagrange.m

function y = myLagrange(x0,y0,x)
%f11 = myLagrange(x11,y11,u_x);
% x11 = -1:0.2:1; %n=10
% y11 = exp(1).^abs(x11);


u_N=1:length(x0); 



y=zeros(size(x));

for i=u_N
    u_ij=find(u_N~=i);   
    y1=1;
    for j=1:length(u_ij)
        y1=y1.*(x-x0(u_ij(j)));  
    end
    
    y2=1;
    for j=1:length(u_ij)
        y2=y2.*(x0(i)-x0(u_ij(j)));  
    end
    
    y=y+y0(i)*y1/y2;

end

 

myRn_solve.m

function u_Rn = myRn_solve(x,x0)

% x11 = -1:0.2:1; %n=10
% y11 = exp(1).^abs(x11);
% f11 = myLagrange(x11,y11,u_x);
% Rn_1 = exp(1)*myRn_solve(u_x,x11);

u_Wn = zeros(1,length(x));
u_Rn = zeros(1,length(x));
for i = 1:length(x)
    u_Wn(i) = prod(x(i)*ones(1,length(x0))-x0);
end

u_Rn = u_Wn*exp(1)/factorial(length(x0));
end