孔祥强
(菏泽学院数学系,山东菏泽 274015)
Mathematica软件在数值分析中的交互可视化设计
孔祥强
(菏泽学院数学系,山东菏泽 274015)
[摘要]本文在Mathematica 9.0.1软件环境下,实现了求解线性方程组的雅可比迭代法和高斯-赛德尔迭代法的动态交互可视化设计,通过实例对两种方法进行了比较,设计了利用插值型公式求数值积分的人机交互化过程,并通过算例验证了程序的正确性。
[关键词]Mathematica软件;交互可视化;迭代法;线性方程组;数值积分
Mathematica软件是Wolfram公司开发的交互式科学计算软件,具有强大的符号计算功能和绘图功能,特别是2013年新推出的Mathematica 9.0.1版本,更加完善了这一软件的交互可视化功能.利用Mathematica软件,可将一些复杂的、晦涩的或者是抽象的数值算法直观、动态地表达出来,而这些表达是可以控制和人工交互的.软件的这些特性在概率论与数理统计和高等数学中已经得到了很好的应用[1-3],但在数值分析课程中应用的还比较少.Mathematica提供了众多可进行可视化设计的函数[4-7],如Animate函数和Manipulate函数等.利用这些函数,可实现对一个或多个参数的动态控制,通过“显示动画按钮”控件实现计算过程的动态展示.
线性方程组和数值积分是数值分析课程的重要内容.求解线性方程组的方法有雅可比(Jacobi)迭代法、高斯-赛德尔(G-S)迭代法和逐次超松弛(SOR)迭代法等.求数值积分的方法有插值型(Newton-Cotes)公式法、Gauss型公式法和复化型公式法等[8].要理解和掌握以上理论知识并不容易,关键问题是计算公式繁杂和计算过程的不可视化.本文给出了Jacobi迭代法和G-S迭代法的可视化设计,在同一界面下选择不同的方法,实现迭代次数和所得近似解的展示,并且可比较出两种方法的优劣.另外,实现了用四种不同的Newton-Cotes公式计算数值积分的动态化,依据选择的步长不同达到选取不同公式的目的,而所得的近似值和误差值实时变化.
1Jacobi迭代法和G-S迭代法的可视化实现
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,则构造的一阶定常迭代法可表示为
(*)
其中,矩阵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迭代法的动态模拟
首先判断出迭代矩阵的谱半径分别为0.506079和0.2,均是小于1的.由迭代法基本定理可知,Jacobi迭代法和G-S迭代法都是收敛的.
比较以上两种迭代法,当A满足一定条件,初始向量x(0)相同,且达到同样精度的情形下,Jacobi迭代法需迭代18次才能达到方程组的精确解{-4,3,2},而G-S迭代法仅需9次,可知G-S迭代法比Jacobi迭代法收敛速度快.
1.3程序实现
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["迭代次数:",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["迭代次数:",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迭代法进行编程,计算出迭代矩阵、谱半径、迭代次数和近似解,然后对G-S迭代法进行编程,最后通过Manipulate函数使得两种方法显示在同一界面内,实现了动态交互的一体化功能. 2Newton-Cotes公式的可视化实现 2.1基本知识 2.2可视化设计 通过界面左下方的“步长选择”控件可选择不同的公式,而选取的公式名称出现在界面上方.也可以通过点击控件右方的“+显示动画控制”,启动动画的自动播放模式,设定按钮向前、向后移动,控制播放的快慢程度等.界面中表格内的“近似值”和“误差”数据随着公式的不同而动态变化,“真实值”保持不变.对于n=4时的Cotes公式,通过表格的形式给出了其系数的取值.界面中粗红线显示的是函数y=1-e-xsin3x的图形,插值节点用红色圆点表示,连接节点的曲线下方用浅红色进行填充,表示通过公式计算所得到的近似值,即曲边梯形的面积.本文中的四个Newton-Cotes公式都不是高阶的,可以证明,当n≤7时,公式是稳定的;当n>8时,稳定性不能保证;当n→∞时,容易产生Runge现象,收敛性也无法保证,因此,一般不采用高阶的Newton-Cotes公式求数值积分[8]. 图2 Newton-Cotes公式求数值积分的动态模拟 2.3程序实现 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公式系数表",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.在上面的程序中,首先定义出函数、积分区间和Cotes系数表等;通过Plot函数及Graphics函数绘制g1,g2,g3的图形,利用Show函数将绘制的图形置于同一界面下;最后通过Manipulate函数实现动态交互的一体化功能. 3结语 将可运行Notebook文档保存为.nb格式文件,既可在Mathematica 9.0.1软件环境下运行,也可以在Wolfram公司提供的Mathematica Player 9.0.1下运行.本文对Jacobi迭代法及G-S迭代法进行了可视化设计,实际上,逐次超松弛(SOR)迭代法的可视化也可参考本文提供的编程方法,但必须通过控件控制松弛因子ω∈(0,2),SOR迭代法才有可能收敛.当ω≤0或者ω≥2时,SOR迭代法一定是不收敛的.同样,利用Gauss型公式和复化型公式求数值积分时,也可以给出类似Newton-Cotes公式的可视化设计. Mathematica是集数学计算、数据的处理与分析为一体的软件,对科学实验和科学研究起到了很大的促进作用,比如该软件在微分方程[9]、微积分教学[10]和实验数据分析[11]中都有广泛应用.利用Mathematica软件,结合编程,可开发出适合数值分析课程的教学交互可视化课件,对提高教学质量和培养学生学习的积极性都是大有裨益的. [参考文献] [1]郑永凡,王艳青.基于Mathematica的交互动态可视化设计及其应用[J].辽宁大学学报:自然科学版,2010,37(4):324-328. [2]谷创业,王春峰.Mathematica在互动式多媒体数学教学中的应用[J].中国现代教育装备,2010(5):76-68. [3]徐安农.计算机模拟在概率论课程教学中的应用[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]朱晓峰,田益民.数学软件在数学教学中的应用——数学实验一例[J].数学的实践与认识,2009,39(13):232-235. [8]易大义,陈道琦.数值分析引论[M].浙江:浙江大学出版社,2002. [9]秦大康,张宝善.Mathematica软件在微分方程中的应用研究[J].徐州师范大学学报:自然科学版,2001,19(3):9-11. [10]马千里.Mathematica软件在微积分教学中的应用[J].陕西师范大学学报:自然科学版,2009,37(s1):135-136. [11]金哲,张凤兰.Mathematica与实验数据分析[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- ),男,讲师,硕士,从事应用数学研究。 [基金项目]菏泽学院重点课题组项目“关于Mathematica软件在‘高等数学’教学中的交互可视化研究”(2015010);山东省教育科学“十二五”规划“高等教育数学教学专项”重点课题“基于数学软件的大学数学教学的交互可视化研究”(ZBS15004)。 [收稿日期]2015-11-13