Mathematica軟件在數(shù)值分析中的交互可視化設計

孔祥強

(菏澤學院數(shù)學系，山東菏澤 274015)

?

Mathematica軟件在數(shù)值分析中的交互可視化設計

孔祥強

(菏澤學院數(shù)學系，山東菏澤 274015)

[摘要]本文在Mathematica 9.0.1軟件環(huán)境下，實現(xiàn)了求解線性方程組的雅可比迭代法和高斯-賽德爾迭代法的動態(tài)交互可視化設計，通過實例對兩種方法進行了比較，設計了利用插值型公式求數(shù)值積分的人機交互化過程，并通過算例驗證了程序的正確性。

[關(guān)鍵詞]Mathematica軟件；交互可視化；迭代法；線性方程組；數(shù)值積分

Mathematica軟件是Wolfram公司開發(fā)的交互式科學計算軟件，具有強大的符號計算功能和繪圖功能，特別是2013年新推出的Mathematica 9.0.1版本，更加完善了這一軟件的交互可視化功能.利用Mathematica軟件，可將一些復雜的、晦澀的或者是抽象的數(shù)值算法直觀、動態(tài)地表達出來，而這些表達是可以控制和人工交互的.軟件的這些特性在概率論與數(shù)理統(tǒng)計和高等數(shù)學中已經(jīng)得到了很好的應用[1-3]，但在數(shù)值分析課程中應用的還比較少.Mathematica提供了眾多可進行可視化設計的函數(shù)[4-7]，如Animate函數(shù)和Manipulate函數(shù)等.利用這些函數(shù)，可實現(xiàn)對一個或多個參數(shù)的動態(tài)控制，通過“顯示動畫按鈕”控件實現(xiàn)計算過程的動態(tài)展示.

線性方程組和數(shù)值積分是數(shù)值分析課程的重要內(nèi)容.求解線性方程組的方法有雅可比(Jacobi)迭代法、高斯-賽德爾(G-S)迭代法和逐次超松弛(SOR)迭代法等.求數(shù)值積分的方法有插值型(Newton-Cotes)公式法、Gauss型公式法和復化型公式法等[8].要理解和掌握以上理論知識并不容易，關(guān)鍵問題是計算公式繁雜和計算過程的不可視化.本文給出了Jacobi迭代法和G-S迭代法的可視化設計，在同一界面下選擇不同的方法，實現(xiàn)迭代次數(shù)和所得近似解的展示，并且可比較出兩種方法的優(yōu)劣.另外，實現(xiàn)了用四種不同的Newton-Cotes公式計算數(shù)值積分的動態(tài)化，依據(jù)選擇的步長不同達到選取不同公式的目的，而所得的近似值和誤差值實時變化.

1Jacobi迭代法和G-S迭代法的可視化實現(xiàn)

1.1基本知識

設Ax=b，A∈Rn×n且|A|≠0，對A進行如下分解：

若將A分裂為A=M-N(其中M為近似于A的矩陣，且|M|≠0，稱為分裂陣)，則Ax=b?Mx=Nx+b，

即x=M-1Nx+M-1b，則構(gòu)造的一階定常迭代法可表示為

(*)

其中，矩陣B為迭代矩陣.

若取M=D，此時A=D-N，可得Ax=b的Jacobi迭代法，此時式(*)中的迭代矩陣B=I-D-1A=D-1(L+U)，f=D-1b.

若取M=D-L，此時A=M-N，可得Ax=b的G-S迭代法，此時式(*)中的迭代矩陣B=I-(D-L)-1A=(D-L)-1U，f=(D-L)-1b[8].

迭代法基本定理[8](1)設有方程組x=Bx+f,B∈Rn×n，(2)有迭代法x(k+1)=Bx(k)+f.對任意的初始向量x(0)，迭代法收斂的充要條件是B的所有特征值λi(B)滿足|λi(B)|<1或譜半徑ρ(B)<1.

1.2可視化設計

圖1　Jacobi迭代法和G-S迭代法的動態(tài)模擬

