数值计算程序大放送-线性代数方程组

数值计算程序大放送-线性代数方程组

选自<<徐世良数值计算程序集(C)>>

每个程序都加上了适当地注释,陆陆续续干了几个月才整理出来的啊。

今天都给贴出来了

#include "stdlib.h"
#include "math.h"
#include "stdio.h"
// 全选主元高斯消去法
//a-n*n 存放方程组的系数矩阵,返回时将被破坏
//b-常数向量
//x-返回方程组的解向量
//n-存放方程组的阶数
//返回0表示原方程组的系数矩阵奇异
int cagaus(double a[],double b[],int n,double x[])
{
int *js,l,k,i,j,is,p,q;
    double d,t;
    js=malloc(n*sizeof(int));
    l=1;
    for (k=0;k<=n-2;k++)
{
  d=0.0;
        for (i=k;i<=n-1;i++)
  {
   for (j=k;j<=n-1;j++)
            {
    t=fabs(a[i*n+j]);
    if (t>d)
    {
     d=t;
     js[k]=j;
     is=i;
    }
            }
  }
  if (d+1.0==1.0)
  {
   l=0;
  }
  else
  {
   if (js[k]!=k)
   {
    for (i=0;i<=n-1;i++)
    {
     p=i*n+k;
     q=i*n+js[k];
     t=a[p];
     a[p]=a[q];
     a[q]=t;
    }
   }
   if (is!=k)
   {
    for (j=k;j<=n-1;j++)
    {
     p=k*n+j;
     q=is*n+j;
     t=a[p];
     a[p]=a[q];
     a[q]=t;
    }
    t=b[k];
    b[k]=b[is];
    b[is]=t;
   }
  }
  if (l==0)
  {
   free(js);
   printf("fail/n");
   return(0);
  }
  d=a[k*n+k];
  for (j=k+1;j<=n-1;j++)
  {
   p=k*n+j;
   a[p]=a[p]/d;
  }
  b[k]=b[k]/d;
  for (i=k+1;i<=n-1;i++)
  {
   for (j=k+1;j<=n-1;j++)
   {
    p=i*n+j;
    a[p]=a[p]-a[i*n+k]*a[k*n+j];
   }
   b[i]=b[i]-a[i*n+k]*b[k];
  }
}
    d=a[(n-1)*n+n-1];
    if (fabs(d)+1.0==1.0)
{
  free(js);
  printf("fail/n");
        return(0);
}
    x[n-1]=b[n-1]/d;
    for (i=n-2;i>=0;i--)
{
  t=0.0;
        for (j=i+1;j<=n-1;j++)
  {
   t=t+a[i*n+j]*x[j];
  }
        x[i]=b[i]-t;
}
    js[n-1]=n-1;
    for (k=n-1;k>=0;k--)
{
  if (js[k]!=k)
        {
   t=x[k];
   x[k]=x[js[k]];
   x[js[k]]=t;
  }
}
free(js);

 

//复系数方程组的全选主元高斯消去法
//ar-n*n存放复系数矩阵的实部,要被破坏
//ai-n*n存放复系数矩阵的虚部,要被破坏
//n-方程组的阶数
//br-存放方程组右端复常数向量的实部,返回时存放解向量的实部
//bi-存放方程组右端复常数向量的虚部,返回时存放解向量的虚部
//返回值0表示系数矩阵奇异
int cbcgaus(double ar[],double ai[],int n,double br[],double bi[])
{
int *js,l,k,i,j,is,u,v;
    double p,q,s,d;
    js=malloc(n*sizeof(int));
    for (k=0;k<=n-2;k++)
{
  d=0.0;
        for (i=k;i<=n-1;i++)
        {
   for (j=k;j<=n-1;j++)
   {
    u=i*n+j;
    p=ar[u]*ar[u]+ai[u]*ai[u];
    if (p>d)
    {
     d=p;
     js[k]=j;
     is=i;
    }
   }
        }
        if (d+1.0==1.0)
  {
   free(js);
   printf("err**fail/n");
            return(0);
  }
        if (is!=k)
  {
   for (j=k;j<=n-1;j++)
   {
    u=k*n+j;
    v=is*n+j;
                p=ar[u];
                ar[u]=ar[v];
                ar[v]=p;
                p=ai[u];
                ai[u]=ai[v];
                ai[v]=p;
   }
            p=br[k];
            br[k]=br[is];
            br[is]=p;
            p=bi[k];
            bi[k]=bi[is];
            bi[is]=p;
  }
        if (js[k]!=k)
        {
   for (i=0;i<=n-1;i++)
            {
    u=i*n+k;
    v=i*n+js[k];
    p=ar[u];
    ar[u]=ar[v];
    ar[v]=p;
    p=ai[u];
    ai[u]=ai[v];
    ai[v]=p;
            }
        }
        v=k*n+k;
        for (j=k+1;j<=n-1;j++)
  {
   u=k*n+j;
            p=ar[u]*ar[v];
   q=-ai[u]*ai[v];
            s=(ar[v]-ai[v])*(ar[u]+ai[u]);
            ar[u]=(p-q)/d;
   ai[u]=(s-p-q)/d;
  }
        p=br[k]*ar[v];
  q=-bi[k]*ai[v];
        s=(ar[v]-ai[v])*(br[k]+bi[k]);
        br[k]=(p-q)/d;
  bi[k]=(s-p-q)/d;
        for (i=k+1;i<=n-1;i++)
  {
   u=i*n+k;
            for (j=k+1;j<=n-1;j++)
   {
    v=k*n+j;
    l=i*n+j;
                p=ar[u]*ar[v];
    q=ai[u]*ai[v];
                s=(ar[u]+ai[u])*(ar[v]+ai[v]);
                ar[l]=ar[l]-p+q;
                ai[l]=ai[l]-s+p+q;
   }
            p=ar[u]*br[k];
   q=ai[u]*bi[k];
            s=(ar[u]+ai[u])*(br[k]+bi[k]);
            br[i]=br[i]-p+q;
   bi[i]=bi[i]-s+p+q;
  }
}
    u=(n-1)*n+n-1;
    d=ar[u]*ar[u]+ai[u]*ai[u];
    if (d+1.0==1.0)
{
  free(js);
  printf("err**fail/n");
        return(0);
}
    p=ar[u]*br[n-1];
q=-ai[u]*bi[n-1];
    s=(ar[u]-ai[u])*(br[n-1]+bi[n-1]);
    br[n-1]=(p-q)/d;
bi[n-1]=(s-p-q)/d;
    for (i=n-2;i>=0;i--)
    {
  for (j=i+1;j<=n-1;j++)
  {
   u=i*n+j;
   p=ar[u]*br[j];
   q=ai[u]*bi[j];
   s=(ar[u]+ai[u])*(br[j]+bi[j]);
   br[i]=br[i]-p+q;
   bi[i]=bi[i]-s+p+q;
  }
}
    js[n-1]=n-1;
    for (k=n-1;k>=0;k--)
{
  if (js[k]!=k)
        {
   p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
   p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
        }
}
free(js);
return(1);
}

//全选主元高斯-约当消去法
//a-n*n 存放方程组的系数矩阵,返回时将被破坏
//b-n*m常数向量,每列为一组,返回时存放m组解向量
//n-方程组的阶数
//m-方程组右端常数向量的个数
//返回0表示原方程组的系数矩阵奇异
int ccgj(double a[],double b[],int n,int m)
{
int *js,l,k,i,j,is,p,q;
    double d,t;
    js=malloc(n*sizeof(int));
    l=1;
    for (k=0;k<=n-1;k++)
{
  d=0.0;
        for (i=k;i<=n-1;i++)
  {
   for (j=k;j<=n-1;j++)
            {
    t=fabs(a[i*n+j]);
    if (t>d)
    {
     d=t;js[k]=j;is=i;
    }
            }
  }
        if (d+1.0==1.0)
  {
   l=0;
  }
        else
  {
   if (js[k]!=k)
   {
    for (i=0;i<=n-1;i++)
                {
     p=i*n+k; q=i*n+js[k]; t=a[p]; a[p]=a[q]; a[q]=t;
                }
   }
   if (is!=k)
   {
    for (j=k;j<=n-1;j++)
    {
     p=k*n+j; q=is*n+j; t=a[p]; a[p]=a[q]; a[q]=t;
    }
    for (j=0;j<=m-1;j++)
    {
     p=k*m+j; q=is*m+j; t=b[p]; b[p]=b[q]; b[q]=t;
    }
   }
  }
        if (l==0)
  {
   free(js); printf("fail/n");
            return(0);
  }
        d=a[k*n+k];
        for (j=k+1;j<=n-1;j++)
  {
   p=k*n+j; a[p]=a[p]/d;
  }
        for (j=0;j<=m-1;j++)
  {
   p=k*m+j; b[p]=b[p]/d;
  }
        for (j=k+1;j<=n-1;j++)
  {
   for (i=0;i<=n-1;i++)
            {
    p=i*n+j;
    if (i!=k)
    {
     a[p]=a[p]-a[i*n+k]*a[k*n+j];
    }
            }
  }
        for (j=0;j<=m-1;j++)
  {
   for (i=0;i<=n-1;i++)
   {
    p=i*m+j;
    if (i!=k)
    {
     b[p]=b[p]-a[i*n+k]*b[k*m+j];
    }
   }
  }
}
    for (k=n-1;k>=0;k--)
{
  if (js[k]!=k)
  {
   for (j=0;j<=m-1;j++)
   {
    p=k*m+j; q=js[k]*m+j; t=b[p]; b[p]=b[q]; b[q]=t;
   }
  }
}
    free(js);
    return(1);
}

//复系数全选主元高斯-约当消去法
//ar-n*n 存放方程组的系数矩阵的实部,返回时将被破坏
//ai-n*n 存放方程组的系数矩阵的虚部,返回时将被破坏
//br-n*m常数向量的实部,每列为一组,返回时存放m组解向量的实部
//bi-n*m常数向量的虚部,每列为一组,返回时存放m组解向量的虚部
//n-方程组的阶数
//m-方程组右端常数向量的个数
//返回0表示原方程组的系数矩阵奇异
int cdcgj(double ar[],double ai[],double br[],double bi[],int n,int m)
{
int *js,l,k,i,j,is,u,v;
    double p,q,s,d;
    js=malloc(n*sizeof(int));
    for (k=0;k<=n-1;k++)
{
  d=0.0;
        for (i=k;i<=n-1;i++)
  {
   for (j=k;j<=n-1;j++)
   {
    u=i*n+j;
    p=ar[u]*ar[u]+ai[u]*ai[u];
    if (p>d)
    {
     d=p;js[k]=j;is=i;
    }
   }
  }
        if (d+1.0==1.0)
  {
   free(js);
   printf("err**fail/n");
            return(0);
  }
        if (is!=k)
  {
   for (j=k;j<=n-1;j++)
   {
    u=k*n+j; v=is*n+j;
                p=ar[u]; ar[u]=ar[v]; ar[v]=p;
                p=ai[u]; ai[u]=ai[v]; ai[v]=p;
   }
            for (j=0;j<=m-1;j++)
   {
    u=k*m+j; v=is*m+j;
                p=br[u]; br[u]=br[v]; br[v]=p;
                p=bi[u]; bi[u]=bi[v]; bi[v]=p;
   }
  }
        if (js[k]!=k)
  {
   for (i=0;i<=n-1;i++)
            {
    u=i*n+k; v=i*n+js[k];
    p=ar[u]; ar[u]=ar[v]; ar[v]=p;
    p=ai[u]; ai[u]=ai[v]; ai[v]=p;
            }
  }
        v=k*n+k;
        for (j=k+1;j<=n-1;j++)
  {
   u=k*n+j;
            p=ar[u]*ar[v]; q=-ai[u]*ai[v];
            s=(ar[v]-ai[v])*(ar[u]+ai[u]);
            ar[u]=(p-q)/d; ai[u]=(s-p-q)/d;
  }
        for (j=0;j<=m-1;j++)
  {
   u=k*m+j;
            p=br[u]*ar[v]; q=-bi[u]*ai[v];
            s=(ar[v]-ai[v])*(br[u]+bi[u]);
            br[u]=(p-q)/d; bi[u]=(s-p-q)/d;
  }
        for (i=0;i<=n-1;i++)
  {
   if (i!=k)
            {
    u=i*n+k;
    for (j=k+1;j<=n-1;j++)
                {
     v=k*n+j; l=i*n+j;
     p=ar[u]*ar[v]; q=ai[u]*ai[v];
     s=(ar[u]+ai[u])*(ar[v]+ai[v]);
     ar[l]=ar[l]-p+q;
     ai[l]=ai[l]-s+p+q;
                }
    for (j=0;j<=m-1;j++)
                {
     l=i*m+j; v=k*m+j;
     p=ar[u]*br[v]; q=ai[u]*bi[v];
     s=(ar[u]+ai[u])*(br[v]+bi[v]);
     br[l]=br[l]-p+q; bi[l]=bi[l]-s+p+q;
                }
            }
  }
}
    for (k=n-1;k>=0;k--)
{
  if (js[k]!=k)
  {
   for (j=0;j<=m-1;j++)
   {
    u=k*m+j;v=js[k]*m+j;
    p=br[u]; br[u]=br[v]; br[v]=p;
    p=bi[u]; bi[u]=bi[v]; bi[v]=p;
   }
  }
}
    free(js);
    return(1);
}

//求解三对角线方程组的追赶法
//存放三对角线矩阵中三对角线上的元素
//存放顺序为a00,a01,a10,a11,a12,a21,a22,a23,...,an-1,n-2,an-1,n-1
//就是按行依次存储
//n-方程组的阶数
//m-b的长度,应为m=3n-2,本函数中用来做检验
//d-存放方程组右端常数向量,返回时存放解向量
int cetrd(double b[],int n,int m,double d[])
{
int k,j;
    double s;
    if (m!=(3*n-2))
{
  printf("err/n"); return(-2);
}
    for (k=0;k<=n-2;k++)
{
  j=3*k; s=b[j];
        if (fabs(s)+1.0==1.0)
  {
   printf("fail/n"); return(0);
  }
        b[j+1]=b[j+1]/s;
        d[k]=d[k]/s;
        b[j+3]=b[j+3]-b[j+2]*b[j+1];
        d[k+1]=d[k+1]-b[j+2]*d[k];
}
    s=b[3*n-3];
    if (fabs(s)+1.0==1.0)
{
  printf("fail/n"); return(0);
}
    d[n-1]=d[n-1]/s;
    for (k=n-2;k>=0;k--)
{
  d[k]=d[k]-b[3*k+1]*d[k+1];
}
    return(2);
}

//一般带型方程组的求解
//n-方程组的阶数
//l-系数矩阵的半带宽
//il-系数矩阵的带宽,应为il=2l+1
//m-方程组右端的常数
int cfband(double b[],double d[],int n,int l,int il,int m)
{
int ls,k,i,j,is,u,v;
    double p,t;
    if (il!=(2*l+1))
{
  printf("fail/n"); return(-2);
}
    ls=l;
    for (k=0;k<=n-2;k++)
{
  p=0.0;
        for (i=k;i<=ls;i++)
  {
   t=fabs(b[i*il]);
            if (t>p)
   {
    p=t; is=i;
   }
  }
        if (p+1.0==1.0)
  {
   printf("fail/n"); return(0);
  }
        for (j=0;j<=m-1;j++)
  {
   u=k*m+j; v=is*m+j;
            t=d[u]; d[u]=d[v]; d[v]=t;
  }
        for (j=0;j<=il-1;j++)
  {
   u=k*il+j; v=is*il+j;
            t=b[u]; b[u]=b[v]; b[v]=t;
  }
        for (j=0;j<=m-1;j++)
  {
   u=k*m+j; d[u]=d[u]/b[k*il];
  }
        for (j=1;j<=il-1;j++)
  {
   u=k*il+j; b[u]=b[u]/b[k*il];
  }
        for (i=k+1;i<=ls;i++)
  {
   t=b[i*il];
            for (j=0;j<=m-1;j++)
   {
    u=i*m+j; v=k*m+j;
                d[u]=d[u]-t*d[v];
   }
            for (j=1;j<=il-1;j++)
   {
    u=i*il+j; v=k*il+j;
                b[u-1]=b[u]-t*b[v];
   }
            u=i*il+il-1; b[u]=0.0;
  }
        if (ls!=(n-1))
  {
   ls=ls+1;
  }
}
    p=b[(n-1)*il];
    if (fabs(p)+1.0==1.0)
{
  printf("fail/n"); return(0);
}
    for (j=0;j<=m-1;j++)
{
  u=(n-1)*m+j; d[u]=d[u]/p;
}
    ls=1;
    for (i=n-2;i>=0;i--)
{
  for (k=0;k<=m-1;k++)
  {
   u=i*m+k;
            for (j=1;j<=ls;j++)
   {
    v=i*il+j; is=(i+j)*m+k;
                d[u]=d[u]-b[v]*d[is];
   }
  }
        if (ls!=(il-1)) ls=ls+1;
}
    return(2);
}

//求解对称方程组的分解法
//a n*n数组,存放方程组的系数矩阵(应为对称矩阵)
//n 方程组的阶数
//m 方程组右端常数向量的个数
//c n*m 存放方程组右端的m组常数向量;返回时存放方程组的m组解向量
//返回值小于0表示程序工作失败,大于0表示正常返回
int cgldl(double a[],int n,int m,double c[])
{
int i,j,l,k,u,v,w,k1,k2,k3;
double p;
if (fabs(a[0])+1.0==1.0)
{
  printf("fail/n");
  return(-2);
}
for (i=1; i<=n-1; i++)
{
  u=i*n;
  a[u]=a[u]/a[0];
}
for (i=1; i<=n-2; i++)
{
  u=i*n+i;
  for (j=1; j<=i; j++)
  {
   v=i*n+j-1;
   l=(j-1)*n+j-1;
   a[u]=a[u]-a[v]*a[v]*a[l];
  }
  p=a[u];
  if (fabs(p)+1.0==1.0)
  {
   printf("fail/n");
   return(-2);
  }
  for (k=i+1; k<=n-1; k++)
  {
   u=k*n+i;
   for (j=1; j<=i; j++)
   {
    v=k*n+j-1;
    l=i*n+j-1;
    w=(j-1)*n+j-1;
    a[u]=a[u]-a[v]*a[l]*a[w];
   }
   a[u]=a[u]/p;
  }
}
u=n*n-1;
for (j=1; j<=n-1; j++)
{
  v=(n-1)*n+j-1;
  w=(j-1)*n+j-1;
  a[u]=a[u]-a[v]*a[v]*a[w];
}
p=a[u];
if (fabs(p)+1.0==1.0)
{
  printf("fail/n");
  return(-2);
}
for (j=0; j<=m-1; j++)
{
  for (i=1; i<=n-1; i++)
  {
   u=i*m+j;
   for (k=1; k<=i; k++)
   {
    v=i*n+k-1;
    w=(k-1)*m+j;
    c[u]=c[u]-a[v]*c[w];
   }
  }
}
for (i=1; i<=n-1; i++)
{
  u=(i-1)*n+i-1;
  for (j=i; j<=n-1; j++)
  {
   v=(i-1)*n+j;
   w=j*n+i-1;
   a[v]=a[u]*a[w];
  }
}
for (j=0; j<=m-1; j++)
{
  u=(n-1)*m+j;
  c[u]=c[u]/p;
  for (k=1; k<=n-1; k++)
  {
   k1=n-k;
   k3=k1-1;
   u=k3*m+j;
   for (k2=k1; k2<=n-1; k2++)
   {
    v=k3*n+k2;
    w=k2*m+j;
    c[u]=c[u]-a[v]*c[w];
   }
   c[u]=c[u]/a[k3*n+k3];
  }
}
return(2);
}

//求解对称正定方程组的平方根法
//用Cholesky分解法(即平方根法)
//a n*n 存放系数矩阵(应为对称正定矩阵),返回时其上三角部分存放分解后的U矩阵
//n 方程的阶数
//m 方程组有短的常数向量的个数
//d n*m 存放方程组右端m组常数向量;返回时,存放方程组的m组解向量
//返回值小于0,表示程序工作失败
int chchol(double a[],int n,int m,double d[])
{
int i,j,k,u,v;
    if ((a[0]+1.0==1.0)||(a[0]<0.0))
{
  printf("fail/n"); return(-2);
}
    a[0]=sqrt(a[0]);
    for (j=1; j<=n-1; j++)
{
  a[j]=a[j]/a[0];
}
    for (i=1; i<=n-1; i++)
{
  u=i*n+i;
        for (j=1; j<=i; j++)
  {
   v=(j-1)*n+i;
            a[u]=a[u]-a[v]*a[v];
  }
        if ((a[u]+1.0==1.0)||(a[u]<0.0))
  {
   printf("fail/n"); return(-2);
  }
        a[u]=sqrt(a[u]);
        if (i!=(n-1))
  {
   for (j=i+1; j<=n-1; j++)
   {
    v=i*n+j;
                for (k=1; k<=i; k++)
    {
     a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
    }
                a[v]=a[v]/a[u];
   }
  }
}
    for (j=0; j<=m-1; j++)
{
  d[j]=d[j]/a[0];
        for (i=1; i<=n-1; i++)
  {
   u=i*n+i; v=i*m+j;
            for (k=1; k<=i; k++)
   {
    d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
   }
            d[v]=d[v]/a[u];
  }
}
    for (j=0; j<=m-1; j++)
{
  u=(n-1)*m+j;
        d[u]=d[u]/a[n*n-1];
        for (k=n-1; k>=1; k--)
  {
   u=(k-1)*m+j;
            for (i=k; i<=n-1; i++)
   {
    v=(k-1)*n+i;
                d[u]=d[u]-a[v]*d[i*m+j];
   }
            v=(k-1)*n+k-1;
            d[u]=d[u]/a[v];
  }
}
    return(2);
}

//求解大型稀疏方程组的全选主元高斯-约当消去法
//只是在消去过程中避免了对零元素的运算,从而可以大大减少运算次数
//本函数不考虑系数矩阵的压缩存储问题
//a n*n数组,存放方程组的系数矩阵,返回时要被破坏
//n 方程组的阶数
//d-存放方程组右端常数向量,返回时存放解向量
//返回值为0表示系数矩阵奇异
int ciggj(double a[],int n,double b[])
{
int *js,i,j,k,is,u,v;
    double d,t;
    js=malloc(n*sizeof(int));
    for (k=0; k<=n-1; k++)
{
  d=0.0;
        for (i=k; i<=n-1; i++)
        {
   for (j=k; j<=n-1; j++)
   {
    t=fabs(a[i*n+j]);
    if (t>d)
    {
     d=t; js[k]=j; is=i;
    }
   }
        }
        if (d+1.0==1.0)
  {
   free(js); printf("fail/n"); return(0);
  }
        if (is!=k)
  {
   for (j=k; j<=n-1; j++)
   {
    u=k*n+j; v=is*n+j;
                t=a[u]; a[u]=a[v]; a[v]=t;
   }
            t=b[k]; b[k]=b[is]; b[is]=t;
  }
        if (js[k]!=k)
        {
   for (i=0; i<=n-1; i++)
            {
    u=i*n+k; v=i*n+js[k];
    t=a[u]; a[u]=a[v]; a[v]=t;
            }
        }
        t=a[k*n+k];
        for (j=k+1; j<=n-1; j++)
  {
   u=k*n+j;
            if (a[u]!=0.0)
            {
    a[u]=a[u]/t;
   }
  }
        b[k]=b[k]/t;
        for (j=k+1; j<=n-1; j++)
  {
   u=k*n+j;
            if (a[u]!=0.0)
   {
    for (i=0; i<=n-1; i++)
    {
     v=i*n+k;
                    if ((i!=k)&&(a[v]!=0.0))
     {
      is=i*n+j;
                        a[is]=a[is]-a[v]*a[u];
     }
    }
   }
  }
        for (i=0; i<=n-1; i++)
  {
   u=i*n+k;
            if ((i!=k)&&(a[u]!=0.0))
            {
    b[i]=b[i]-a[u]*b[k];
   }
  }
}
    for (k=n-1; k>=0; k--)
    {
  if (k!=js[k])
        {
   t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;
  }
}
    free(js);
    return(1);
}

//求解Toeplitz型方程组的Levinson递推算法
//t 长度为n 的一维数组,存放对称T型矩阵中的元素t0,t1,...tn-1
//n 方程组的阶数
//b 长度为n 的一维数组,存放方程组右端的常数向量
//x 长度为n 的一维数组,返回方程组的解
//返回值小于0表示程序工作失败
int cjtlvs(double t[],int n,double b[],double x[])
{
int i,j,k;
    double a,beta,q,c,h,*y,*s;
    s=malloc(n*sizeof(double));
    y=malloc(n*sizeof(double));
    a=t[0];
    if (fabs(a)+1.0==1.0)
{
  free(s); printf("fail/n"); return(-1);
}
    y[0]=1.0; x[0]=b[0]/a;
    for (k=1; k<=n-1; k++)
{
  beta=0.0; q=0.0;
        for (j=0; j<=k-1; j++)
  {
   beta=beta+y[j]*t[j+1];
            q=q+x[j]*t[k-j];
  }
        if (fabs(a)+1.0==1.0)
  {
   free(s); printf("fail/n"); return(-1);
  }
        c=-beta/a; s[0]=c*y[k-1]; y[k]=y[k-1];
        if (k!=1)
  {
   for (i=1; i<=k-1; i++)
   {
    s[i]=y[i-1]+c*y[k-i-1];
   }
  }
        a=a+c*beta;
        if (fabs(a)+1.0==1.0)
  {
   free(s); printf("fail/n"); return(-1);
  }
        h=(b[k]-q)/a;
        for (i=0; i<=k-1; i++)
  {
   x[i]=x[i]+h*s[i]; y[i]=s[i];
  }
        x[k]=h*y[k];
}
    free(s); free(y);
    return(1);
}

 

展开阅读全文

没有更多推荐了,返回首页

©️2019 CSDN 皮肤主题: 技术工厂 设计师: CSDN官方博客
应支付0元
点击重新获取
扫码支付

支付成功即可阅读