基于Mathematica軟件在常微分方程初值問題中的可視化

孔祥強

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

?

基于Mathematica軟件在常微分方程初值問題中的可視化

孔祥強

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

基于Mathematica軟件的功能特點，針對常微分方程數值解的3種常見算法，本文通過編程實現了在同一界面下、不同算法之間的動態可視化操作，直觀地驗證了3種方法的優劣，并通過具體的數值算例驗證了操作的可行性。

Mathematica軟件；可視化；數值解；初值問題；常微分方程

Mathematica是一款有關科學計算的軟件，它將符號計算、圖形描繪和數值計算有效地結合起來.該軟件有如下特點：(1)簡單易學，操作方便；(2)強大的繪圖和動態顯示功能；(3)編程語言精練，效率高；(4)與其它數學軟件可互相調用[1].Mathematica軟件廣泛應用于圖形繪制、動態模擬、隨機模擬和物理現象模擬等諸多方面[2].常微分方程是一門較抽象的學科，其中的一些算法不好理解，比如求微分方程數值解的幾種算法，利用Mathematica軟件可使這些復雜的問題簡單化、直觀化.通過Mathematica軟件，可實現方程數值解的可視化操作，使多種不同的計算方法在同一界面下顯示出來；通過參數的調整，由圖像直觀觀察出數值解與解析解的誤差，從而對算法的優劣進行判定，達到抽象內容具體化、教學內容直觀化和提高教學質量的目的.

1 常微分方程初值問題的3種常用解法

一般地有yn+1=yn+hf(xn,yn),n=0,1,2,…，N-1.

稱該式為Euler法公式.Euler法的局部截斷誤差的階為O(h2)，為一階精度的方法[4].Euler法有明顯的幾何意義.以f(x0,y0)為斜率，過點(x0,y0)做直線，它與直線x=x1的交點為y1，依次類推，yn+1是以f(xn,yn)為斜率，過點(xn,yn)的直線與直線x=xn+1的交點，故Euler法也稱為折線法[5].

稱上式為梯形法公式，梯形法為二階精度的方法[3].

Runge-Kutta法是一種非線性的高階單步法，該方法的導出基于Taylor展開，故要求所求的解有較高的光滑度，對于光滑性不太好的解，最好采用低階算法而將步長取小.經典的四階Runge-Kutta法的計算公式為

其局部截斷誤差為O(h5)，為四階精度的方法，階數越高，方法越精確[5].

2 3種算法的動態可視化實現

Euler法、梯形法和Runge-Kutta法均可求方程的數值解，下面通過數值算例動態演示算法的解題過程，從中了解3種算法的優劣.

若用3種方法單獨求每一個方程的初值問題比較繁瑣，計算量大，且不易比較3種方法的優缺點.下面通過Mathematica軟件編程，可動態觀察3種算法的計算過程，從而實現以下功能：(1)通過下拉菜單選擇不同的方程,再進一步動態選擇Euler法、梯形法或Runge-Kutta法；(2)作出微分方程解析解的圖像及數值解的折線圖；(3)對每一個方程，可以選擇不同的步長h(不妨設0.1≤h≤2)，且能在圖像上反映出節點個數，并顯示出h的具體值；(4)通過同一窗口下直觀比較數值解和解析解的誤差，對3種方法進行比較；(5)通過微調x軸、y軸，選擇不同的視角觀察圖形.

實現以上功能的源程序：

f1[x_,y_]:=y-2.01x;f2[x_,y_]:=-1.01y^2;

f3[x_,y_]:=-1.01x^2y;f4[x_,y_]:=x^2-2.01y;

f5[x_,y_]:=Sin[1.01x]y;f6[x_,y_]:=(x-x^2)y;

p[x_,y_]:={x,y};

q[x_List,y_List]:=MapThread[p,{x,y}];

jie[a_,s_,y0_,t_]:=DSolve[{y’[t]s[t,y[t]],y[a]y0},y,t];

dian={ImageSize450,PlotRangePadding0.41};

s1[{a_,b_},s_,y0_,n_]:=Module[{h,x,y,u},h=(b-a)/n;

x=Table[a+(j-1)h,{j,1,n+1}];

y=Table[0,{n+1}];y[[1]]=y0;

Do[y[[j]]=y[[j-1]]+hs[x[[j-1]],y[[j-1]]],{j,2,n+1}];u=q[x,y]];

t1[{y0_,x0_,h_,s_},{D_,b_,d_}]:=Module[{g,z},g=Length[b];

z=Table[0,{j,1,g}];z[[1]]=s[x0,y0];

Do[z[[j]]=s[x0+d[[j]]h,y0+hSum[D[[j-1,k]]z[[k]],{k,1,j-1}]],{j,2,g}];

yn=y0+hSum[b[[j]]z[[j]],{j,1,g}]];

t2[{a_,r_},s_,y0_,n_,{D_,b_,d_}]:=Module[{h,x,y,v,j},h=(r-a)/n;

x=Table[a+(j-1)h,{j,1,n+1}];

y=Table[0,{n+1}];y[[1]]=y0;

Do[y[[j]]=t1[{y[[j-1]],x[[j-1]],h,s},{D,b,d}],{j,2,n+1}];

v=q[x,y]];

s2[{a_,r_},s_,y0_,n_]:=Module[{D,b,d},D={{1/2}};

b={0,1};d={0,1/2};

t2[{a,r},s,y0,n,{D,b,d}]];

s3[{a_,r_},s_,y0_,n_]:=Module[{D,b,d},D={{1/2},{0,1/2},{0,0,1}};

b={1/6,1/3,1/3,1/6};d={0,1/2,1/2,1};t2[{a,r},s,y0,n,{D,b,d}]];

