當前位置:首頁 » 編程語言 » lu分解法c語言

lu分解法c語言

發布時間: 2022-09-18 23:02:12

① 對矩陣x進行QR分解和LU分解,QR分解和LU分解是什麼意思呢

你好!
LU分解是矩陣的三角分解,產生一個上三角矩陣和一個下三角矩陣。
QR分解是矩陣的正交分解。
如果對你有幫助,望採納。

② 矩陣的LU分解該怎麼具體做。親,。舉個例子吧

LU分解在本質上是高斯消元法的一種表達形式。實質上是將A通過初等行變換變成一個上三角矩陣,其變換矩陣就是一個單位下三角矩陣。
(A,E)~r2-2r1
1 -2 3 1 0 0
0 8 -4 -2 1 0
0 1 1 0 0 1
r3-1/8r2
1 -2 3 1 0 0
0 8 -4 -2 1 0
0 0 3/2 1/4 -1/8 1
該矩陣記做(U,P)
求出矩陣P=
1 0 0
-2 1 0
1/4 -1/8 1
的逆P-1=
1 0 0
2 1 0
0 1/8 1
因為PA=U,所以A=P-1U=LU,L=P-1

③ 怎樣用LU分解法解線性方程組

Ax=B,改寫成Ly=B,Ux=y的方程組。就相當於將A=LU分解成了兩個矩陣。稱為矩陣A的三角分解,或LU分解。如果L為單位下三角陣,則叫Doolittle分解,若U為單位上三角陣,則叫Crout分解。只要A的各順序主子式不為零,則A可唯一分解成一個單位下三角陣L與一個上三角陣U的乘積。

•設Ax=b,A=LU,則Ax=LUx=b

於是令Ux=y,則Ly=b

這樣原來方程能化為兩個簡單方程組

在線性代數中, LU分解(LU Decomposition)是矩陣分解的一種,可以將一個矩陣分解為一個單位下三角矩陣和一個上三角矩陣的乘積(有時是它們和一個置換矩陣的乘積)。LU分解主要應用在數值分析中,用來解線性方程、求反矩陣或計算行列式。

(3)lu分解法c語言擴展閱讀:

相關演算法

LU分解在本質上是高斯消元法的一種表達形式。實質上是將A通過初等行變換變成一個上三角矩陣,其變換矩陣就是一個單位下三角矩陣。

這正是所謂的杜爾里特演算法:從下至上地對矩陣A做初等行變換,將對角線左下方的元素變成零,然後再證明這些行變換的效果等同於左乘一系列單位下三角矩陣,這一系列單位下三角矩陣的乘積的逆就是L矩陣,它也是一個單位下三角矩陣。這類演算法的復雜度一般在(三分之二的n三次方) 左右。

④ c語言能解決方程嗎

我學過《數值分析》,李慶揚老師主編,施普辛格出版社。老師上課時雖然沒有講用C語言實現,但是給了流程圖。我學過的用C語言解方程主要是牛頓迭代法,主要解非線性方程;解線性方程組的主要是LU分解法和列主元高斯消去法,這些東西都是我們的前輩早就研究出來的,但是現在可以用計算機方便地實現。

⑤ 計算方法LU分解法求解線性方程組 在線等

其中L為下三角矩陣,U為上三角矩陣。
例如,4×4矩陣A的情況,(1)式如下:
 (2)
可以用如(1)式分解來求解線性方程組
                                  (3)
首先求解向量y使得
                                                                               (4)
然後再來求解
                                                      (5)
此拆分方法的優點在於求解一個三角形方程組相當容易,這樣,(4)式可用向前替代過程求解,如下:
                                                (6)
(5)式可用回代過程求解,這與(2)式~(3)式一樣,
                        (7)
(6)式和(7)式共需執行N2次內層循環(對每個右端項b),每個內層循環包括一次乘法和一次加法。如果有N個右端項,它們是單位列向量(在求矩陣逆時就是這種情況),考慮這些零元素可把(6)式的總執行次數從N3/2減少到N3/6,而(7)式的執行次數不變,仍為N3/2。
       注意:一點對A進行了LU分解,就可以一次求解所有要解的右端項。
 
