td?xxb?2xa?2xd??td???2x1xd?1xc??1xd?21源时蓝鲸数目会减少，其减少速度与当时蓝鲸的数目成线性关系，即
dx1??cx1(t) （1） . dt当有食物来源时，蓝鲸的数目会增加。增加的速度和它捕食的数目有关，即
dx1= dx1(t) x2(t) （2） . dt合并（1）和（2），得到蓝鲸变化速度满足的微分方程
dx1??cx1(t)? dx1(t) x2(t). (3) dt同样，在没有蓝鲸时，磷虾的增加速度满足
dx2=ax(t). (4) dt考虑到被捕食情况，则磷虾的数目满足
dx2=ax(t)－bx1(t) x2(t). (5) dt合并（3）和（5），得到著名的Lotka-Volterra方程
（6）
其中a,b,c,d均为正常数。
（6）是一个非线性常微分方程组，不可能有解析解。
实验要求：假设a?1.2,b?0.6,c?0.8,d?0.3,而且初始值为x1(0)=2,
x2(0)=1.
1）
用不同的数值方法求解（6）。把x1(t) 和x2(t)画在同一张
图上并给以解释。
6
2） 把（6）的两个方程相除，得到
dx1?cx1?dx1x2 （7） ?dx2ax2?bx1x2用不同的数值方法解出x2~x1之间的函数关系。并把它画在以x1,x2为坐标的图上，对所得结果加以解释。
附1：部分数值分析常用实验程序
1）Guass列主元消元法解线性方程组：
//高斯列主元消去法求解非线性方程组源程序清单 # define NUMBER 20 # include \
float A[NUMBER][NUMBER+1] ,ark; int flag,n;
exchange(int r,int k) { int i;
for(i=1;i<=n+1;i++) A[0][i]=A[r][i];
for(i=1;i<=n+1;i++) A[r][i]=A[k][i];
for(i=1;i<=n+1;i++) A[k][i]=A[0][i];} float max(int k) { int i;
float temp=0;
for(i=k;i<=n;i++)
if(fabs(A[i][k])>temp) { temp=fabs(A[i][k]); flag=i;} return temp; } main()
{ float x[NUMBER],b[NUMBER]; int r,k,i,j;
me: printf(\用Gauss列主元消元法解线性方程组\ printf(\ scanf(\
printf(\现在输入系数矩阵A和向量b:\\n\\n\
7
for(i=1;i<=n;i++) for(j=1;j<=n+1;j++)
{ if(j==n+1) printf(\请输入b%d:\ else printf(\请输入a%d%d:\ scanf(\ } k=1; label:
/****************************************************/ ark=max(k);
/*if(fabs(A[1][k])>=fabs(A[2][k])) { ark=fabs(A[1][k]);flag=1;} else { ark=fabs(A[2][k]);flag=2;} if(fabs(A[3][k])>ark)
{ ark=fabs(A[3][k]);flag=3; } */
/**************************************************/ if(ark==0)
{ printf(\不是合法的方程组!\ else if(flag!=k) exchange(flag,k); for(i=k+1;i<=n;i++) for(j=k+1;j<=n+1;j++)
A[i][j]=A[i][j]-A[i][k]*A[k][j]/A[k][k]; if(k { k++; goto label;} else { x[n]=A[n][n+1]/A[n][n]; for( k=n-1;k>=1;k--) {float me=0; for(j=k+1;j<=n;j++) { me=me+A[k][j]*x[j];} x[k]=(A[k][n+1]-me)/A[k][k]; } } for(i=1;i<=n;i++) printf(\ goto me; } /////////////////////////////////////////////////////////////////// 2）//插值源程序 # define f(x) 1/(1+25*x*x) # define f1(x) -50*x/((1+25*x*x)*(1+25*x*x)) # include \float xi[11],M[11]; 8 int N=10; float h=0.2; float L(int i,float x) { float ll=1.0;int j; for( j=0;j<=N;j++) { if(i==j) continue; else ll=ll*(x-xi[j]); } for(j=0;j<=N;j++) { if(i==j) continue; else ll=ll/(xi[i]-xi[j]); } return ll; } float S(int N,float x) { float ss=0; int j; for( j=0;j<=N;++j) ss=ss+f(xi[j])*L(j,x); return ss; } float S1(float x0,float x1,float x) { float ss; ss=(x-x1)/(x0-x1)*f(x0)+(x-x0)/(x1-x0)*f(x1); return ss;} float chashan(float x0,float x1,float x2) { return ((f(x2)-f(x0))/(x2-x0)-(f(x1)-f(x0))/(x1-x0))/(x2-x1); } float S3( float xj,float xjp1 ,float Mj,float Mjp1,float x) { float ss1,ss2,ss3,ss4,ss; ss1=(xjp1-x)*(xjp1-x)*(xjp1-x)/(6*h)*Mj; ss2=(x-xj)*(x-xj)*(x-xj)/(6*h)*Mjp1; ss3=(x-xj)/h*(f(xjp1)-h*h/6*Mjp1); ss4=(xjp1-x)/h*(f(xj)-h*h/6*Mj); ss=ss1+ss2+ss3+ss4; return ss; } main() { int i,j;int control;float help1,help2,helpe,helpb; float d[11]; int cchar=0; float x0,xo1[10],xo2[10]; x0=0; for( i=0;i<=N;i++) xi[i]=-1.0+i*2.0/N; for( j=0;j<=9;++j) 9 { xo1[j]=0.06+0.1*j; xo2[j]=-0.06-0.1*j; } label: printf(\现象的克服\printf(\多项式插值\printf(\分段线性插值\ printf(\三次样条函数插值\printf(\退出\printf(\请选择:\ scanf(\switch(control) {case 1: //下面是多项式插值 for(i=0;i<=9;i++){ printf(\ printf(\ printf(\ printf(\ help1=f(xo1[i]);help2=S(N,xo1[i]); helpe=help1-help2; printf(\helpb=helpe/help1*100.0; printf(\相对误差=%9.6f%% \ } break; case 2: // 下面是分段线性插值 for(j=0;j<=9;j++) for(i=0;i<=N;i++) {if(xo1[j]>xi[i] && xo1[j] help2=S1(xi[i],xi[i+1],xo1[j]); helpe=help1-help2; helpb=helpe/help1*100.0; printf(\ printf(\ printf(\printf(\相对误差=%9.6f%% \ } if(xo2[j]>xi[i] && xo2[j] help2=S1(xi[i],xi[i+1],xo2[j]); printf(\ printf(\ } 10