正在为您跳转...
二十四橋明月夜
正文
2018-08-06 17:20:33 | 分类:默认分类 | 标签:结构专业

Doolittle直接分解算法

算法伪码:  A[n][n]—系数矩阵、b[n] —常数项矩阵  L[n][n]、U[n][n]、y[n]、x[n]—L、U矩阵与解向量

Step1 L赋值:对角线元素赋值1,上三角元素赋值0;  Step2 U赋值:下三角元素赋值为0;  Step3 For k=1 To n Do Step4、Step5 /开始分解 /  Step4  For j=k To n Do /求 U中第 k 行元素/  u[k][j] = a[k][j]  For r = 1 To k-1 Do  u[k][j] = u[k][j] - l[k][r]*u[r][j]  EndFor r  EndFor j  Step5  For i=k+1 To n Do /求 L 中第 k 列元素/  l[i][k] = a[i][k];  For r = 1 To k-1 Do  l[i][k] = l[i][k] - l[i][r]*u[r][k];  EndFor r  l[i][k] = l[i][k] / u[k][k]  EndFor i  EndFor k /结束分解 /  Step6  For i = 1 To n Do /由 Ly=b 求出y[i] /  y[i] = b[i]  For j = 1 To i-1 Do  y[i] = y[i] - l[i][j]*y[j];  EndFor j  EndFor i  Step7  For i = n Downto 1 Do /由 Ux=y 求出 x[i] /  x[i] = y[i];  For j = i+1 To n Do  y[i] = y[i] - u[i][j]*x[j]  EndFor j  x[i] = y[i]/u[i][i]  EndFor i


void Doolittle(int n,double *A,double *b)//n为阶数 A为系数矩阵 b为常数矩阵

{

    double *L = new double[n*n];//开辟L矩阵空间

    double *U = new double[n*n];//开辟U矩阵空间

    double *y = new double[n];//开辟y矩阵空间

    for (int i = 0; i < n; i++)

    {

        for (int j = 0; j < n; j++)

        {

            *(U + i*n + j) = 0;//暂时全部赋值为0

            if (i==j)

            {

                *(L + i*n + j) = 1;//对角线赋值为1

            }

            else

            {

                *(L + i*n + j) = 0;//其他暂时赋值为0

            }

        }

    }

    for (int k = 0; k < n; k++)//计算u和l矩阵的值

    {

        for (int j = k; j < n; j++)

        {

            *(U + k*n + j) = *(A + k*n + j);//第一行

            for (int r = 0; r < k; r++)//接下来由L的前一列算u的下一行

            {

                *(U + k*n + j) = *(U + k*n + j) - (*(L + k*n + r)*(*(U + r*n + j)));

            }

        }

        for (int i = k+1; i < n; i++)//计算L的列

        {

            *(L + i*n + k) = *(A + i*n + k);

            for (int r = 0; r < k; r++)

            {

                *(L + i*n + k) = *(L + i*n + k) - (*(L + i*n + r)*(*(U + r*n + k)));

            }

            *(L + i*n + k) = *(L + i*n + k) / (*(U + k*n + k));

        }

    }

    for (int i = 0; i < n; i++)//由Ly=b算y

    {

        *(y + i) = *(b + i);

        for (int j = 0; j < i; j++)

        {

            *(y + i) = *(y + i) - *(L + i*n + j)*(*(y + j));

 

        }

    }

    for (int i = n-1; i >= 0; i--)//由Ux=y算x

    {

        *(x + i) = *(y + i);

        for (int j = i+1; j < n; j++)

        {

            *(y + i) = *(y + i) - *(U + i*n + j)*(*(x + j));

        }

        *(x + i) = *(y + i) / (*(U + i*n + i));

    }

    cout << "解:\n";//得出解

    for (int i = 0; i < n; i++)

    {

        cout <<"x"<<i+1<<":"<< *(x + i) << endl;

    }

    delete[]L;//释放空间

    delete[]U;

    delete[]y;

}



追赶法——三对角线方程组求解

注:三对角线方程组形如 算法伪码:

数据说明:  定义一维数组 a[n],b[n],c[n],l[n],u[n],d[n],x[n],y[n]


Step1  输入数组a、b、c、d各元素  Step2  a[1]=0 ; l[1]=b[1];  u[1]=c[1]/l[1] ; y[1]=d[1]/l[1];  Step3  For i = 1 To n Do /* 追的过程,求 yi */  l[i] = b[i]-a[i]*u[i-1]  u[i] = c[i]/l[i]  y[i] = (d[i]-a[i]*y[i-1])/l[i]  EndFor i  Step4  x[n] = y[n] /* 赶的过程,求 xi */  For i = n-1 Downto 1 Do  x[i] = y[i]-u[i]*x[i+1]  EndFor i


void Catch_Solution(int n,double *a,double *b,double *c,double *d)//参数依次为:阶数 系数a 系数b 系数c 常数列d

