内容发布更新时间 : 2025/10/31 19:51:59星期一 下面是文章的全部内容请认真阅读。
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