首先判斷出迭代矩陣的譜半徑分別為0.506079和0.2，均是小于1的.由迭代法基本定理可知，Jacobi迭代法和G-S迭代法都是收斂的.

比較以上兩種迭代法，當A滿足一定條件，初始向量x(0)相同，且達到同樣精度的情形下，Jacobi迭代法需迭代18次才能達到方程組的精確解{-4,3,2}，而G-S迭代法僅需9次，可知G-S迭代法比Jacobi迭代法收斂速度快.

1.3程序?qū)崿F(xiàn)

Manipulate[

If[t"Jacobi迭代法",A={{5,2,1},{-1,4,2},{2,-3,10}};

b={-12,20,3};x0={0,0,0};epi=0.0001;

F=DiagonalMatrix[{A[[1,1]],A[[2,2]],A[[3,3]]}];

J=DiagonalMatrix[{1,1,1}];

B=J-Inverse[F].A;

P=Max[Abs[Eigenvalues[N[B]]]];

If[P<1,f=Inverse[F].b;

jk=1;x1=B.x0+f;S[1]=x1;

While[Norm[x0-x1]>epi,x0=x1;

x1=B.x0+f;S[jk+1]=x1;jk=jk+1;

If[jk>50,Break[]]];

Pane[Column[{

Style[Row[{"",MatrixForm[A],"X=",MatrixForm[b]}],15,Bold],,

,Style[Row[{ Style["Jacobi迭代矩陣:",Red,Bold

,15],MatrixForm[B]}],15],,

,Framed@Style[Column[{

Row[{Style["譜半徑:",Blue,Bold,15],P}]

,Row[{Style["迭代次數(shù):",Blue,Bold,15],jk,

,Tooltip[

Graphics[{Red,Disk[{0,0},1,{Pi/4,3Pi/4}]},ImageSize 15],

TableForm[Table[ S[i]//N,{i,jk}]

,TableHeadings{Table[Row[{"X", i}] ,{i,jk}],{"近似解"}}]]}]}],15]

,Framed@Style[Row[{Style["近似解",Blue,Bold,15],N[x1]}],15]}]]],

If[t"G—S迭代法",A={{5,2,1},{-1,4,2},{2,-3,10}};

b={-12,20,3};];x0={0,0,0};epi=0.0001;

x1=M.x0+f1;

F=DiagonalMatrix[{A[[1,1]],A[[2,2]],A[[3,3]]}];

L=-Table[If[i>j,A[[i,j]],0],{i,3},{j,3}];

U=-Table[If[i

f1=Inverse[F-L].b;M=Inverse[F-L].U;

Q=Max[Abs[Eigenvalues[N[M]]]];

T[1]=x1;k=1;

While[Norm[x0-x1]>epi,x0=x1;

x1=M.x0+f1;T[k+1]=x1;k=k+1;

If[k>50,Break[]]];

Pane[Column[{

Style[Row[{"",MatrixForm[A],"X=",MatrixForm[b]}],15,Bold],,

,Style[Row[{Style["G-S迭代矩陣:",Red,Bold,15],MatrixForm[M]}]

,15],,

,Framed@Style[Column[{

Row[{Style["譜半徑:",Blue,Bold,15],Q}]

,Row[{Style["迭代次數(shù):",Blue,Bold,15],k,

,Tooltip[

Graphics[{Red,Disk[{0,0},1,{Pi/4,3Pi/4}]},ImageSize 15],

TableForm[Table[T[i]//N,{i,k}]

,TableHeadings{Table[Row[{"X",i}],{i,jk}],{"近似解"}}]]}]}],15]

,Framed@Style[Row[{Style["近似解",Blue,Bold,15],N[x1]}],15]}]]]

,{{t,"Jacobi迭代法",Row[{Style["迭代方法選擇",Blue,Bold,15]}]}

,{"Jacobi迭代法","G-S迭代法"}}

,ControlPlacement Bottom]

在上面的程序中，首先對Jacobi迭代法進行編程，計算出迭代矩陣、譜半徑、迭代次數(shù)和近似解，然后對G-S迭代法進行編程，最后通過Manipulate函數(shù)使得兩種方法顯示在同一界面內(nèi)，實現(xiàn)了動態(tài)交互的一體化功能.

2Newton-Cotes公式的可視化實現(xiàn)

2.1基本知識

2.2可視化設計

通過界面左下方的“步長選擇”控件可選擇不同的公式，而選取的公式名稱出現(xiàn)在界面上方.也可以通過點擊控件右方的“+顯示動畫控制”，啟動動畫的自動播放模式，設定按鈕向前、向后移動，控制播放的快慢程度等.界面中表格內(nèi)的“近似值”和“誤差”數(shù)據(jù)隨著公式的不同而動態(tài)變化，“真實值”保持不變.對于n=4時的Cotes公式，通過表格的形式給出了其系數(shù)的取值.界面中粗紅線顯示的是函數(shù)y=1-e-xsin3x的圖形，插值節(jié)點用紅色圓點表示，連接節(jié)點的曲線下方用淺紅色進行填充，表示通過公式計算所得到的近似值，即曲邊梯形的面積.本文中的四個Newton-Cotes公式都不是高階的，可以證明，當n≤7時，公式是穩(wěn)定的；當n>8時，穩(wěn)定性不能保證；當n→∞時，容易產(chǎn)生Runge現(xiàn)象，收斂性也無法保證，因此，一般不采用高階的Newton-Cotes公式求數(shù)值積分[8].

圖2　Newton-Cotes公式求數(shù)值積分的動態(tài)模擬

2.3程序?qū)崿F(xiàn)

Manipulate[

f[x_]:=1-Exp[-x]* Sin[3 x];

a=0;b=3;n=k;h=3/n;

c={{1/2,1/2,,,},{1/6,4/6,1/6,,}

,{1/8,3/8,3/8,1/8,},{7/90,16/45,2/15,16/45,7/90}};

data=Table[{a+j h,f[a+j h]}, {j,0,n,1}];

table=Grid[{

{Text@Style["Cotes公式系數(shù)表",15,Bold,Blue]}

,{Style[Text@Grid[c,Frame {None,All}],Brown,17]}}];g1=Plot[Evaluate[

Fit[data,Table[x^j,{j,0,n}],x]],{x,0,3}

,PlotRange {{-0.11,3.03},{0,3.1}}

,PlotStyle Magenta,PlotLabel Switch[k,

1,Style["梯形公式(n=1時)",15,Bold,Blue],

2,Style["Simpson公式(n=2時)",15,Bold,Blue],

3,Style["3/8 Simpson公式(n=3時)",15,Bold,Blue],

4,Style["Cotes公式(n=4時)",15,Bold,Blue]]

,AxesStyle Arrowheads[0.02]

,FillingStyle LightRed, Filling Axis

, Epilog {Text[Style[Grid[{

{Style["真實值",Red,14],zs}

,{Style["近似值",Red,14],js}

,{Style["誤差值",Red,14],zs-js}}

,Frame All],14],{2.2,2.3},{0,0}]} ] ;

g2=Graphics[{PointSize[0.023],Red

,Table[Point[data[[j]]],{j,1,n+1,1}]}];

g3=Plot[f[x],{x,0,3},PlotStyle {Thick,Red}];

Grid[{{Show[{g1,g2,g3},ImageSize 500],table}}]

,{{k,1,Style["步長選擇:",15,Bold,Blue]},1,4,1}

,ControlPlacement Bottom]

選擇n=2時的Simpson公式，可得積分的近似值為2.28398255.在上面的程序中，首先定義出函數(shù)、積分區(qū)間和Cotes系數(shù)表等；通過Plot函數(shù)及Graphics函數(shù)繪制g1,g2,g3的圖形，利用Show函數(shù)將繪制的圖形置于同一界面下；最后通過Manipulate函數(shù)實現(xiàn)動態(tài)交互的一體化功能.

3結(jié)語

將可運行Notebook文檔保存為.nb格式文件，既可在Mathematica 9.0.1軟件環(huán)境下運行，也可以在Wolfram公司提供的Mathematica Player 9.0.1下運行.本文對Jacobi迭代法及G-S迭代法進行了可視化設計，實際上，逐次超松弛(SOR)迭代法的可視化也可參考本文提供的編程方法，但必須通過控件控制松弛因子ω∈(0,2)，SOR迭代法才有可能收斂.當ω≤0或者ω≥2時，SOR迭代法一定是不收斂的.同樣，利用Gauss型公式和復化型公式求數(shù)值積分時，也可以給出類似Newton-Cotes公式的可視化設計.

Mathematica是集數(shù)學計算、數(shù)據(jù)的處理與分析為一體的軟件，對科學實驗和科學研究起到了很大的促進作用，比如該軟件在微分方程[9]、微積分教學[10]和實驗數(shù)據(jù)分析[11]中都有廣泛應用.利用Mathematica軟件，結(jié)合編程，可開發(fā)出適合數(shù)值分析課程的教學交互可視化課件，對提高教學質(zhì)量和培養(yǎng)學生學習的積極性都是大有裨益的.

[參考文獻]

[1]鄭永凡,王艷青.基于Mathematica的交互動態(tài)可視化設計及其應用[J].遼寧大學學報:自然科學版,2010,37(4):324-328.

[2]谷創(chuàng)業(yè),王春峰.Mathematica在互動式多媒體數(shù)學教學中的應用[J].中國現(xiàn)代教育裝備,2010(5):76-68.

[3]徐安農(nóng).計算機模擬在概率論課程教學中的應用[J].桂林電子科技大學學報,2008,28(2):153-156.

[4]Ruskeep?? Heikki. Mathematica Navigator:Mathematics,Statistics,and Graphics[M].3rd Edition,Elsevier Inc.,2009.

[5]Sal Mangano.Mathematica Cookbook[M].O’Reilly Media,2010.

[6]吳飛.Mathematica項目演示筆記[M].北京:清華大學出版社,2010.

[7]朱曉峰,田益民.數(shù)學軟件在數(shù)學教學中的應用——數(shù)學實驗一例[J].數(shù)學的實踐與認識,2009,39(13):232-235.

[8]易大義,陳道琦.數(shù)值分析引論[M].浙江:浙江大學出版社,2002.

[9]秦大康,張寶善.Mathematica軟件在微分方程中的應用研究[J].徐州師范大學學報:自然科學版,2001,19(3):9-11.

[10]馬千里.Mathematica軟件在微積分教學中的應用[J].陜西師范大學學報:自然科學版,2009,37(s1):135-136.

[11]金哲,張鳳蘭.Mathematica與實驗數(shù)據(jù)分析[J].延邊大學學報:自然科學版,2003,29(1):36-39.

Mathematica Based Design of Interactive Visualization in Numerical Analysis

KONG Xiang-qiang

(Department of Mathematics, Heze University, Heze Shandong 274015,China)

Abstract:Based on the Mathematica 9.0.1 software,the dynamic interactive visual design to solve linear equations by the Jacobi iteration and Gauss-Seidel is realized.And two methods were compared by a arithmetic example ,the human-computer interaction process of solving numerical integration by using the type of interpolation formula is designed , and the program is validated the correctness by example.

Key words:Mathematica software;interactive visualization;iterative method;linear equations;numerical integration

[中圖分類號]O245

[文獻標識碼]A

[文章編號]2095-7602(2016)04-0043-06

[作者簡介]孔祥強(1983- )，男，講師，碩士，從事應用數(shù)學研究。

[基金項目]菏澤學院重點課題組項目“關(guān)于Mathematica軟件在‘高等數(shù)學’教學中的交互可視化研究”(2015010)；山東省教育科學“十二五”規(guī)劃“高等教育數(shù)學教學專項”重點課題“基于數(shù)學軟件的大學數(shù)學教學的交互可視化研究”(ZBS15004)。

[收稿日期]2015-11-13