Doolitter分解 三对角矩阵分解 拟三对角分解

最新推荐文章于 2022-12-08 17:25:22 发布
最新推荐文章于 2022-12-08 17:25:22 发布
阅读量1.2k 收藏
点赞数
文章标签： matlab 
  1 #include <cstdio>
  2 #include <cstdlib>
  3 #include <algorithm>
  4 #include <cmath>
  5 #include <cassert>
  6 #include <vector>
  7 #include <ctime>
  8 
  9 
 10 class MclVector
 11 {
 12 public:
 13     int n;
 14     double *Mat;
 15     /**
 16       type=0: 列向量(默认)
 17       type=1: 行向量
 18     **/
 19     int type;
 20 
 21     MclVector() { Mat=NULL; n=0; }
 22     MclVector(int len,double initVal=0.0)
 23     {
 24         n=len;
 25         Mat=new double[n+1];
 26         for(int i=0;i<=n;i++) Mat[i]=initVal;
 27         type=0;
 28     }
 29 
 30 
 31     double operator[](int id) const
 32     {
 33         return Mat[id];
 34     }
 35 
 36 
 37     double& operator[](int id)
 38     {
 39         return Mat[id];
 40     }
 41 
 42     double length() const
 43     {
 44         double sum=0;
 45         for(int i=1;i<=n;i++) sum+=Mat[i]*Mat[i];
 46         return sqrt(sum);
 47     }
 48 
 49     MclVector operator*(double val) const
 50     {
 51         MclVector ans=MclVector(n);
 52         for(int i=1;i<=n;i++) ans[i]=Mat[i]*val;
 53         return ans;
 54     }
 55 
 56     MclVector operator/(double val) const
 57     {
 58         MclVector ans=MclVector(n);
 59         for(int i=1;i<=n;i++) ans[i]=Mat[i]/val;
 60         return ans;
 61     }
 62 
 63     MclVector operator+(const MclVector &newVector) const
 64     {
 65         MclVector ans=MclVector(n);
 66         for(int i=1;i<=n;i++) ans[i]=Mat[i]+newVector[i];
 67         return ans;
 68     }
 69 
 70 
 71     MclVector operator-(const MclVector &newVector) const
 72     {
 73         MclVector ans=MclVector(n);
 74         for(int i=1;i<=n;i++) ans[i]=Mat[i]-newVector[i];
 75         return ans;
 76     }
 77 
 78     MclVector operator*=(double val)
 79     {
 80         for(int i=1;i<=n;i++) Mat[i]=Mat[i]*val;
 81         return *this;
 82     }
 83 
 84     MclVector operator/=(double val)
 85     {
 86         for(int i=1;i<=n;i++) Mat[i]=Mat[i]/val;
 87         return *this;
 88     }
 89 
 90 
 91     MclVector operator+=(const MclVector &newVector)
 92     {
 93         for(int i=1;i<=n;i++) Mat[i]+=newVector[i];
 94         return *this;
 95     }
 96 
 97     MclVector operator-=(const MclVector &newVector)
 98     {
 99         for(int i=1;i<=n;i++) Mat[i]-=newVector[i];
100         return *this;
101     }
102 
103     MclVector GetTranspose() const
104     {
105         MclVector ans=*this;
106         ans.type=1;
107         return ans;
108     }
109 
110 
111     void print() const
112     {
113         for(int i=1;i<=n;i++) printf("%8.3lf ",Mat[i]);
114         puts("");
115     }
116 
117 
118 };
119 
120 class MclMatrix
121 {
122 public:
123     int row,col;
124     MclVector *Mat;
125 
126     MclMatrix() {Mat=NULL;}
127     MclMatrix(int _row,int _col,double initVal=0.0)
128     {
129         row=_row;
130         col=_col;
131         Mat=new MclVector[row+1];
132         for(int i=0;i<=row;i++) Mat[i]=MclVector(col,initVal);
133     }
134 
135     void setIdentityMatrix()
136     {
137         for(int i=1;i<=row;i++)
138         {
139             for(int j=1;j<=col;j++)
140             {
141                 if(i==j) Mat[i][j]=1;
142                 else Mat[i][j]=0;
143             }
144         }
145     }
146 
147     MclMatrix GetTranspose() const
148     {
149         MclMatrix ans=MclMatrix(col,row);
150         for(int i=1;i<=ans.row;i++)
151         {
152             for(int j=1;j<=ans.col;j++)
153             {
154                 ans[i][j]=Mat[j][i];
155             }
156         }
157         return ans;
158     }
159 
160     void print() const
161     {
162         for(int i=1;i<=row;i++) Mat[i].print();
163         puts("");
164     }
165 
166     MclVector& operator[](int id) const
167     {
168         return Mat[id];
169     }
170 
171 
172     MclVector& operator[](int id)
173     {
174         return Mat[id];
175     }
176 
177     MclMatrix operator*(const MclMatrix &Matrix) const
178     {
179         MclMatrix ans=MclMatrix(row,Matrix.col);
180         for(int i=1;i<=row;i++)
181         {
182             for(int j=1;j<=Matrix.col;j++)
183             {
184                 for(int k=1;k<=col;k++)
185                 {
186                     ans[i][j]+=Mat[i][k]*Matrix[k][j];
187                 }
188             }
189         }
190         return ans;
191     }
192 
193     MclMatrix operator+(const MclMatrix &Matrix) const
194     {
195         MclMatrix ans=MclMatrix(row,Matrix.col);
196         for(int i=1;i<=row;i++)
197         {
198             for(int j=1;j<=Matrix.col;j++)
199             {
200                 ans[i][j]=Mat[i][j]+Matrix[i][j];
201             }
202         }
203         return ans;
204     }
205 
206     MclMatrix operator-(const MclMatrix &Matrix) const
207     {
208         MclMatrix ans=MclMatrix(row,Matrix.col);
209         for(int i=1;i<=row;i++)
210         {
211             for(int j=1;j<=Matrix.col;j++)
212             {
213                 ans[i][j]=Mat[i][j]-Matrix[i][j];
214             }
215         }
216         return ans;
217     }
218 
219     MclVector GetCol(int colId) const
220     {
221         MclVector ans=MclVector(row);
222         for(int i=1;i<=row;i++) ans[i]=Mat[i][colId];
223         return ans;
224     }
225     MclVector GetRow(int rowId) const
226     {
227         MclVector ans=MclVector(row);
228         for(int i=1;i<=col;i++) ans[i]=Mat[rowId][i];
229         return ans;
230     }
231 
232 
233     MclMatrix operator*=(const MclMatrix &Matrix)
234     {
235         return *this=*this*Matrix;
236     }
237     MclMatrix operator+=(const MclMatrix &Matrix)
238     {
239         return *this=*this+Matrix;
240     }
241     MclMatrix operator-=(const MclMatrix &Matrix)
242     {
243         return *this=*this-Matrix;
244     }
245 
246     MclMatrix operator*(double x) const
247     {
248         MclMatrix ans=*this;
249         for(int i=1;i<=row;i++)
250         {
251             for(int j=1;j<=col;j++)
252             {
253                 ans[i][j]*=x;
254             }
255         }
256         return ans;
257     }
258 
259 };
260 
261 MclMatrix vectorMulVector(const MclVector &A,const MclVector& B)
262 {
263     if(A.type==0)
264     {
265         MclMatrix ans=MclMatrix(A.n,B.n);
266         for(int i=1;i<=A.n;i++)
267         {
268             for(int j=1;j<=B.n;j++)
269             {
270                 ans[i][j]+=A[i]*B[j];
271             }
272         }
273         return ans;
274     }
275     else
276     {
277         assert(A.n==B.n);
278         MclMatrix ans=MclMatrix(1,1);
279         for(int i=1;i<=A.n;i++)
280         {
281             ans[1][1]+=A[i]*B[i];
282         }
283         return ans;
284     }
285 }
286 
287 int sgn(double x)
288 {
289     if(x<-0.000001) return -1;
290     if(x>0.000001) return 1;
291     return 0;
292 }
293 
294 
295 
296 /**
297    矩阵的 Doolittle分解:
298       [1] Mat是方阵
299       [2] Mat的前n-1阶主子式行列式不为0
300       [3] 分解的L为单位下三角阵
301       [4] 分解的U为上三角阵
302       [5] 返回值为<L,R>
303 **/
304 std::pair<MclMatrix,MclMatrix> DoolittleSplit(const MclMatrix &Mat)
305 {
306     int n=Mat.row;
307     MclMatrix L=MclMatrix(n,n);
308     MclMatrix U=MclMatrix(n,n);
309     for(int k=1;k<=n;k++)
310     {
311         for(int j=k;j<=n;j++)
312         {
313             U[k][j]=Mat[k][j];
314             for(int t=1;t<=k-1;t++) U[k][j]-=L[k][t]*U[t][j];
315         }
316         if(k==n) continue;
317 
318         for(int i=k+1;i<=n;i++)
319         {
320             L[i][k]=Mat[i][k];
321             for(int t=1;t<=k-1;t++) L[i][k]-=L[i][t]*U[t][k];
322             L[i][k]/=U[k][k];
323         }
324     }
325     for(int i=1;i<=n;i++) L[i][i]=1;
326     return std::make_pair(L,U);
327 }
328 
329 
330 
331 /**
332    三角矩阵分解:
333       [1] Mat是方阵
334       [2] j<i且i-j>r时 Mat[i][j]=0
335       [2] j>i且j-i>s时 Mat[i][j]=0
336 **/
337 std::pair<MclMatrix,MclMatrix> TriangleSplit(const MclMatrix &Mat,int r,int s)
338 {
339     int n=Mat.row;
340     MclMatrix L=MclMatrix(n,n);
341     MclMatrix U=MclMatrix(n,n);
342     for(int k=1;k<=n;k++)
343     {
344         for(int j=k;j<=n;j++)
345         {
346             U[k][j]=Mat[k][j];
347             for(int t=std::max(1,std::max(k-r,j-s));t<=k-1;t++) U[k][j]-=L[k][t]*U[t][j];
348         }
349         if(k==n) continue;
350 
351         for(int i=k+1;i<=n;i++)
352         {
353             L[i][k]=Mat[i][k];
354             for(int t=std::max(1,std::max(i-r,k-s));t<=k-1;t++) L[i][k]-=L[i][t]*U[t][k];
355             L[i][k]/=U[k][k];
356         }
357     }
358     for(int i=1;i<=n;i++) L[i][i]=1;
359     return std::make_pair(L,U);
360 }
361 
362 /**
363   拟三对角矩阵分解
364   对n=5 矩阵A样子如下：
365       a1  c1  0   0  d1
366       d2  a2  c2  0  0
367       0   d3  a3  c3 0
368       0   0   d4  a4 c4
369       c5  0   0   d5 a5
370    即输入为三个长度为n的向量
371 
372    A=LU
373    L样子如下：
374      p1  0  0  0  0
375      d2  p2 0  0  0
376      0   d3 p3 0  0
377      0   0  d3 p4 0
378      r1  r2 r3 r4 r5
379 
380    U样子如下:
381      1  q1  0  0  s1
382      0  1   q2 0  s2
383      0  0   1  q3 s3
384      0  0   0  1  s4
385      0  0   0  0  1
386 
387     即将返回p,q,s,r四个向量(所有的向量长度都是n)
388     vector[0]=p
389     vector[1]=q
390     vector[2]=s
391     vector[3]=r
392 **/
393 
394 std::vector<MclVector> QuasiDiagonalSplit(const MclVector &a,const MclVector &c,const MclVector &d)
395 {
396     int n=a.n;
397     assert(c.n==n);
398     assert(d.n==n);
399     assert(n>0);
400 
401     MclVector p=MclVector(n);
402     MclVector q=MclVector(n);
403     MclVector s=MclVector(n);
404     MclVector r=MclVector(n);
405 
406     p[1]=a[1];
407     for(int i=1;i<=n-2;i++)
408     {
409         q[i]=c[i]/p[i];
410         p[i+1]=a[i+1]-d[i+1]*q[i];
411     }
412 
413     s[1]=d[1]/p[1];
414     for(int i=2;i<=n-2;i++) s[i]=-d[i]*s[i-1]/p[i];
415     s[n-1]=(c[n-1]-d[n-1]*s[n-2])/p[n-1];
416 
417     r[1]=c[n];
418     for(int j=2;j<=n-2;j++) r[j]=-r[j-1]*q[j-1];
419     r[n-1]=d[n]-r[n-2]*q[n-2];
420     r[n]=a[n];
421     for(int j=1;j<=n-1;j++) r[n]=r[n]-r[j]*s[j];
422 
423     std::vector<MclVector> ans;
424     ans.push_back(p);
425     ans.push_back(q);
426     ans.push_back(s);
427     ans.push_back(r);
428 
429     return ans;
430 }

 