演算法實現:
首先,寫出(1)式或(2)式的第i,j分量。它總是一個和式,開始部分形式如下:
和式中的項數依賴於i和j中較小的數。事實上有三種形式:
                         (8,9,10)
顯然,(8)~(10)式共有N2個方程,而要求N2+N個未知的α和β(因對角線的未知元素有兩套),既然未知數的個數比方程個數多,就人為指定N各位指數,然後再來求解其他的未知數。事實上,總是令
                                                                          (11)
有一個演算法稱為Crout演算法,它僅按某種次序排列方程,就能容易的求出(8)式~(11)式的N2+N各方程中的所有α和β。步驟如下:

    ,即(11)式
對每個j=0,1,2,...,N-1進行以下兩步:
第一步,對每個i=0,1,...,j用(8)式、(9)式和(11)式來解βij,即
                                                                          (12)
第二步,對每個i=j+1,j+2,...,N-1用(10)式來求解αij,即
                              (13)
在求解下一個j之前要保證進行了以上兩步。
 
如果按上述過程進行幾次迭代後,就會發現(12)式和(13)式右端的α和β在需要時已經得到,還會發現,每一個aij僅被使用一次就不再使用了。這意味著分解是「同址」進行的。簡言之Crout演算法得到的矩陣是混合矩陣,對本例排列如下:

註:不是把矩陣A分解成LU形式,而是將其按行置換的方式分解。
 
Crout演算法的精妙之處:
l         (12)式,在i=j(最後一次應用)時,與(13)式(除後者還要做一次除法外)是完全一樣的,這兩種情況要求和的上線都是k=j-1(=i-1)。這意味著,不必費心去考慮對角線元素βjj是否會正落在對角線上,也不必考慮該列中,它下面的某個元素(未做除法的)αij,i=j+1,j+2,...,N-1是否會提升成為對角線元素β。
l         它首先找到每行的最大元素,而後(在找最大主元時)乘以一個比例系數,這就實現了「隱式主元法」。
 
運行示例:
Origin coefficient matrix:
 | 0.0 2.0 0.0 1.0 |
 | 2.0 2.0 3.0 2.0 |
 | 4.0 -3.0 0.0 1.0 |
 | 6.0 1.0 -6.0 -5.0 |
-----------------------------------------------
LU mixed matrix:
 | 6.0 1.0 -6.0 -5.0 |
 | 0.0 2.0 0.0 1.0 |
 | 0.3333333333333333 0.8333333333333334 5.0 2.833333333333333 |
 | 0.6666666666666666 -1.8333333333333333 0.8 3.8999999999999995 |
-----------------------------------------------
Origin left-hand vector b:
 | 0.0 |
 | -2.0 |
 | -7.0 |
 | 6.0 |
-----------------------------------------------
Final solution vector:
 | -0.5000000000000003 |
 | 1.0000000000000002 |
 | 0.33333333333333337 |
 | -2.0000000000000004 |
-----------------------------------------------
 
示常式序:
package com.nc4nr.chapter02.lu;

public class LU ...{

    // 4 * 4 coefficient matrix a
    double[][] a = ...{
            ...{0.0, 2.0, 0.0, 1.0},
            ...{2.0, 2.0, 3.0, 2.0},
            ...{4.0, -3.0, 0.0, 1.0},
            ...{6.0, 1.0, -6.0, -5.0}
    };
   
    // 4 * 1 coefficient matrix b
    double[] b = ...{
            0.0,
            -2.0,
            -7.0,
            6.0
    };
   
    int anrow = 4;
    int[] indx = new int[anrow];
    int parity = 1;
   