Manipulate[Show[{Plot[Evaluate[y[t]/.jie[a,s,d,t]],{t,a,b+0.1},

PlotStyle{Orange,Thick},AxesOrigin{0,0},PlotRangeAll,

PlotLabelRow[{

Style["步長:",16,Bold,Red],Style[N[(b-a)/qq,12],16,Brown]}]],

ListPlot[Tooltip[Check[fangfa[{a,b},s,d,If[MemberQ[{f3,f5,f6},s],N[qq],qq]],{a}]],PlotStyle{Blue,PointSize[0.016]}],

ListLinePlot[Check[fangfa[{a,b},s,d,If[MemberQ[{f3,f5,f6},s],

N[qq],qq]],{a,d}],PlotStyle{Blue,Dashed},PlotRangeAll]},

ImageSize600,AxesStyleArrowheads[0.02],

Epilog

{Text[Framed[Column[{Row[{Style["-",Red,Bold,18],

Style["方程的解析解",14,Bold,Magenta]}],

Row[{Style["-",Blue,Bold,18],Style["方程的數值解",14,Bold,Magenta]}]}]],{1,0.8},{1,0.8}]}],{{s,f1,

Style["選擇方程y'=",13,Blue,Bold]},{f1"y-2.01x",

f2"-1.01y2",f3"-1.01x2y",f4"x2-2.01y",

f5"Sin[1.01x]y",f6"(x-x2)y"},

ControlTypePopupMenu,ControlPlacementBottom},

{{fangfa,s1,Style["選擇方法",13,Brown,Bold]},{s1"Euler法",

s2 "梯形法",s3 "Runge-Kutta法"},ControlPlacementBottom},

Delimiter,{{qq,5,

Style["選擇步長",13,Blue,Bold]},1,20,1,ControlPlacementBottom,

Appearance"Labeled"},Delimiter,{{a,0,

Style["向左微調x軸",13,Bold,Green]},-1/3,1/3,1/5,

Appearance"Labeled",ControlPlacementBottom},{{b,2,

Style["向右微調x軸",13,Bold,Green]},1,3.5,1/4,

ControlPlacementBottom,Appearance"Labeled"},

{{d,1,Style["上下微調y軸",13,Bold,Green]},1/2,5/2,1/5,

ControlPlacementBottom,Appearance"Labeled"},

TrackedSymbolsManipulate,SaveDefinitionsTrue]

輸出結果：

圖1 用Euler法求初值問題y′=y-2.01x,y(0)=1

圖1中，連續曲線為方程解析解的圖形，折線為數值解圖像；步長為0.400000000000，通過單擊“動畫控件”按鈕，可實現步長的自動播放功能.通過窗口最下方的3個微調按鈕，可選擇合適的視角觀察圖形.

圖2 用梯形法求初值問題，y′=y-2.01x,y(0)=1

圖2選擇的是梯形法，與圖1比較，在同步長的情形下，梯形法所得數值解產生的誤差明顯小于Euler法，梯形法是比Euler法更精確的方法.在調節步長的同時，節點也隨時變化.隨著步長的減小，圖形的節點動態增加.節點個數越多，所得數值解折線和解析解曲線越吻合.

圖3采用Runge-Kutta法解同一個初值問題，與圖1和圖2比較可知，該法所得數值解曲線和解析解曲線基本一致，是3種方法里誤差最小的，是比Euler法和梯形法都要精確的方法.

圖4 用Euler法求初值問題y′=x2-2.01y,y(0)=1

圖4中選擇的方程為y′=x2-2.01y，通過下拉菜單也可以選擇y′=-1.01y2或y′=-1.01x2y或y′=sin(1.01x)·y或y′=(x-x2)y，這樣求解速度快，大大提高了做題效率，實用性強.案例中的微分方程，完全可以更換為其它方程，只需將源程序中的函數表達式改為其它表達式即可，因此程序具有一般性.

3 結語

通過案例可看出，利用Mathematica軟件編程，可達到將復雜問題簡單化、直觀化的目的[6].在教學中，能夠突破教學難點，激發學生的學習興趣，使他們更好地理解和認識所學的各種算法，以達到提高解決實際問題能力的目的.Mathematica軟件的使用為數學教學注入了生命力，使數學軟件和數學教學有機地結合起來，成為一種新型的教學方法和手段.

[1]王紹恒,王藝靜.Mathematica軟件在大學數學課程教學中的應用[J].教育理論與實踐,2013,33(21):39-40.

[2]李士恒.Mathematica在數學實驗中的應用[J].中央民族大學學報:自然科學版,2014,23(3):14-17.

[3]關治,陳景良.數值計算方法[M].北京:清華大學出版社,2005:236-242.

[4]李榮華,劉播.微分方程數值解法[M].北京:高等教育出版社,2011:1-37.

[5]張韻華,奚梅成,陳效群編.數值計算方法與算法[M].北京:科學出版社,2006:168-186.

[6]程春蕊,朱軍輝.Mathematica軟件在常微分方程教學中的應用[J].高等函授學報:自然科學版,2012,25(1):21-23.

Mathematica-based Visual Insight in the Initial Value Problems of Ordinary Differential Equation

KONG Xiang-qiang

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

Based on the functions and characteristics of the Mathematica software, in view of the numerical solution of three common algorithms of ordinary differential equations,we realized the dynamic visualization operation under the same interface between different algorithms through the programming. The advantages and disadvantages of three kinds of method were verified, and we verified the feasibility of operation by detailed numerical examples.

Mathematica software; visualization; numerical solution; initial value problems;ordinary differential equation

2015-07-17

菏澤學院重點課題組項目(201311)。

孔祥強(1983- )，男，山東菏澤人，菏澤學院數學系講師，碩士，從事計算數學研究。

O245

A

2095-7602(2015)10-0020-06