可以写成某个参数s的函数。另外因为F(u,K)=-u3+Ku连续,Fu在除原点之外的各点都非奇异,所以(3)式的解在除原点之外的各点附近存在且唯一,并且解曲线局部上可以用K作为自然参数表示,这与我们已经求出了3个解u(i)(K)(i=1,2,3)是一致的。y(i)(K)={(u(i)(K),K)|KIR1}(i=1,2,3)是(3)式的解曲线,并且在y(1)(0)=(0,0)的任何E方形邻域内都存在不落在y(1)(K)上的(1)的解点,因为(-E/2,E2/4)Iy(2)(K)和(-E/2,E/2)Iy(2)(K)必有其一在这个邻域内。根据定义2,原点是(3)式的分岔点。因为原点是3条解曲线的交点,那么根据定义3,原点不是(3)式的简单分岔点。因为原点两侧都有稳定的解点,所以根据定义4,原点处的分岔是超临界叉形分岔。313分岔点的判别准则分岔点的定义无疑是一条判别准则,但却并非每个非线性特征值问题(1)式都会简单到像(3)式那样能方便地应用定义判别分岔点。实际计算中常用到如下一条比较实用的判别准则:当p=1时,如果y(s)是(1)式的光滑解曲线且行列式d(s)=detFc(y(s))Ûy(s)T在s=s*处变号,则y(s*)必是(1)式的分岔点。但这仅仅是/y(s*)是分岔点0的充分条件,单纯应用此准则是可能遗漏分岔点的!比如对(3)式,(0,0)就不能通过这条准则被确认为分岔点。首先,3条解曲线中y(1)(K)和y(2)(K)都在原点处不光滑(导数不存在);其次,即使是针对在原点处光滑的解曲线y(3),d(s)也在原点处失去意义。我们有理由认为,在判别分岔点时,应充分注意一些可疑点,比如:解曲线上罕见的不光滑点、恰好使某些判别准则失效的点等等,要对可疑点附近的性质进行数值分析,进而判别它是否为分岔点。4非线性特征值问题的正则解及步长控制411正则解的定义及两种算法对于(1)式的解曲线上某一点y(s0),若F0y=Fy(y(s0))的秩为n,则y(s0)称为(1)式的正则解。再考虑方程(3)。因为n=1,Fy(u,K)=(-3u2+K,u),所以解曲线上除原点之外的每一点都是(3)式的正则解。方程(1)的解是连续曲线而不是孤立的点,数值计算中就用曲线上适当分布的点来刻画该曲线。41111预估-校正算法假设(u0,K0)是方程(1)的解且Fu(u0,K0)非奇异,由隐函数定理知在K0的一个邻域内,(1)式存在一条唯一的解曲线y(K)。计算这类正则解,预估-校正算法是最简单直观的。预估就是通过(u0,K0)的信息构造出y(K)上(u0,K0)附近一点(u1,K1)的近似值(u1,K1),$K=K1-K0称为延拓步长;校正则指对固定的K1,以u1为初始值构造出收敛于y(K)上的点(u1,K1)的迭代过程。41112连续法假设(u0,K0)是方程(1)的解,当Fu(u0,K0)奇异时,预估-校正算法将无法进行。而基于隐函数定理的连续法却可以计算所有使(4)式成立的正则点。连续法的基本思想是引入其它参数,以避免因为采用自然参数K遇到Fu(u,K)的奇异点。412步长控制求解(1)式的数值方法若对一个实际问题而言是可行的,那么计算其解曲线上的点越少,计算的代价就越低。因此衡量算法优劣的标准之一就是延拓步长$K能否取得尽可能大,但步长太大就会影响收敛的速度。适当的步长与校正迭代的最佳步数、校正算法、要求的精度等有关。5数值计算实例即使我们未求出(3)式的解,在4.1中我们也已经知道(3)式的0绝大部分0解是正则解,所以我们可以采用预估-校正算法来求(3)式的解曲线。以(ui,Ki)标记求解过程中出现的第i个解点。511割线法预估-校正第一步找到(3)式的两个正则解(u0,K0)、(u1,K1),作为本算法的预估起点;第二步根据(ui,Ki)和(ui+1,Ki+1),以$Ki为步长进行延拓,进行割线预估,找到新解点(ui+2,Ki+2)的近似点(ui+2,Ki+2):Ki+2=Ki+1+$Ki,(i=0,1,2,,);ui+2=ui+1+(ui+1-ui)#Ki+2-Ki+1Ki+1-Ki=ui+1+$Ki#ui+1-uiKi+1-Ki,(i=0,1,2,,);95第2期姚俐等:动态人口模型分岔问题与数值计算u(t)=au0exp[a(t-t0)]a-bu0+bu0exp[a(t-t0)],其中a,b为正的常数,limty]u(t)=ab。本文将主要讨论K=B1-uu12情形下的广义人口模型:dudt=B1-uu12u,其中B0为常数,u1表示该地区人口最大承载量。为论述简明起见,令K=u21,S=t#Bu21,易知dudS=-u3+Ku¦F(u,K)(Ó)这是一个变量可分离的微分方程,容易求出其解为:u(t)=u0K[u20+(K-u20)exp(2K(t0-t))]u20+(K-u20)exp(2K(t0-t)),K的实际意义是该地区人口的最大承载量的平方。下面针对模型(Ó)的平衡解给出其稳定性分析。事实上,该地区人口相对稳定状态对应微分方程(Ó)的平衡解(steadystates),即下述方程之解:F(u,K)=-u3+Ku=0(3)显然,(3)式是一个简单的含参数非线性特征值问题。在不考虑K的实际意义时,认为KIR1,分情况讨论如下。情形1K[0此时uS0是(3)式的唯一解,且dudS0,u0,0,u0。情形2K0此时uS0和u=?K是(3)式的3个解,且dudS0,u-K或0uK,0,-Ku0或uK。因此,(Ó)的3个平衡解就是:u(1)(K)=0,K[0,K,K0;u(2)(K)=0,K[0,-K,K0;u(3)(K)S0。它们的稳定性如图1所示。图中箭头表示各解在相应区域的变化趋势。可以看出,非零解都是稳定的平衡解。图1解的流形3分岔现象通过前一节的论述,可以直观地看出:在某些相对稳定的假设条件下,人口模型的稳态解可能不唯一,并且可能在一点处相交。因此,有必要对这种现象加以讨论。隐函数定理告诉我们,只要函数F定义4如果分岔点y*的左右两侧都有稳定的解点,那么y*处的分岔称为超临界叉形分岔(supercriticalpitchforkbifurcation)。以上基本概念,在文献[3,4]中都有论述。312实际模型的分岔现象就(3)式而言,n=p=1,Fu(u,K)=-3u2+K,FK(u,K)=u,那么除原点之外的任意点都满足(4)式。因此除原点之外各点附近方程(3)的解总94信息工程大学学报2003年收稿日期:2003-01-15基金项目:河南省杰出青年科学基金资助课题(0016);河南省高校杰出科研人才创新工程项目。作者简介:姚俐(1980-),女,四川开江人,信息工程大学硕士研究生,主要研究方向为数值分析与应用。动态人口模型分岔问题与数值计算姚俐1,江成顺1,刘广彦2(11信息工程大学信息工程学院,河南郑州450002;21信息工程大学理学院,河南郑州450001)摘要:考虑动态人口模型。主要探讨其平衡解的稳定性以及在适当条件下解的分岔问题。研究非线性特征值问题的正则解的一般控制计算方法。针对具体模型,给出了一系列的计算实例,最后给出了研究结论。关键词:人口模型;平衡解;分岔;特征值问题中图分类号:O177191文献标识码:A文章编号:1671-0673(2003)02-0093-04BifurcationforSomeDynamicPopulationModelsandNumericalComputationYAOLi,JIANGCheng-shun,LIUGuang-yan(InstituteofInformationEngineering,InformationEngineeringUniversity,Zhengzhou450002,China)Abstract:Thispaperisconcernedwithsomedynamicpopulationmodels.Theauthorsdiscussthestabilityofsteadystatesandthebifurcationproblem.Somecontrollingandnumericalcomputationmethodsforsomenonlineareigenvalueproblemsarestudied.Aimedatthegivenmodel,aseriesofcomputationexperimentsarepresentedandsomeconclusionsmade.Keywords:populationmodel;steadystate;bifurcation;eigenvalueproblem1引言动态人口系统数学模型的参数反演计算问题往往归结为求解非线性问题,FBDRn@RpyRn,F(u,K)=0,uIRn,KIRp(1)即所谓的含参数K的非线性特征值问题。这类问题的分析与计算,一般都相当复杂。本文针对Malthus人口模型的推广模型,论述其解的定性性质,着重讨论模型的分岔问题和数值方法。2基本模型及其稳定性分析首先考虑最简单的动态人口模型。设u(t)表示某地区在t时刻的人口总数,u0是t0时刻的人口总数,可以建立如下模型:dudt=F(u,K)u(t0)=u0(2)其中KIRp为与t无关的模型参数。后文中各模型的初值条件也都默认为u(t0)=u0。当F(u,K)=Ku并且K为一正的常数时,(2)式即著名的Malthus人口模型[1~3]:dudt=Ku(Ñ)其解为u(t)=u0exp[K(t-t0)],K的实际意义是该地区人口的净增长率。当F(u,K)=-u2+Ku并且K为一正的常数时,(2)式可演化为著名的Logistic人口模型[1~3]:dudt=(a-bu)u(Ò)其解为第4卷第2期信息工程大学学报Vol.4No.22003年6月JournalofInformationEngineeringUniversityJun.2003第三步对于固定的Ki+2,以ui+2为初始值构造出收敛于y(K)上的点(ui+2,Ki+2)。如算法所述,延拓的步长是应该实时调整的,而简单起见,我们现在假定延拓的步长恒定。此算法可以用MATLAB编程实现,下面的实例中,每一步的结果都精确到小数点后4位。实例对于起点(u0,K0)=(0.1,0101)和(u1,K1)=(011414,0102),分别采用$K=011,10各延拓10次,结果如表1所示。表1割线法实例i$K=011$K=10KiuiuiKiuiui0010100/011000010100/0110001010200/011414010200/011414201120001596801346410102004115828-31165430.22000.59240.46902010200-614789-41474440.32000.73880.56573010200-710949-51479150.42000.87490.64814010200-814951-61326160.52000.99410.72115010200-917162-71072570.62001.10090.78746010200-1018057-71747380.72001.19840.84857010200-1117973-81367890.82001.28860.90558010200-1217128-819454100.92001.37300.95929010200-1315673-914879111.02001.45251.010010010200-1413714-1010010512切线法预估-校正第一步找到(3)式的一个正则解(u0,K0),作为本算法的预估起点;第二步根据(ui,Ki),以$Ki为步长进行延拓,进行切线预估,找到新解点(ui+1,Ki+1)的近似点(ui+1,Ki+1):Ki+1=Ki+$Ki,(i=0,1,2,,);ui+1=ui