    private void lucmp() ...{
        final double tiny = 1.0e-20;
        int imax = 0, n = anrow;
        double big, m, sum, temp;
        double[] vv = new double[n];
       
        System.out.println("Origin coefficient matrix:");
        output(a,4);
       
        for (int i = 0; i < n; i++) ...{
            big = 0.0;
            for (int j = 0; j < n; j++) ...{
                if ((temp = Math.abs(a[i][j])) > big) big = temp;
            }
            if (big == 0.0) System.out.println("lu: singular matrix in lucmp.");
            vv[i] = 1.0/big;
        }
       
        for (int j = 0; j < n; j++) ...{
            for (int i = 0; i < j; i++) ...{
                sum = a[i][j];
                for (int k = 0; k < i; k++) sum -= a[i][k]*a[k][j];
                a[i][j] = sum;
            }
            big = 0.0;
            for (int i = j; i < n; i++) ...{
                sum = a[i][j];
                for (int k = 0; k < j; k++)    sum -= a[i][k]*a[k][j];
                a[i][j] = sum;
                if ((m = vv[i]*Math.abs(sum)) >= big) ...{
                    big = m;
                    imax = i;
                }
            }
            if (j != imax) ...{
                for(int k = 0; k < n; k++) ...{
                    m = a[imax][k];
                    a[imax][k] = a[j][k];
                    a[j][k] = m;
                }
                parity = -parity;
                m = vv[imax];
                vv[imax] = vv[j];
                vv[j] = m;
            }
            indx[j] = imax;
            if (a[j][j] == 0.0) a[j][j] = tiny;
            if (j != n - 1) ...{
                m = 1.0/a[j][j];
                for (int i = j+1; i < n; i++) a[i][j] *= m;
            }
        }
       
        System.out.println("LU mixed matrix:");
        output(a,4);
    }
   
    private void lubksb() ...{
        double sum;
        int n = anrow, ii = 0;
       
        System.out.println("Origin left-hand vector b:");
        output(b,4);
       
        for (int i = 0; i < n; i++) ...{
            int ip = indx[i];
            sum = b[ip];
            b[ip] = b[i];
            if (ii != 0)
                for (int j = ii - 1; j < i; j++) sum -= a[i][j]*b[j];
            else if (sum != 0.0)
                ii = i + 1;
            b[i] = sum;
        }
        for (int i = n-1; i >= 0; i--) ...{
            sum = b[i];
            for(int j = i + 1; j < n; j++) sum -= a[i][j]*b[j];
            b[i] = sum / a[i][i];
        }
        System.out.println("Final solution vector:");
        output(b,4);
    }
   
    private void output(double a[][], int anrow) ...{
        for (int i = 0; i < anrow; i++) ...{
            System.out.println(" | " + a[i][0] + " " +
                    a[i][1] + " " +
                    a[i][2] + " " +
                    a[i][3] + " | ");
        }
        System.out.println("-----------------------------------------------");
    }
   
    private void output(double[] b, int bnrow) ...{
        for (int i = 0; i < bnrow; i++) ...{
            System.out.println(" | " + b[i] + " | ");
        }
        System.out.println("-----------------------------------------------");
    }
   
    public LU() ...{
        lucmp(); // 分解
        lubksb(); // 回代
    }
   
    public static void main(String[] args) ...{
        new LU();
    }

}

⑥ 急球c語言!矩陣直接分解法(lu分解法)

