差分方程实验

差分方程实验第一页，共二十九页，编辑于2023年，星期日4.2日常生活中的经济问题银行存款与利率假如你在银行开设了一个1000元的存款账户，银行的年利率为7%.用an表示n年后你账户上的存款额，那么下面的数列就是你每年的存款额：

a0,a1,a2,a3,…,an,…

设r为年利率，由于an+1=an+ran,因此存款问题的数学模型是：

a0=1000,an+1=(1+r)an,n=1,2,3,…第二页，共二十九页，编辑于2023年，星期日家庭教育基金从1994年开始，我国逐步实行了大学收费制度.为了保障子女将来的教育费用，小张夫妇从他们的儿子出生时开始，每年向银行存入x元作为家庭教育基金.若银行的年利率为r，试写出第n年后教育基金总额的表达式.预计当子女18岁入大学时所需的费用为30000元，按年利率10%计算，小张夫妇每年应向银行存入多少元?设n年后教育基金总额为an，每年向银行存入x元，依据复利率计算公式，得到家庭教育基金的数学模型为：

a0=x,an+1=(1+r)an+x,n=0,1,2,3,…第三页，共二十九页，编辑于2023年，星期日抵押贷款小李夫妇要购买二居室住房一套，共需10万元.他们已经筹集4万元，另外6万元申请抵押贷款.若贷款月利率为1%，还贷期限为25年，问小李夫妇每月要还多少钱？设贷款额为a0，每月还贷额为x，月利率为r，第n个月后的欠款额为an，则

a0=60000,a1=(1+r)a0-x,a2=(1+r)a1-x,……an=(1+r)an-1-x,n=1,2,3,…第四页，共二十九页，编辑于2023年，星期日分期付款小王看到一则广告：商场对电脑实行分期付款销售.一台售价8000元的电脑，可分36个月付款，每月付300元即可.同时他收到了银行提供消费贷款的消息：10000元以下的贷款，可在三年内还清，年利率为15%.那么，他买电脑应该向银行贷款，还是直接向商店分期付款？经过分析可知，分期付款与抵押贷款模型相同.设第n个月后的欠款额为an，则

a0=8000,an+1=(1+r)an-300,n=0,1,2,3,…

贷款模型

a0=8000,an+1=(1+0.15/12)an-x,n=0,1,2,3,…第五页，共二十九页，编辑于2023年，星期日一阶线性差分方程在上述模型中，给出了an+1与an之间的递推公式.将它们写成统一的形式：

a0=c,an+1=an+b,n=0,1,2,3,…称此类递推关系为一阶线性差分方程.当b=0时称为齐次差分方程，否则称为非齐次差分方程.

定义1

对任意数列A={a1,a2,…,an,…}，其差分算子定义如下：a1=a2-a1,a2=a3-a2,…an=an+1-an,…

定义2

对数列A={a1,a2,…,an,…}，其一阶差分的差分称为二阶差分,记为2A=(A).即：2an=an+1-an=(an+2-an+1)-(an+1-an)=an+2-2an+1+an一般地，可以定义n阶差分.第六页，共二十九页，编辑于2023年，星期日例1用计算机计算存款模型的各阶差分.(*首先计算20年内的存款清单*)r=0.07;a[0]=1000;a[n_]:=(1+r)a[n-1];money1=Table[{n,a[n]},{n,0,20}];TableForm[Join[{{年份,存款额}},money1]](*其次计算各阶差分*)da[n_]:=a[n+1]-a[n];d2a[n_]:=da[n+1]-da[n];d3a[n_]:=d2a[n+1]-d2a[n];diff=Table[{n,a[n],da[n],d2a[n],d3a[n]},{n,0,9}];TableForm[Join[{{N,An,Dan,D2an,D3an}},diff]]dif1=Transpose[diff];TableForm[{dif1[[3]]/dif1[[2]],dif1[[4]]/dif1[[3]],dif1[[5]]/dif1[[4]]}]第七页，共二十九页，编辑于2023年，星期日差分方程an=b的解由an+1-an=b,n=0,1,2,…,得an-a0=nb.

如果a0=c,则有an=nb+c.

一般地,差分方程kan=b的解是:an=cknk+ck-1nk-1+……+c1n+c0,其中ck=b/k!.验证如下：a[n_]:=c[4]n^4+c[3]n^3+c[2]n^2+c[1]n+c[0];da[n_]:=a[n+1]-a[n];d2a[n_]:=da[n+1]-da[n];d3a[n_]:=d2a[n+1]-d2a[n];d4a[n_]:=d3a[n+1]-d3a[n];d3a[n]//Simplifyd4a[n]//Simplify第八页，共二十九页，编辑于2023年，星期日差分方程an+1=an+b的解

定理1

一阶线性差分方程an+1=an+b的通解是：

定理2

对一阶线性差分方程an+1=an+b,

若||<1,则an无限趋近于平衡解b/(1-)(收敛型不动点);