{

    double *l = new double[n];

    double *u = new double[n];

    double *x = new double[n];

    double *y = new double[n];//开空间存储数据

    a[0] = 0;

    l[0] = b[0];

    u[0] = c[0] / l[0];

    y[0] = d[0] / l[0];

    for (int i = 1; i < n; i++)

    {

        l[i] = b[i] - a[i] * u[i - 1];

        u[i] = c[i] / l[i];

        y[i] = (d[i] - a[i] * y[i - 1]) / l[i];

    }

    x[n - 1] = y[n - 1];

    for (int i = n-1; i >= 0; i--)//求x

    {

        x[i] = y[i] - u[i] * x[i + 1];

    }

    for (int i = 0; i < n; i++)

    {

        cout << x[i] << endl;

    }

    delete[]x;

    delete[]y;

    delete[]l;

    delete[]u;

}


LDLT分解法

注:要求系数矩阵为对称正定矩阵


算法伪码:

数据说明  a[n][n] ——存放系数矩阵A的系数;  b[n] ——存放方程组右端常数项;  x[n],y[n],z[n] ——即LTx=y,Dy=z,Lz=b;  L[n][n] ——存放L矩阵中的元素;  D[n] ——存放D矩阵中的元素;


Step1  输入系数矩阵、常数项矩阵元素a[i][j],b[i]  Step2  按下面的公式计算(k=1,…,n):  /* L对角线元素赋值为1*/  Step2-1  For i = 1 To n Do  l[i][i]=1  EndFor i  Step2-2  For k = 1 To n Do  d[k] = a[k][k]  For j=1 To k-1 Do  d[k] = d[k]-l[k][j]*l[k][j]*d[j]  EndFor j  For i = k+1 To n Do  l[i][k]=a[i][k]  For j = 1 To k-1 Do  l[i][k] = l[i][k]-l[i][j]*l[k][j]*d[j]  EndFor j  l[i][k] = l[i][k]/d[k]  EndFor i  EndFor k  Step3  For i = 1 To n Do /* 求zi */  z[i] = b[i]  For j = 1 To i-1 Do  z[i] = z[i]-l[i][j]*z[j]  EndFor j  EndFor i  Step4  For i = 1 To n Do /* 求 yi */  y[i]=z[i]/d[i]  EndFor i  Step5  For i = n Downto 1 Do /* 求 xi */  x[i] = y[i];  For j = i+1 To n Do  x[i] = x[i]-l[j][i]*x[j]  EndFor j  EndFor i


void LDL(int n, double *a, double *b)//LDL法,参数依次:阶数 系数矩阵a 常数矩阵b

{

    double *U = new double[n*n];

    double *y = new double[n];

    double *z = new double[n];

    double *L = new double[n*n];

    double *D = new double[n];

for (int i = 0; i < n; i++)//用LU先算出L U

    {

        for (int j = 0; j < n; j++)

        {

            *(U + i*n + j) = 0;//暂时全部赋值为0

            if (i == j)

            {

                *(L + i*n + j) = 1;//对角线赋值为1

            }

            else

            {

                *(L + i*n + j) = 0;//其他暂时赋值为0

            }

        }

    }

 

    for (int k = 0; k < n; k++)//计算u和l矩阵的值

    {

        for (int j = k; j < n; j++)

        {

            *(U + k*n + j) = *(a + k*n + j);//第一行

            for (int r = 0; r < k; r++)//接下来由L的前一列算u的下一行

            {

                *(U + k*n + j) = *(U + k*n + j) - (*(L + k*n + r)*(*(U + r*n + j)));

            }

        }

        for (int i = k + 1; i < n; i++)//计算L的列

        {

            *(L + i*n + k) = *(a + i*n + k);

            for (int r = 0; r < k; r++)

            {

                *(L + i*n + k) = *(L + i*n + k) - (*(L + i*n + r)*(*(U + r*n + k)));

            }

            *(L + i*n + k) = *(L + i*n + k) / (*(U + k*n + k));

        }

    }

 

 

    for (int i = 0; i < n; i++)//把D赋值

    {

        *(D + i) = *(U + i*n + i);

    }

 

 

 

    for (int i = 0; i < n; i++)//由Lz=b算z

    {

        *(z + i) = *(b + i);

        for (int j = 0; j < i; j++)

        {

            *(z + i) = *(z + i) - *(L + i*n + j)*(*(z + j));

 

        }

    }

 

 

 

 

    for (int i = 0; i < n; i++)//算y

    {

        *(y + i) = *(z + i) / (*(D + i));

    }

 

    double *temp = new double[n*n];

    for (int i = 0; i < n; i++)//这里实现对L的转置

    {

        for (int j = 0; j < n; j++)

        {

            *(temp + i*n + j) = *(L + j*n + i);

        }

 

    }

    for (int i = 0; i < n; i++)

    {

        for (int j = 0; j < n; j++)

        {

            *(L + i*n + j) = *(temp + i*n + j);

        }

 

    }

    delete[]temp;//释放

 

    for (int i = n-1; i >= 0; i--)//最后算x

    {

        *(x + i) = *(y + i);

        for (int j = i+1; j <n; j++)

        {

            *(x + i) = *(x + i) - *(L + i*n + j)*(*(x + j));

        }

    }

    for (int i = 0; i < n; i++)

    {

        cout << "解为:\n";

        cout << *(x + i) << endl;

    }

    delete[]U;

    delete[]y;

    delete[]z;

    delete[]L;

    delete[]D;

}