以四維為例。系數在M中,常數項在N中
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
void main()
{
int i,j,m,n;
double TM=0,TMm=0,TN=0,TNn=0;
int NN = 4;
double L[4][4]= {{1,0,0,0}, {0,1,0,0}, {0,0,1,0},{0,0,0,1}, };
// M coef 4x4, N const 4x1
double M[4][4]={4,-1,0,2,-1,4,-1,0,0,-1,4,-1,2,0,-1,4};
double N[4][1]={-1,-7,9,0};

for(i=1;i<NN;i++) {
for(j=0;j<i;j++) {
for(m=0;m<j;m++) {
TM=M[i][m] * L[j][m] + TMm;
TMm=TM;
TM=0;
};
M[i][j]=M[i][j]-TMm;
TMm=0;
L[i][j]=M[i][j]/M[j][j];
};
for(n=0;n<i;n++) {
TN=M[i][n]*L[i][n]+TNn;
TNn=TN;
TN=0;
}
M[i][i]=M[i][i]-TNn;
TNn=0;
};

// solve
for(i=0;i<NN;i++) {
for(m=0;m<i;m++) {
TM=L[i][m]*N[m][0]+TMm;
TMm=TM;
TM=0;
};
N[i][0]=N[i][0]-TMm;
TMm=0;
};
for(i=0;i<NN;i++) {
N[i][0]=N[i][0]/M[i][i];
};
for(i=NN-1;i>=0;i--) {
for(m=i+1;m<NN;m++) {
TM=L[m][i]*N[m][0]+TMm;
TMm=TM;
TM=0;
};
N[i][0]=N[i][0]-TMm;
TMm=0;
};
for(i=0;i<NN;i++) printf("%0.2f\n",N[i][0]);
}

⑦ 求C語言課程設計:用高斯列主元消元法解線性方程組

這里向你推薦一下克魯特演算法(其實就是對高斯列主元消元法進行優化,使之更適合於計算機編程),首先將矩陣A進行LU分解(將系數矩陣分解成一個上三角矩陣和一個下三角矩陣),分解的過程中用到了隱式的主元尋找法,同時利用克魯特演算法可以將兩個n*n矩陣壓縮到一個n*n矩陣中,大大節省了存儲空間提高了計算速度。
方程可化為L*U*x=B,令U*x=y --->L*y=B
然後利用回代先求y,再利用y求x
因為該方法在求解過程中不涉及增廣矩陣所以矩陣B幾乎不參與什麼運算,所以它的計算速度應該能夠達到高斯列主元消元法的三倍,但原理與其基本一致。
而且我在程序中使用了動態數組方便你今後進行擴展。

以下程序按照《矩陣論第二版》和《C語言數值計演算法方法大全》編寫,LU分解部分程序主要參考了《C語言數值計演算法方法大全》第二章的程序
如果你需要詳細的理論講解我可以將這兩本書和源程序發給你,上面的論述相當詳細足夠你答辯用的了,我的郵箱[email protected]

計算結果:
A矩陣:
2 2 5
3 4 7
1 3 3
B矩陣:
5
6
5
解矩陣:
x 1=-7
x 2=0.333333
x 3=3.66667

Press any key to continue

