孔祥強(qiáng)

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

?

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

孔祥強(qiáng)

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

基于Mathematica軟件的功能特點(diǎn)，針對(duì)常微分方程數(shù)值解的3種常見算法，本文通過編程實(shí)現(xiàn)了在同一界面下、不同算法之間的動(dòng)態(tài)可視化操作，直觀地驗(yàn)證了3種方法的優(yōu)劣，并通過具體的數(shù)值算例驗(yàn)證了操作的可行性。

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

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

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

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

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

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

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

其局部截?cái)嗾`差為O(h5)，為四階精度的方法，階數(shù)越高，方法越精確[5].

2 3種算法的動(dòng)態(tài)可視化實(shí)現(xiàn)

Euler法、梯形法和Runge-Kutta法均可求方程的數(shù)值解，下面通過數(shù)值算例動(dòng)態(tài)演示算法的解題過程，從中了解3種算法的優(yōu)劣.

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

實(shí)現(xiàn)以上功能的源程序：

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["步長(zhǎng):",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["方程的數(shù)值解",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["選擇步長(zhǎng)",13,Blue,Bold]},1,20,1,ControlPlacementBottom,

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

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

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

Style["向右微調(diào)x軸",13,Bold,Green]},1,3.5,1/4,

ControlPlacementBottom,Appearance"Labeled"},

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

ControlPlacementBottom,Appearance"Labeled"},

TrackedSymbolsManipulate,SaveDefinitionsTrue]

輸出結(jié)果：

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

圖1中，連續(xù)曲線為方程解析解的圖形，折線為數(shù)值解圖像；步長(zhǎng)為0.400000000000，通過單擊“動(dòng)畫控件”按鈕，可實(shí)現(xiàn)步長(zhǎng)的自動(dòng)播放功能.通過窗口最下方的3個(gè)微調(diào)按鈕，可選擇合適的視角觀察圖形.

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

圖2選擇的是梯形法，與圖1比較，在同步長(zhǎng)的情形下，梯形法所得數(shù)值解產(chǎn)生的誤差明顯小于Euler法，梯形法是比Euler法更精確的方法.在調(diào)節(jié)步長(zhǎng)的同時(shí)，節(jié)點(diǎn)也隨時(shí)變化.隨著步長(zhǎng)的減小，圖形的節(jié)點(diǎn)動(dòng)態(tài)增加.節(jié)點(diǎn)個(gè)數(shù)越多，所得數(shù)值解折線和解析解曲線越吻合.

圖3采用Runge-Kutta法解同一個(gè)初值問題，與圖1和圖2比較可知，該法所得數(shù)值解曲線和解析解曲線基本一致，是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，這樣求解速度快，大大提高了做題效率，實(shí)用性強(qiáng).案例中的微分方程，完全可以更換為其它方程，只需將源程序中的函數(shù)表達(dá)式改為其它表達(dá)式即可，因此程序具有一般性.

3 結(jié)語

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

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

[2]李士恒.Mathematica在數(shù)學(xué)實(shí)驗(yàn)中的應(yīng)用[J].中央民族大學(xué)學(xué)報(bào):自然科學(xué)版,2014,23(3):14-17.

[3]關(guān)治,陳景良.數(shù)值計(jì)算方法[M].北京:清華大學(xué)出版社,2005:236-242.

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

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

[6]程春蕊,朱軍輝.Mathematica軟件在常微分方程教學(xué)中的應(yīng)用[J].高等函授學(xué)報(bào):自然科學(xué)版,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

菏澤學(xué)院重點(diǎn)課題組項(xiàng)目(201311)。

孔祥強(qiáng)(1983- )，男，山東菏澤人，菏澤學(xué)院數(shù)學(xué)系講師，碩士，從事計(jì)算數(shù)學(xué)研究。

O245

A

2095-7602(2015)10-0020-06