若||>1,则an逐渐远离平衡解b/(1-)(发散型不动点).第九页，共二十九页，编辑于2023年，星期日家庭教育基金模型由a0=x,an+1=(1+r)an+x,n=0,1,2,3,…

得通解:将a0=x,=1+r,b=x代入,得c=x(1+r)/r,因此方程的特解是:将a18=30000,r=0.1代入计算出x=586.41.第十页，共二十九页，编辑于2023年，星期日购房抵押贷款模型由a0=60000,an+1=(1+r)an-x,n=0,1,2,3,…将=1+r,b=-x代入得到方程的特解:若在第N个月还清贷款，令aN=0,得:将a0=60000,r=0.01,N=25*12=300代入计算出x=631.93.第十一页，共二十九页，编辑于2023年，星期日分期付款模型若小王采取分期付款方式，每月要付300元.如果采用贷款方式，类似于上一模型，将a0=8000,r=0.15/12,N=36代入计算出x=277.32.比较两种支付方式，他应该选择消费贷款方式。第十二页，共二十九页，编辑于2023年，星期日4.3Fibonacci数列问题

13世纪意大利著名数学家Fibonacci在他的著作《算盘书》中记载着这样一个有趣的问题：一对刚出生的幼兔经过一个月可长成成兔，成兔再经过一个月后可以繁殖出一对幼兔.若不计兔子的死亡数，问一年之后共有多少对兔子？月份

01234567…幼兔

10112358…成兔

011235813…总数

1123581321…第十三页，共二十九页，编辑于2023年，星期日将兔群总数记为fn,n=0,1,2,…，经过观察可以发现，数列{fn}满足下列递推关系：

f0=f1=1,fn+2=fn+1+fn,n=0,1,2,…这个数列称为Fibonacci数列.Fibonacci数列是一个十分有趣的数列，在自然科学和数学领域中都有着广泛的应用.

Fibonacci数列的一些实例.1.蜜蜂的家谱

2.钢琴音阶的排列

3.树的分枝

4.杨辉三角形第十四页，共二十九页，编辑于2023年，星期日观察Fibonacci数列(*计算Fibonacci数列的前20项，并作图*)F[0]=F[1]=1;F[n_]:=F[n-1]+F[n-2];fib=Table[F[i],{i,0,20}]tu1=ListPlot[fib,PlotStyle->PointSize[0.018]];(*取对数后再观察，可以发现图像近似一条直线.*)lgf=Log[fib];tu2=ListPlot[lgf,PlotStyle->PointSize[0.018]];(*使用线性函数对数据进行拟合*)f[x_]=Fit[lgf,{1,x},x]tu3=Plot[f[x],{x,0,21},PlotStyle->RGBColor[0,0,1]];Show[tu3,tu2]通过计算可知，fn0.465577e0.478438n.第十五页，共二十九页，编辑于2023年，星期日Fibonacci数列的通项公式

Fibonacci数列满足递推关系fn+2=fn+1+fn，称为二阶线性差分方程.通过前面的计算，可以猜测fn具有指数形式.

不妨设fn=n,代入差分方程，得2--1=0.其解记为1,2.得到差分方程的通解为:

fn=C11n+C22n.

r=Solve[x^2-x-1==0,x];a=x/.r[[1]];b=x/.r[[2]];F1[n_]:=c1a^n+c2b^n;cc=Solve[{F1[0]==1,F1[1]==1},{c1,c2}]//SimplifyF1[n]/.cc[[1]]//Simplify第十六页，共二十九页，编辑于2023年，星期日生成函数对给定数列a0,a1,…,an,…，以{an}为系数构造一个形式幂级数：

G(x)=a0+a1x+a2x2+…+anxn+…称为数列{an}的生成函数(也称为母函数).

例1.有限数列的生成函数是G(x)=(1+x)n.

例2.无穷数列的生成函数是G(x)=ex.

例3.以G(x)=为生成函数的数列是an=2n-1.第十七页，共二十九页，编辑于2023年，星期日Fibonacci数列的生成函数设Fibonacci数列的生成函数是:F(x)=f0+f1x+f2x2+…+fnxn+…,其中fn+2=fn+1+fn.由，得.从而得.

再由f0=f1=1,得:

第十八页，共二十九页，编辑于2023年，星期日4.5分叉与混沌Logistic方程在受环境制约的情况下，生物种群的增长变化行为比较复杂.例如在池塘内，环境可供1000条鱼生存.在鱼的数量远远低于此数时，鱼群的增长接近于指数增长.但当鱼的数量接近生存限时，由于生态环境逐渐恶化，鱼群的增长逐渐变慢，几乎停止增长.如果鱼群数量超过了生存限，由于环境不堪重负，鱼群会出现负增长.这种现象可以用logistic方程进行刻画.pn+1-pn=kpn(N-pn)第十九页，共二十九页，编辑于2023年，星期日

例1

池塘中鱼的数量满足差分方程

pn+1-pn=0.001pn(1000-pn)选择不同的初值，观察鱼群数量的变化趋势.