n阶三取对角阵的LU分解及利用其解方程组
12-24
高校计算方法上机作业之对n阶三对角阵的LU分解及利用其解方程组
(追赶法求三对角矩阵、LU分解)
Giannis_34的博客
08-08 1万+
(追赶法求三对角矩阵、LU分解) 应用场景 数模中的微分方程问题中偶尔回涉及到使用向后差分法求解问题，不论是六点法向后差分还是四点法向后差分都会涉及到求解三对角矩阵的情况。 以四点法向后差分为例 四点法得到的代数矩阵为三对角线性方程组AU=b，一般可以用matlab求逆命令求解U，但因本题中系数矩阵规模较大，且主对角线和两个次对角线之外，其余元素均为0，因而求逆解方程组会导致浪费内存、程序运行缓慢。因此本文利用追赶法求解三对角线性方程组。三对角线性方程组的追赶法求解过程如下图。 计算过程 三对角矩阵的形式如
对称矩阵的三对角分解(Lanzos分解算法）-MINRES算法预热
最新发布
forrestguang的博客
12-08 1721
首先介绍Lanczos分解，Lanzos把对称矩阵转换为一个三对角对称矩阵。考虑三对角对称矩阵如下，考虑正交分解 T=QTAQT = Q^T A QT=QTAQ T=(α1β10⋯00β1α2β20⋯00β2α3β3⋯00⋯⋯⋯⋯00⋯0βn−2αn−1βn−100⋯0βn−1αn)T=\left(\begin{array}{cccccc} \alpha_1 & \beta_1 & 0 & \cdots & 0 & 0\\ \beta_1 & \alpha_2 & \beta_2 & 0 & \cdots
矩阵分解
phily123的博客
11-27 402
矩阵分解的方法有：三角分解法（LU分解）、OR分解和奇异分解法。 LU分解：讲矩阵分解为一个下三角矩阵与一个上三角矩阵相乘，matlab实现：lu(A) QR分解：将矩阵分解为正定矩阵和上三角矩阵相乘的形式，matlab实现：qr(A) 奇异分解：matlab实现：svd(A) 矩阵的LU分解 ...
三对角矩阵lu分解matlab代码-Mech105-projectCode:这是Mech105（机械工程问题解决）的所有项目文件
06-16
三对角矩阵划分matlab代码Mech105-projectCode 这是Mech 105（机械工程问题解决）的所有项目文件 特殊矩阵 该程序将创建一个具有两个输入数字的特殊矩阵。 我操作了二维数组的各种单元格，同时过滤整个数组以访问以前的数组值。 例如，如果我的函数是在输入为n = 5和m = 4的情况下调用的，则将输出以下矩阵： 1 2 3 4 2 4 7 11 3 7 14 25 4 11 25 50 5 16 41 91 从示例中您可以看到我的函数创建了一个 5 行 4 列的数组。 第一个元素从一个开始。 对于数组的其余部分，元素是通过将其上方的元素与其左侧的元素相加来定义的。 ND 此文件计算输入的月数和天数时的天数。 该程序能够解释闰年。 如果计算闰年的天数，程序将在最终答案中添加 1 天以说明闰年。 用户必须输入是否为闰年。 例如，如果我的函数在输入为mo = 2 、 da = 1 、 leap = 0 的情况下被调用，则结果将输出为 60 天。 作业 7 问题 1 此代码显示了矩阵代数中使用的左除法。 基本公式如下： A\B ; 作业 15 此代码显示 LU 分解。 作
三对角矩阵_一个挂角的三对角矩阵_
10-02
对于三对角矩阵，尤其是带挂角元素的情况，可以考虑使用高斯消元法、LU分解、或者更高效的Thomas算法（适用于纯三对角矩阵）。 高斯消元法是线性代数的基础，通过一系列行变换将系数矩阵转化为上三角形或下三角形，...
LDL矩阵分解：将矩阵分解为下三角矩阵L和对角矩阵D-matlab开发
06-01
MATLAB 包括 LDL 分解，但是它返回块对角矩阵 D 而不是对角矩阵 D。这些实现返回真正的对角矩阵 D。 该软件包包括两个 LDL 实现。 一种实现采用对称矩阵 A 并返回 [L, D] ： L*D*L' = ldl(A)。 每当您在 Gramian ...
SOR 方法：输入一个方阵。 将矩阵分解为对角矩阵、下三角矩阵和上三角矩阵。-matlab开发
05-29
%将方阵A分解为三个矩阵：对角矩阵(D)； 严格下三角矩阵（L）； 严格上三角矩阵(U) D = 诊断（诊断（A））； L =-tril(A,-1); U = -triu(A,1); a = (D-欧米茄*L); 对于 i=1:N x = a\(((1-omega)*D + omega*U)*x_0) ...
基于 Lapack 接口的 Hermitian 或对称矩阵的三对角化：代码通过调用 Lapack 例程计算 Hermitian 矩阵的三对角分解。-matlab开发
06-01
此代码基于 Lapack 子程序 ZHETRD 和 ZUNGTR 提供了可靠的三对角矩阵分解程序。 Lapack 是数值分析中最可靠的 Fortran 程序之一。 在 Lapack 中，三对角化用于计算 Hermitian 矩阵的特征值分解。 请参阅...
基于Lapack接口的矩阵双对角化：代码通过调用Lapack例程计算矩阵的双对角矩阵分解。-matlab开发
06-01
此代码提供了基于 Lapack 子例程的可靠双对角分解例程。 Lapack 是数值分析中最可靠的 Fortran 程序之一。 在 Lapack 中，双对角化用于通过例程 zgebrd 和 zungbr 计算矩阵的 SVD。 请参阅...
三对角阵特征解法
05-25
三对角矩阵的特征值最佳解法，讲解详细，通俗易懂。
深度学习花书 笔记1-矩阵对角分解、矩阵可逆与最小二乘
ruotianxia的博客
04-03 244
深度学习花书 笔记1-矩阵对角分解、矩阵可逆与最小二乘
矩阵对角化，SVD分解
Ford程的博客
09-20 5823
- [矩阵对角化](#矩阵对角化) - [SVD分解](#svd分解) - [参考链接](#参考链接) 矩阵对角化 矩阵的相似 设 A\boldsymbol{A}A、 B\boldsymbol{B}B 为两个nnn阶矩阵，若存在可逆矩阵 P\boldsymbol{P}P，使得 P−1AP=B \boldsymbol{P}^{^{-1}}\boldsymbol{A}\boldsymbol{P...
矩阵对角化和SVD分解
J的博客
03-03 461
矩阵对角化的概念与意义： SVD： SVD还可以极大减少算法的复杂度：
【线性代数】矩阵的特征值分解（对角化、谱分解）
热门推荐
sxl的博客
06-12 1万+
   矩阵的特征分解又可以称作矩阵的对角化、谱分解。其在机器学习和图机器学习中有非常广泛的应用。本节主要介绍矩阵的特征分解的解法，意义，实际应用。除此之外，矩阵的特征分解与矩阵的特征值和特征向量有关联，之前在【线性代数】理解特征值和特征向量文章中，对于特征值和特征向量的一些相关概念没有涉及到的，会在此进行补充。   内容为自己的学习总结，其中多有借鉴他人的地方，最后一并给出链接。  特征向量所在直线上的向量都是特征向量。 这样的直线可能不止一条。 我们称矩阵A所有特征向量对应的直线的集合为特征空间。  对角
线性代数笔记(8) 各种矩阵分解：对角、三角、满秩、QR、奇异值分解
Cake_C的博客
01-04 3071
已经考完了，考的很开心，下次还考（） 记录一下各种分解的做法： 三角分解 LDR又称为Cholesky分解 LR分解就直接设分解下来的下、上三角矩阵各未知数，根据矩阵运算解方程即可。最简单了。 对角分解 正规矩阵总能对角分解： 分解流程： 满秩分解 最后一步其实就相当于P^-1取前r列，Q^-1取前r行。其中r为A的秩。 行最简型就是各行第一个为1，且所在列其他行为0，列同理。 QR分解 QR分解就是分解成酉矩阵（单位正交基）*上三角矩阵的形式 例题： 奇异值分解 PQ均为酉矩阵。 下
【Matlab】Matlab实现方阵LU分解——Doolittle方法
fengyuanshang的博客
03-19 8483
1、消元过程：求解L/U三角方阵，并在原方阵A中存储 2、回代过程 A = [5 7 9 10; 6 8 10 9; 7 10 8 7; 5 7 6 5]; b = [1;1; 1; 1]; xc = A\b n = 4 %求解三角方阵 for r = 1 : n for j = r : n if r &gt; 1; ...
三对角阵的LU分解和三对角方程组的求解（C语言）
03-30 7303
/*三对角阵的LU分解和三对角方程组的求解 -------------A=LU的分解算法------- 参考教材：《数值分析》李乃成，梅立泉，科学出版社     《计算方法教程》第二版 凌永祥，陈明逵 */ #include #include #include int main(void) {  int i,j,n;  int N;  printf("请输入 N(10,
三对角、五对角追赶法求解线性方程组
wxkhturfun的博客
05-28 6183
前言 谨以此文纪念2022-5-27 日19：00-21：30学业生涯最后一次考试《数值分析》 参考《五对角线性方程组追赶法》——王礼广等；以及《求解五对角方程组追赶法》——夏爱生等；以及《数值分析》课本第二版 一.三对角 1.1问题引入 问题引入： A={613861...........................86186}，A= \begin{Bmatrix}6& 1 & 3 &\\ 8 & 6 & 1& \\ & ... &
MATLAB解决挂角三对角矩阵求解技巧
对于标准的三对角矩阵，可以使用`tridiag`函数创建一个三对角矩阵，然后使用矩阵分解方法，如LU分解，或者特殊的三对角矩阵求解器，例如`triql`函数（MATLAB 6.5版本之前）或`mldivide` (`\`) 操作符（对角矩阵可求...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值