exercise1.cpp

数值分析大作业第1题

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
int i,j,k,n;

main ()
{
double a[501],b,c,m,t,n,det,g,cond,l[2],f[2],lambda,uu[501],lb[502],lc[503],ub[501],l1,l501;
for (k=1;k<=501;k++)
  a[k-1]=(1.64-0.024*k)*sin(0.2*k)-0.64*exp(0.1/k);
b=0.16;
c=-0.064;
double u[501],y[501], v[501],z[501],w[501];
u[0]=2;n=0;l[1]=0; i=1;v[0]=2;y[0]=0; z[0]=0;w[0]=0;
for (k=1;k<=500;k++)
{u[k]=1;v[k]=1;y[k]=0;z[k]=0;w[k]=0; }
while (n==0)
  {
   l[0]=l[1];
    l[1]=0;
   for (m=0,k=1;k<=501;k++)
     m=m+(u[k-1])*(u[k-1]);
   t=sqrt(m);

   for (k=1;k<=501;k++)
      y[k-1]=u[k-1]/t;

    u[0]=a[0]*y[0]+b*y[1]+c*y[2];
    u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3];

   for (j=3;j<=499;j++)
       u[j-1]=c*y[j-3]+b*y[j-2]+a[j-1]*y[j-1]+b*y[j]+c*y[j+1];
    u[499]=c*y[497]+b*y[498]+a[499]*y[499]+b*y[500];
    u[500]=c*y[498]+b*y[499]+a[500]*y[500];
    for (l[1]=0,k=1;k<=501;k++)
      l[1]=l[1]+y[k-1]*u[k-1];

    if (((l[1]-l[0])/l[1])<=(1e-12)&&((l[1]-l[0])/l[1])>=(-1e-12))
      n=1;
    i++;
  }
   if (l[1]>0)
   printf ("lambda(501)：%12.11E\n",l[1]);
   if (l[1]<0)
   printf ("lambda(1)：%12.11E\n",l[1]);
   l1=l[1];

    uu[0]=a[0]; uu[1]=a[1]-b*b;
    ub[0]=b;
    lb[1]=b/uu[0]; lc[2]=c/uu[0];
    for (k=2;k<=501;k++)
      {
            ub[k-1]=b-lb[k-1]*c;
            uu[k-1]=a[k-1]-lc[k-1]*c-lb[k-1]*ub[k-2];
            lb[k]=(b-lc[k]*ub[k-2])/uu[k-1];
            lc[k+1]=c/uu[k-1];
      }
     for (k=0,det=1;k<=500;k++)
      det=det*uu[k];
      printf("det(A)：%12.11E\n",det);

v[0]=1;n=0;f[1]=a[0]; i=0;
while (n==0)
  {
   f[0]=f[1];
   for (m=0,k=1;k<=501;k++)
     m=m+(v[k-1])*(v[k-1]);
   t=sqrt(m);
 for (k=0;k<=500;k++)
      z[k]=v[k]/t;

        w[0]=z[0];
      w[1]=z[1]-w[0]*lb[1];

      for (j=2;j<=i,j<=500;j++)
        w[j]=z[j]-lc[j]*w[j-2]-lb[j]*w[j-1];



   if (i<=500)
      { v[i]=w[i]/uu[i];
       v[i-1]=(w[i-1]-ub[i-1]*v[i])/uu[i-1];
      for (j=i;j>=2;j--)
       v[j-2]=(w[j-2]-ub[j-2]*v[j-1]-c*v[j])/uu[j-2];
      }
 if (i>500)
   { v[500]=w[500]/uu[500];
       v[499]=(w[499]-ub[499]*v[500])/uu[499];
      for (j=500;j>=2;j--)
       v[j-2]=(w[j-2]-ub[j-2]*v[j-1]-c*v[j])/uu[j-2];
      }

    for (f[1]=0,k=0;k<=500;k++)
      f[1]=f[1]+z[k]*v[k];

    if (fabs(1-f[1]/f[0])<=(1e-12))
      {n=1;lambda=1/f[1];}
    i++;
   }
printf ("lambda(s)：%12.11E\n",lambda);
cond=fabs(l[1]/lambda);
printf("cond（A）:%12.11E\n",cond);


n=0;v[0]=2;y[0]=0; z[0]=0;w[0]=0;
for (k=1;k<=500;k++)
{v[k]=1;y[k]=0;z[k]=0;w[k]=0; }

 for (k=1;k<=501;k++)
  a[k-1]=(1.64-0.024*k)*sin(0.2*k)-0.64*exp(0.1/k)+l[1];
  uu[0]=a[0]; uu[1]=a[1]-b*b;
    ub[0]=b;
    lb[1]=b/uu[0]; lc[2]=c/uu[0];
    for (k=2;k<=501;k++)
      {
            ub[k-1]=b-lb[k-1]*c;
            uu[k-1]=a[k-1]-lc[k-1]*c-lb[k-1]*ub[k-2];
            lb[k]=(b-lc[k]*ub[k-2])/uu[k-1];
            lc[k+1]=c/uu[k-1];
      }


  v[0]=1;n=0;f[1]=a[0]; i=0;