p[x_]:=2x-0.001x^2;pict[a_]:=Module[{data1},data1=NestList[p,a,30];ListPlot[data1,PlotStyle->{RGBColor[0,0,1],PointSize[0.018]}]]pict[0];pict[1];pict[500];pict[1000];pict[1500];第二十页，共二十九页，编辑于2023年，星期日

例2

学校有两名同学在星期一返校时患了流感，假设流感的传染率为0.002，问两周之后全校400名学生中会有多少人感染过流感？记an

为到第n天时感染过流感的学生人数.假定流感患者的增加速度与流感患者同尚未感染流感的接触次数an(400-an)成正比.因此，an满足logistic方程

an+1-an=0.002an(400-an)

p1[x_]:=x+0.002x(400-x);pict1[a_]:=Module[{data1},data1=NestList[p1,a,14];ListPlot[data1,PlotStyle->{RGBColor[0,0,1],PointSize[0.018]}]]pict1[2];第二十一页，共二十九页，编辑于2023年，星期日Logistic方程的迭代

logistic方程是非线性方程，其标准形式为:

an+1=ran(1-an),下面通过实验观察迭代数列的收敛性.

logistic[r_,a_,n_]:=Module[{p,data,tu1,tu2},p[x_]:=rx(1-x);data=NestList[p,a,n];tu1=ListPlot[data,PlotStyle->PointSize[0.018],DisplayFunction->Identity];tu2=ListPlot[data,PlotJoined->True,PlotStyle->RGBColor[0,0,1],DisplayFunction->Identity];Show[tu1,tu2,DisplayFunction->$DisplayFunction]];第二十二页，共二十九页，编辑于2023年，星期日

(*初值r=0.7,a0=0.2*)logistic[0.7,0.2,30];

容易看出，迭代数列单调收敛于0.

(*初值r=2.9,a0=0.2*)logistic[2.9,0.2,30];

迭代数列上下振荡，趋向于不动点(r-1)/r.

(*初值r=3.4,a0=0.2*)logistic[3.4,0.2,30];

经过一段时间的调整，迭代数列开始接近在0.42和0.82之间振荡.这类振荡称为2-循环.第二十三页，共二十九页，编辑于2023年，星期日

(*初值r=3.55,a0=0.2*)logistic[3.55,0.2,30];

出现了周期为4的振荡，称为4-循环.通过以上的观察可以发现，当参数r变化时，相应的迭代数列从收敛到唯一的不动点(1-循环)到2-循环再到4-循环，这样的分裂行为称为分叉(bifurcation).

(*初值r=3.7,a0=0.2*)logistic[3.7,0.2,30];

此时没有稳定的周期性.迭代数列在区间(0,1)内振荡，而且表现出对初始条件非常敏感的依赖性，这种状态称为混沌(chaos).第二十四页，共二十九页，编辑于2023年，星期日Feigenbaum图设f(x)是定义在实数域上的实值函数，如果存在x*，使得f(x*)=x*，则称x*为f(x)不动点.如果所有附近的点在迭代过程中都趋于某个不动点，则称该不动点为吸引点，或称为稳定点；如果所有附近的点在迭代过程中都远离它而去，则称该项点为排斥点(不稳定点).如果f(a1)=a2,f(a2)=a3,…,f(ak)=a1,并且aja1,j=2,3,…,k,则a1,a2,…,ak构成一个k-循环.a1称为k-周期点，a1,a2,…,ak称为一个k-周期轨道.为了观察r

对迭代格式an+1=r

an(1-an)的影响，将区间(0,4]以步长r离散化.对每个离散的r

值进行迭代，忽略前50个迭代值，把点(r,a51),(r,a52),…,(r,a100)显示在坐标平面上.这样形成的图形称为Feigenbaum图,它反映了混沌与分叉的基本特性.第二十五页，共二十九页，编辑于2023年，星期日

Feig[n_,x0_]:=Module[{plist={},a,i,temp,pilist={}},For[a=1,a<=n,a++,temp=x0;plist={};For[i=1,i<=50,i++,temp=4*a*temp(1-temp)/n];For[i=51,i<=100,i++,temp=4*a*temp(1-temp)/n;AppendTo[plist,{4a/n,temp}]];AppendTo[pilist,ListPlot[plist,PlotStyle->{RGBColor[1,0,0],PointSize[0.008]},DisplayFunction->Identity]]];Show[pilist,DisplayFunction->$DisplayFunction]];Feig[500,0.2];第二十六页，共二十九页，编辑于2023年，星期日练习二对迭代格式

an+1=4an(1-an),n=1,2,…使用初值0.21进行迭代，记录前100次的迭代数据.把[0,1]区间十等分，统计迭代数列中落在各个小区间内的项数，做出统计表.迭代数列在[0,1]区间内分布均匀吗？任取区间(0,1)内的一些初值重复这一实验，总结实验结果.对迭代格式

an+1=3.45an(1-an),n=1,2,…重复上述实验，实验结果有区别吗？你能解释这个现象吗？第二十七页，共二十九页，编辑于2023年，