高斯列主元法

算法伪码

数据说明:  N—所能求解方程组的最大阶数;  A[N][N]—表示系数矩阵;b[N]—表示常数项矩阵;  x[N]—表示解向量;n—方程组实际阶数;  Aug[N][N+1]—表示由A和b构成的增广矩阵;  Mtemp[N+1]—用于交换行;r—主元所在行;Pe—主元


Step1  由A[n][n]和b[n]生成增广矩阵Aug[n][n+1]  Step2 For k=1 To n-1 Do  /* 开始消元 */  Step3  r=k ; Pe=fabs(Aug[k][k]) /* 寻找主元 */  Step4  For i=k+1 To n Do  If Pe<  fabs(Aug[i][k]) Then  Pe=fabs(Aug[i][k]) ; r=i  EndIf  EndFor  Step5  If r< >k Then  Swap(r , k) /交换r,k两行 /  EndIf  Step6  If Aug[k][k]=0 Then /主元为0,失败/  Output(“Method Failed!”) ; Stop  EndIf  /* 将第 k 列中从 Aug[k+1][k] 开始的元素消为 0 */  Step7  For i=k+1 To n Do  m=Aug[i][k]/Aug[k][k] /* 取倍数m */  For j=k To n+1 Do  Aug[i][j]=Aug[i][j]-m*Aug[k][j]  EndFor j  EndFor i  EndFor k /* 消元完成*/  Step8  从 x[n] 开始逐步回代,求出解向量 x  x[n]=Aug[n][n+1]/Aug[n][n]  For i=n-1 Downto 1 Do  sum=0.0  For j=i+1 To n Do  sum=sum+Aug[i][j]*x[j]  EndFor j  x[i]=(Aug[i][n]-sum)/Aug[i][i]  EndFor i 


void Swap(int a,int b,double *Aug,int n)//交换矩阵行 参数依次:a,b行(交换行) Aug交换的矩阵 n阶数

{

    double *temp = new double[n + 1];

    for (int i = 0; i <n+1; i++)

    {

        *(temp + i) = *(Aug + a*(n + 1) + i);

    }

    for (int i = 0; i <n + 1; i++)

    {

        *(Aug + a*(n + 1) + i) = *(Aug + b*(n + 1) + i);

    }

    for (int i = 0; i <n + 1; i++)

    {

        *(Aug + b*(n + 1) + i) = *(temp + i);

    }

    delete[]temp;

}

void Guss(int n,double *Aug,double *x)//高斯列主元法,n为阶数,Aug为增广矩阵(由系数矩阵和常数矩阵合并),x为解向量

{

    for (int k = 0; k < n-1; k++)

    {

        double pe;

        int r;

        r = k;

        pe = fabs(*(Aug + k*(n + 1) + k));//假定k行k列为最大主元

        for (int i = k+1; i < n; i++)//寻找最大主元

        {

            if (pe<fabs(*(Aug + i*(n + 1) + k)))

            {

                pe = abs(*(Aug + i*(n + 1) + k)); r = i;

            }

        }

 

        if (r!=k)//后面的行为最大主元时,交换,使主元始终在前

        {

            Swap(r, k,Aug,n);

        }

        if (pe==0)

        {

            cout << "方法失败!"; break;

        }

        for (int i = k+1; i < n; i++)//开始化简矩阵

        {

            double m;

            m = *(Aug + i*(n + 1) + k) / pe;

            for (int j = k; j < n+1; j++)

            {

                *(Aug + i*(n + 1) + j) = *(Aug + i*(n + 1) + j) - m*(*(Aug + k*(n + 1) + j));

            }

        }

    }

 

 

    *(x + n) = *(Aug + (n-1)*(n + 1) + n + 1) / (*(Aug + (n-1)*(n + 1) + n));//计算出最后一个未知数的值,再往回带

    for (int i = n-1; i > 0; i--)//回带公式的实现

    {

        double sum = 0;

        for (int j = i+1; j <= n; j++)

        {

            double m = (*(x + 1));

            sum +=( *(Aug + (i - 1)*(n + 1) + j-1))*(*(x + j - 1));

        }

        *(x + i-1) = (*(Aug + (i-1)*(n + 1) + n) - sum) / (*(Aug + (i-1)*(n + 1) + (i-1)));

    }

    cout << "增广矩阵阵为:\n";

    for (int i = 0; i < n; i++)

    {

        for (int j = 0; j < n+1; j++)

        {

            cout << *(Aug + i*(n + 1) + j) << '\t';

        }

        cout << endl;

    }

    cout << "解向量为:\n";

    for (int i = 0; i < n; i++)

    {

        cout << "x"<<i+1<<":"<<*(x + i) << endl;

    }

}

 


<< 上一篇:《C++隐藏程序主界面》
评论
遵守文明上网宣言,净化网络环境
声明 本平台提供的所有内容仅供学习、分享与交流,平台无法保证内容一定是正确可靠的。通过使用本平台内容而带来的风险与本平台无关,请知悉。