#include <cmath>
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <functional>
#include <vector>
#include <algorithm>
using namespace std;
#define TINY 1.0e-20 //A small number.
#define N 3
void ludcmp(vector<vector<float> > &a, int n, vector<int> &indx, float &d);//對矩陣進行LU分解
void lubksb(vector<vector<float> > &a, int n, vector<int> &indx, vector<float> &b);//對矩陣進行前向和後向回代
void root(vector<vector<float> > &x,vector<float> &col);//解方程結果保存在y中
void iniv(vector<vector<float> > &x,vector<float> line,int n);//對二維動態數組進行初始化
void main()
{
int i,j,n=N;//輸入矩陣的維數
float A[N][N]={{2,2,5},{3,4,7},{1,3,3}};//左邊A矩陣
float B[N]={5,6,5};//右邊B矩陣
vector<vector<float> > x;//建立動態二維數組存放A,保證你的程序進行擴展時只改A,B,N
vector<float> line;
vector<float> y(n);//建立動態數組存放B
iniv(x,line,n);
y.clear();
for(i=0;i<n;i++)//將A賦給x,B賦給y
{
y.push_back(B[i]);
for(j=0;j<n;j++)
{
x[i].push_back(A[i][j]);
}
}
cout<<"A矩陣:"<<endl;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<x[i][j]<<" ";
}
cout<<endl;
}
cout<<"B矩陣:"<<endl;
for(i=0;i<n;i++)
{
cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<y[i]<<endl;
}
root(x,y);//求根
cout<<"解矩陣:"<<endl;
for(i=0;i<n;i++)
{
cout<<setw(2)<<setiosflags(ios::left)<<"x"<<i+1<<"="<<y[i]<<endl;
}
cout<<endl;
}
void root(vector<vector<float> > &x,vector<float> &col)
{
int n=x.size(),i=0,j=0;
vector<int> index(n);//用於記錄尋找主元素過程中對矩陣的初等變換
index.clear();
float m=1.0;//記錄變換方式,此程序中無用
ludcmp(x,n,index,m);//進行LU分解
lubksb(x,n,index,col);//根據分解結果進行回帶
}
//以下程序按照《矩陣論第二版》和《C語言數值計演算法方法大全》編寫,LU分解部分程序主要參考了《C語言數值計演算法方法大全》第二章的程序
//如果你需要詳細的理論講解我可以將這兩本書和源程序發給你,我的郵箱[email protected]
void ludcmp(vector<vector<float> > &a, int n, vector<int> &indx, float &d)
{

int i,imax,j,k;
float big=0,m=0,sum=0,temp=0;
vector<float> vv(n);
vv.clear();
d=1.0;
for (i=0;i<n;i++)
{
big=0.0;
for (j=0;j<n;j++)
if ((temp=fabs(a[i][j])) > big)
big=temp;
vv[i]=1.0/big;
}
for (j=0;j<n;j++)
{
for (i=0;i<j;i++)
{
sum=a[i][j];
for (k=0;k<i;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
}
big=0.0;
for (i=j;i<n;i++)
{
sum=a[i][j];
for (k=0;k<j;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
if ( (m=vv[i]*fabs(sum)) >= big)
{
big=m;
imax=i;
}
}
if (j != imax)
{
for (k=0;k<n;k++)
{
m=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=m;
}
d = -(d);
vv[imax]=vv[j];
}
indx[j]=imax;
if (a[j][j] == 0.0)
a[j][j]=TINY;
if (j != n)
{
m=1.0/(a[j][j]);
for (i=j+1;i<n;i++)
a[i][j] *= m;
}
}
}
void lubksb(vector<vector<float> > &a, int n, vector<int> &indx, vector<float> &b)
{
int i,ii=0,ip,j;
float sum;
for(i=0;i<n;i++)//按LU分解時尋找主元所進行的初等變換進行反邊變換。
{
ip=indx[i];
sum=b[ip];
b[ip]=b[i];
b[i]=sum;
}
sum=0;
for (i=1;i<n;i++)
{
sum=0;
for(j=0;j<i;j++)
{
sum+=a[i][j]*b[j];
}
b[i]=b[i]-sum;
}
b[n-1]=b[n-1]/a[n-1][n-1];
for (i=n-2;i>=0;i--)
{
sum=0;
for(j=i+1;j<n;j++)
{
sum+=a[i][j]*b[j];
}
b[i]=(b[i]-sum)/a[i][i];
}
}
void iniv(vector<vector<float> > &x,vector<float> line,int n)
{
int i,j;
for(i=0;i<n;i++)
{
x.push_back(line);
for(j=0;j<n;j++)
{
x[i].clear();
}
}
}

熱點內容
c封裝dll文件夾 發布:2024-12-21 20:01:01 瀏覽:206
活塞空氣壓縮機 發布:2024-12-21 19:59:17 瀏覽:31
勤哲伺服器如何用WPS 發布:2024-12-21 19:52:55 瀏覽:320
c語言b是什麼意思 發布:2024-12-21 19:52:51 瀏覽:916
人渣需要什麼樣的配置 發布:2024-12-21 19:52:44 瀏覽:558
unity腳本製作 發布:2024-12-21 19:52:02 瀏覽:562
伺服器和雲伺服器區別 發布:2024-12-21 19:26:38 瀏覽:650
linuxhttp文件伺服器搭建 發布:2024-12-21 19:26:35 瀏覽:676
android與c通信 發布:2024-12-21 19:16:37 瀏覽:334
電子秤腳本 發布:2024-12-21 19:12:33 瀏覽:276