while (n==0)
  {
   f[0]=f[1];
   for (m=0,k=1;k<=501;k++)
     m=m+(v[k-1])*(v[k-1]);
   t=sqrt(m);
 for (k=0;k<=500;k++)
      z[k]=v[k]/t;

        w[0]=z[0];
      w[1]=z[1]-w[0]*lb[1];

      for (j=2;j<=i,j<=500;j++)
        w[j]=z[j]-lc[j]*w[j-2]-lb[j]*w[j-1];



   if (i<=500)
      { v[i]=w[i]/uu[i];
       v[i-1]=(w[i-1]-ub[i-1]*v[i])/uu[i-1];
      for (j=i;j>=2;j--)
       v[j-2]=(w[j-2]-ub[j-2]*v[j-1]-c*v[j])/uu[j-2];
      }
 if (i>500)
   { v[500]=w[500]/uu[500];
       v[499]=(w[499]-ub[499]*v[500])/uu[499];
      for (j=500;j>=2;j--)
       v[j-2]=(w[j-2]-ub[j-2]*v[j-1]-c*v[j])/uu[j-2];
      }

    for (f[1]=0,k=0;k<=500;k++)
      f[1]=f[1]+z[k]*v[k];
     // printf ("%12.11E\n",((1/f[1])-l[1]));
    if (fabs((1/f[1]-1/f[0])/(1/f[1]-l[1]))<=(1e-12))
      {n=1; lambda=(1/f[1])-l[1]; }

    i++;
   }
    lambda=(1/f[1])-l[1];
    l501=lambda;
  if (n>0)
   printf ("lambda(501)：%12.11E\n",lambda);
   if (n<0)
   printf ("lambda(1)：%12.11E\n",lambda);


 double miu[39];int p;
 for  (p=0;p<=38;p++)
 { miu[p]=l1+(p+1)*(l501-l1)/40;
 n=0;v[0]=2;y[0]=0; z[0]=0;w[0]=0;
for (k=1;k<=500;k++)
{v[k]=1;y[k]=0;z[k]=0;w[k]=0; }

 for (k=1;k<=501;k++)
  a[k-1]=(1.64-0.024*k)*sin(0.2*k)-0.64*exp(0.1/k)-miu[p];
  uu[0]=a[0]; uu[1]=a[1]-b*b;
    ub[0]=b;
    lb[1]=b/uu[0]; lc[2]=c/uu[0];
    for (k=2;k<=501;k++)
      {
            ub[k-1]=b-lb[k-1]*c;
            uu[k-1]=a[k-1]-lc[k-1]*c-lb[k-1]*ub[k-2];
            lb[k]=(b-lc[k]*ub[k-2])/uu[k-1];
            lc[k+1]=c/uu[k-1];
      }


  v[0]=1;n=0;f[1]=0; i=0;
while (n==0)
  {
   f[0]=f[1];
    f[1]=0;
   for (m=0,k=1;k<=501;k++)
     m=m+(v[k-1])*(v[k-1]);
   t=sqrt(m);
 for (k=0;k<=500;k++)
      z[k]=v[k]/t;

        w[0]=z[0];
      w[1]=z[1]-w[0]*lb[1];

      for (j=2;j<=i,j<=500;j++)
        w[j]=z[j]-lc[j]*w[j-2]-lb[j]*w[j-1];



   if (i<=500)
      { v[i]=w[i]/uu[i];
       v[i-1]=(w[i-1]-ub[i-1]*v[i])/uu[i-1];
      for (j=i;j>=2;j--)
       v[j-2]=(w[j-2]-ub[j-2]*v[j-1]-c*v[j])/uu[j-2];
      }
 if (i>500)
   { v[500]=w[500]/uu[500];
       v[499]=(w[499]-ub[499]*v[500])/uu[499];
      for (j=500;j>=2;j--)
       v[j-2]=(w[j-2]-ub[j-2]*v[j-1]-c*v[j])/uu[j-2];
      }

    for (f[1]=0,k=0;k<=500;k++)
      f[1]=f[1]+z[k]*v[k];

    if (fabs((1/f[1]-1/f[0])/(1/f[1]+miu[p]))<=(1e-12))
    n=1;
    i++;
  }
    lambda=(1/f[1])+miu[p];
   printf ("lambda(i%d)：%12.11E\n",p+1,lambda);
 }
}