清华数学实验2-差分方程和数值微分

数学实验 李骥聪 2007011861 实验实验 2 2 差分方程和数值微分差分方程和数值微分 化学工程系 化 71 李骥聪 2007011861 【实验目的实验目的】 1 掌握用 MA TLAB优化工具箱的基本用法，对不同算法进行初步分析、比较。 2 练习用无约束优化方法建立和求解实际问题模型（包括非线性最小二乘拟合） 。 【实验内容实验内容】 1.1.题目题目 2.2. 一老人 60 岁时将养老金 10 万元存入基金会，月利率 0.4%，他每月取 1000 元作为生活 费，建立差分方程计算他每岁末尚有多少钱？多少岁时将基金用完？如果想用到 80 岁，问 60 岁时应存入多少钱？ 模型及其求解模型及其求解 假定基金利息为复利计算，设老人 60 岁时存入 0 x元，月利率为r， 每月月末取款为b，第k月月末的余额为 k x。基于以上假定，则第k月末的余额 k x只和第 1k 月末的余额 1k x 及月利息r有关，可以得到如下关系式： 1 (1) kk xr xb 记(1)ar，则： 1kk xaxb 代入题目中给定的数据，在 MA TLAB中计算。源程序为： function x=exf1(x0,m,r,b) a=1+r; x=x0; for k=1:m x(k+1)=a*x(k)-b; end k=(0:240); x1=exf1(100000,240,0.004,1000); y=1:20; M(y+1)=x1(12*y); M(1)=100000; figure, plot(0:20,M,r*,0:20,M,-);%每年年末余额变化(x0=100000) title(图 1.每年年末余额) xlabel(年份) ylabel(余额) gtext(x0=100000); hold on, z=0*(0:20); plot(0:20,z),hold off r=0.004; a=1+r; 数学实验 李骥聪 2007011861 b=1000; X=b*(1-a240)/a240/(1-a) Y=b/r x2=exf1(X,240,0.004,1000); x3=exf1(200000,240,0.004,1000); x4=exf1(Y,240,0.004,1000); x5=exf1(300000,240,0.004,1000); figure, plot(k,x1,r,k,x2,b,k,x3,g,k,x4,k,x5, m) title(图 2.最初存入金额不同的余额变化曲线（按月）); xlabel(月) ylabel(月末余额) gtext(x1=100000); gtext(x2=X); gtext(x3=200000); gtext(x4=Y); gtext(x5=300000); （1 1）由 MA TLAB计算得到 1 至 20 年的年末余额，如下表： 表表 1 1. .存入基金后每年年末余额存入基金后每年年末余额 年份年份 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 1010 余额余额 93266 85575 77507 69043 60163 50848 41075 30824 20068 8785.7 年份年份 1111 1212 1313 1414 1515 1616 1717 1818 1919 2020 余额余额 -3050.8 -15468 -28495 -42160 -56497 -71537 -87314 -1.04E+05 -1.21E+05 -1.39E+05 表表 2.第第 10 年每月月末余额年每月月末余额 月份月份 1 1 2 2 3 3 4 4 5 5 6 6 余额余额 7820.8 6852.1 5879.5 4903 3922.6 2938.3 月份月份 7 7 8 8 9 9 1010 1111 1212 金额金额 1950.1 957.89 -38.278 -1038.4 -2042.6 -3050.8 从表 2 得出在存入后第 10 年的 8 月，即老人 70 岁零 8 个月时，基金将用完。 MA TLAB 绘制出曲线如图 1： 数学实验 李骥聪 2007011861 （2 2）由差分方程 1kk xaxb 逐步递推得到： 21 00 1 (1) 1 k kkk k a xa xbaaaa xb a 按照题目要求，若在老人 80 岁时仍有余额，即 240 0 x，所以 240 0 240 1 (1) a xb aa 由 MA TLAB 计算得到，第 20 年年末余额刚好为零时，初始存入金额X为 154093.3 元。 结果分析结果分析 在建立本题目的模型过程中，首先作了一些基本假定： 1. 基金利息以复利计算，并且利率不变。 2. 老人每月月末取用下月的生活费，并且每月取用金额相同。 在此基础上，建立了形如： 1kk xaxb 的差分方程数学模型，并对问题求解。 利用 MA TLAB画出了老人在 60 岁时存入不同金额时， 每月月末余额 （并用曲线连接） ， 如图 2： 数学实验 李骥聪 2007011861 （1 1）从图上曲线的变化趋势可以很清楚的看到，在r和b相同的条件下，初始存入的 基金金额 0 x是影响每月月末余额的重要因素。 当 0 x小于X时， 余额在第 20 年年末为负值； 当 0 x等于X时， 余额在第 20 年年末为零； 当 0 x大于X时， 余额在第 20 年年末将为正值。 （2 2）从图上曲线也可以看出：当 0 x小于Y时，余额数将减小，这是因为每月产生的利 息小于月末的取款；当 0 x等于Y时，值数将不变；当 0 x大于某值Y时，因为每月产生的利 息大于月末的取款，老人的基金将累积，这时 0 x满足： 0 x rb 不等式取等号时 b Y r ，由 MA TLAB解得初始金额Y为 250000 元。 由 0 1 (1) 1 k kk k ab xa xbcr ar ，其中 0 () b cx r ，可知， b Y r 是差分方程的 平衡点，由于11r ，故平衡点 b Y r 是不稳定的。这与图中曲线在k 发散是相符 合的。 （3 3）在现实中，基金利率r和老人每月取用生活费的时间和金额b都会随时间发生变 化，需要在差分方程中考虑对利率r和金额b的修正，若取用时间不在月初月末，那么还要 考虑基金的日利率 r的情况。方程将变得很复杂。 数学实验 李骥聪 2007011861 2.2.题目题目 8 8 在某种环境下猫头鹰的主要食物来源是田鼠， 设田鼠的年平均增长率为 1 r， 猫头鹰的存 在引起的田鼠增长率的减少与猫头鹰的数量成正比，比例系数为 1 a；猫头鹰的年平均减少 率为 2 r；田鼠的存在引起的猫头鹰减少率的增加与田鼠的数量成正比，比例系数为 2 a，建 立差分方程模型描述田鼠和猫头鹰共处时的数量变化规律，对以下情况作图，给出 50 年的 变化过程。 （1） 设 1212 0.2,0.3,0.001,0.002,rraa开始时有100只田鼠和50只猫头鹰。 （2） 1 r 2 r 1 a 2 a同上，开始时有 100 只田鼠和 200 头猫头鹰。 （3） 适当改变参数 1 a 2 a（初始值同上） 。 （4） 求差分方程的平衡点，它们稳定吗？ 模型及其求解模型及其求解 如题假设，若田鼠和猫头鹰的初始数量分别设为 0 x， 0 y。第k年年 末田鼠和猫头鹰的数量分别为 k x， k y。按照题目中的假定，可以写出如下关于 k x， k y的 方程组： 111 (1) kkk xra yx 122 (1) kkk yra xy 代入题目中给定的数据，在 MA TLAB中计算。源程序为： function z=exf2(x0,y0,r1,r2,a1,a2) x=x0; y=y0 for k=1:50 x(k+1)=(1+r1-a1*y(k)*x(k); y(k+1)=(1-r2+a2*x(k)*y(k); end z=x,y; z1=exf2(100,50,0.2,0.3,0.001,0.002); k=0:50; plot(k,z1(:,1),r,k,z1(:,2),b); title(图 1.(1)条件下 50 年内的变化); xlabel(年) ylabel(数量) gtext(猫头鹰); gtext(田鼠); z2=exf2(100,200,0.2,0.3,0.001,0.002); k=0:50; plot(k,z2(:,1),r,k,z2(:,2),b); title(图 2.(2)条件下 50 年内的变化); xlabel(年) ylabel(数量) gtext(猫头鹰); gtext(田鼠); z3=exf2(100,200,0.2,0.3,0.001,0.0025); z4=exf2(100,200,0.2,0.3,0.0015,0.002); 数学实验 李骥聪 2007011861 z5=exf2(100,200,0.2,0.3,0.001,0.0015); z6=exf2(100,200,0.2,0.3,0.0005,0.002); k=0:50; subplot(2,2,1), plot(k,z3(:,1),r,k,z3(:,2),b); title(3)a1=0.001,a2=0.0025 条件下 50 年内的变化); xlabel(年) ylabel(数量) gtext(猫头鹰); gtext(田鼠); subplot(2,2,2), plot(k,z4(:,1),r,k,z4(:,2),b); title(3)a1=0.0015,a2=0.002 条件下 50 年内的变化); xlabel(年) ylabel(数量) gtext(猫头鹰); gtext(田鼠); subplot(2,2,3), plot(k,z5(:,1),r,k,z5(:,2),b); title(3)a1=0.001,a2=0.0015 条件下 50 年内的变化); xlabel(年) ylabel(数量) gtext(猫头鹰); gtext(田鼠); subplot(2,2,4), plot(k,z6(:,1),r,k,z6(:,2),b); title(3)a1=0.0005,a2=0.002 条件下 50 年内的变化); xlabel(年) ylabel(数量) gtext(猫头鹰); gtext(田鼠); z7=exf2(100,200,0.2,0.3,0.001,0.0025); z8=exf2(100,200,0.2,0.3,0.001,0.0027); zph=exf2(100,200,0.2,0.3,0.001,0.003); z81=exf2(100,200,0.2,0.3,0.001,0.0032); k=0:50; plot(k,z7(:,1),r,k,z7(:,2),b,k,z8(:,1) ,r-.,k,z8(:,2),b-.) hold on, k=0:50; plot( k,zph(:,1),r:,k,zph(:,2),b:, k,z81(:,1),r-,k,z81(:,2),b-), hold off title(图 4.(3)平衡点趋于初始值时，50 年内的变化); xlabel(年) ylabel(数量) gtext(猫头鹰 a1=0.001,a2=0.0025); gtext(田鼠 a1=0.001,a2=0.0025); 数学实验 李骥聪 2007011861 gtext(a1=0.001,a2=0.0027); gtext(a1=0.001,a2=0.0027); gtext(a1=0.001,a2=0.003); gtext(a1=0.001,a2=0.003); gtext(a1=0.001,a2=0.0032); gtext(a1=0.001,a2=0.0032); z9=exf2(150,200,0.2,0.3,0.001,0.002); z10=exf2(151,200,0.2,0.3,0.001,0.002); z11=exf2(150,201,0.2,0.3,0.001,0.002); k=0:50; plot(k,z9(:,1),r,k,z9(:,2),b,k,z10(:,1) ,r-.,k,z10(:,2),b-. ,k,z11( :,1),r:,k,z11(:,2),b:) title(图 5.初始值在平衡点附近，50 年内的变化); xlabel(年) ylabel(数量) gtext(猫头鹰 200); gtext(田鼠 150); gtext(猫头鹰 201); gtext(田鼠 151); function z=exf3(x0,y0,r1,r2,a1,a2) x=x0; y=y0; for k=1:100 x(k+1)=(1+r1-a1*y(k)*x(k); y(k+1)=(1-r2+a2*x(k)*y(k); end z=x,y; z12=exf3(150,200,0.2,0.3,0.001,0.002); z13=exf3(151,200,0.2,0.3,0.001,0.002); z14=exf3(150,201,0.2,0.3,0.001,0.002); k=0:100; plot(k,z12(:,1),r,k,z12(:,2),b,k,z13(:,1) ,r-.,k,z13(:,2),b-.) hold on, k=0:100; plot( k,z14(:,1),r:,k,z14(:,2),b:), hold off title(图 6.初始值在平衡点附近，100 年内的变化); xlabel(年) ylabel(数量) gtext(猫头鹰 200); gtext(田鼠 150); gtext(猫头鹰 201); gtext(田鼠 151); （1 1）在 1212 0.2,0.3,0.001,0.002,rraa开始时有 100 只田鼠和 50 只猫头鹰的 条件下，由 MA TLAB绘图得到田鼠和猫头鹰共处时 50 年中的数量的变化过程，如图 1： 数学实验 李骥聪 2007011861 图中，田鼠和猫头鹰的数量交替增加、减少。将田鼠和猫头鹰将田鼠和猫头鹰各自各自每相邻两次数量达每相邻两次数量达 到最小值之间的时间段叫做一个增长周期。到最小值之间的时间段叫做一个增长周期。 从图中可以看到刚开始猫头鹰的数量减少，田鼠的数量增加。在前 25 年里，当田鼠增当田鼠增 长最快时，猫头鹰的数量达到最低。 当田鼠数量达到最大时， 猫头鹰的增长速度达到最快。长最快时，猫头鹰的数量达到最低。 当田鼠数量达到最大时， 猫头鹰的增长速度达到最快。 这 一 点 可 由 方 程 122 (1) kkk yra xy 的 性 质 分 析 得 出 ： 猫 头 鹰 年 增 长 的 数 量 12222 (1)() kkkkkkk yyyra xyyra xy ，田鼠的数量 k x越大，猫头鹰的增 长速度就越快。 图中猫头鹰的数量曲线有较尖锐的峰，这是由于图中猫头鹰的数量曲线有较尖锐的峰，这是由于 2121 ,aa rr，猫头鹰的数量变化的，猫头鹰的数量变化的 幅度大于田鼠。幅度大于田鼠。 图中第 20 至第 25 年， 田鼠的数量维持在很低的水平， 这显然是不符合的现实繁殖规律 的，因为田鼠种群可持续繁殖的种群数量存在最小极限值，当种群数量低于这个值时，在实 际情况下，田鼠的种群就要灭绝，从而导致猫头鹰的数量不断减少，最终灭绝。 （2 2）在 1212 0.2,0.3,0.001,0.002,rraa开始时有 100 只田鼠和 200 只猫头鹰 的条件下，由 MA TLAB绘图得到田鼠和猫头鹰共处时 50 年中的数量的变化过程，如图 2： 数学实验 李骥聪 2007011861 从图中可以看出，刚开始时猫头鹰的数量是（1）条件时的四倍，这大大减弱了田鼠在 第一个增长周期时的数量峰值，从而也抑制了猫头鹰在第一个增长周期内的数量峰值。 由 MA TLAB 计算得到第 35 年时田鼠的数量为 110 只， 猫头鹰的数量为 50 只， 这与 （1） 的初始条件近似，所以（2）条件下第 35 至 50 年的变化情况与（1）条件下前 15 年的变化 情况类似。 田鼠和猫头鹰的数量随时间波动， 图中波动的中心值不明显 （即平衡点不能直接读出） ，田鼠和猫头鹰的数量随时间波动， 图中波动的中心值不明显 （即平衡点不能直接读出） ， 波动范围越来越大（即田鼠和猫头鹰数量不收敛） 。波动范围越来越大（即田鼠和猫头鹰数量不收敛） 。 （3 3）在 12 0.2,0.3,rr开始时有 100 只田鼠和 200 只猫头鹰的条件下，对 12 ,a a作 适当改变，在保持其中一个不变的情况下，另一个参数增加或减少 0.0005，由 MA TLAB绘 图得到田鼠和猫头鹰共处时 50 年中的数量的变化过程，如图 3： 从第一幅图中可以看出 a1=0.001,a2=0.0025 时， 在第一增长周期中猫头鹰的数量在 200 附近波动， 田鼠的数量在 120 附近波动， 波动的波动的平衡点较明显， 这是因为初始值 （平衡点较明显， 这是因为初始值 （100100， 200200） 与图中读出的平衡点值（与图中读出的平衡点值（120120，200200）接近，所以能在短时间内保持一定的稳定性，随着时）接近，所以能在短时间内保持一定的稳定性，随着时 间推移，波动就变得很剧烈了。间推移，波动就变得很剧烈了。 在第三幅图中，猫头鹰的数量在 300 附近波动，由于于初始值 200 相差较大，波动在第 二增长周期就已经相当剧烈，到第 40 年时，猫头鹰几乎已经没绝，田鼠疯狂增长。 第四幅图，在 a1=0.0005,a2=0.002 时，由 MA TLAB计算可以达到，在第 48 年时，田 鼠和猫头鹰的数量为 268 和 2435，第 49 年是变为了-5 和 3008。这显然是不符合实际情况 的，即在第即在第 49 年时，田鼠已经灭绝，其后的图像不再具有实际意义，而应该是猫头鹰数量年时，田鼠已经灭绝，其后的图像不再具有实际意义，而应该是猫头鹰数量 逐渐减少为零。逐渐减少为零。 在 12 0.2,0.3,rr开始时有 100 只田鼠和 200 只猫头鹰的条件下，由由中得出的，中得出的， 数学实验 李骥聪 2007011861 当初始值与平衡点当初始值与平衡点接近时，波动较稳定，而平衡点接近时，波动较稳定，而平衡点与与 1 r 2 r 1 a 2 a有关有关， （这一点将在结果分， （这一点将在结果分 析中得到证明）析中得到证明） ， 下面通过改变， 下面通过改变 12 ,a a使得平衡点接近初值使得平衡点接近初值。 令 12 ,a a分别为 （0.001， 0.0025） ， （0.001，0.0027） ， （0.001，0.003） ，由 MA TLAB绘图得到田鼠和猫头鹰共处时 50 年中的 数量的变化过程，如图 4： 从图中可以看出，当 12 ,a a为（0.001，0.003）时，平衡点与初始值一样为（100，200） ， 此时田鼠和猫头鹰的数量将保持恒定。当 2 a发生微小变化时，平衡点也发生微小变化，平 衡点将不同于初始值，两物种的数量均开始波动，并且 2 a偏离 0.003 越多，波动的越厉害。 当 2 0.003a 时，田鼠的平衡点值小于 100，即平衡点小于初始值，所以田鼠的数量会 下降接近平衡点。 如图 4 中田鼠 （红线） a1=0.001， a2=0.0032， 开始时下降靠近平衡点 93.75。 当到达平衡点后，数量并不能达到稳定，而是继续减少， “冲过”平衡点，然后再上升靠近当到达平衡点后，数量并不能达到稳定，而是继续减少， “冲过”平衡点，然后再上升靠近 平衡点。如此往复，并且“冲过”平衡点。如此往复，并且“冲过”程度程度越来越大，曲线发散。越来越大，曲线发散。 当 2 0.003a 时，田鼠的平衡点值大于 100，即平衡点大于初始值，所以田鼠的数量会 下降接近平衡点。如图 4 中田鼠（红线）a1=0.001，a2=0.0027 和 a1=0.001，a2=0.0025， 开始时上升靠近平衡点 111.1 和 120。并且 2 a越小，波动范围越大。 （4 4） 分析图 1 至图 4 得出： 以上各种以上各种非平衡非平衡条件下， 田鼠和猫头鹰的个数随时间波动，条件下， 田鼠和猫头鹰的个数随时间波动， 且波动范围越来越大，即且波动范围越来越大，即 k x， k y不收敛。不收敛。在在 1212 0.2,0.3,0.001,0.002,rraa条件条件 下，下，平衡点及其稳定性平衡点及其稳定性为（为（150150，200200） ，平衡点不稳定） ，平衡点不稳定（将在结果分析中给出证明或分析） 。 结果分析结果分析 建立本模型的基本假设： 除题目假设条件外， 田鼠和猫头鹰的数量不受其 它条件的影响， 如田鼠的食物总量对田鼠数量的限制以及猫头鹰的其它食物来源对猫头鹰数 量的增加。由题设建立了如下非线性差分方程组数学模型： 111 (1) kkk xra yx 122 (1) kkk yra xy 数学实验 李骥聪 2007011861 （1 1）令方程组中 * 1kk xxx ， * 1kk yyy ，则方程组化为： * 11 (1)xra y x * 22 (1)yra x y 解此方程组得到该方程组的平衡点： * 2 2 r x a ， * 1 1 r y a ，或 * 0 x ， * 0y （此种情 况不符合实际情况，故不再分析） 。代入数据求出 * (,)(150,200)xy 。 （2）若有 * (,)(,) kkk xyx y ，则 * (,)xy是稳定平衡点。否则不稳定。 为分析题目中平衡点的稳定性，考虑初始值为 * (,)xx y和 * (,)xyy时，田鼠和 猫头鹰数量(,) kk xy的变化。 用 MA TLAB作图得到在 1212 0.2,0.3,0.001,0.002,rraa的条件下，初始值为 （150，200） ， （150，201） ， （151，200）时田鼠和猫头鹰 50 年内的数量变化，如图 5： 从图中可以看出，当初始条件与平衡点有微小偏差（此处仅有 1 个单位的偏差）时， (,) kk xy都会在平衡点附近波动，由于偏差很小，波动的振幅开始时并不大，但是随着时间 的推移，波动振幅越来越大，显示出(,) kkk xy 不收敛于 * (,)xy，所以该平衡点为不稳 定平衡点。 田鼠和猫头鹰数量波动的“动力”是初始数量和平衡点的偏差，偏差越大， “动力”越田鼠和猫头鹰数量波动的“动力”是初始数量和平衡点的偏差，偏差越大， “动力”越 大，并且随时间增加而增大。大，并且随时间增加而增大。在实际情况中，当一种物种的数量小于其在该系统中能够繁在实际情况中，当一种物种的数量小于其在该系统中能够繁 殖的最小值时，该物种将灭绝。如果田鼠先灭绝，那么猫头鹰的数量将逐渐减少，最终灭殖的最小值时，该物种将灭绝。