数值分析常见基本算法及MATLAB代码总结

2023-10-14 04:55| 来源: 网络整理| 查看: 265

刚考完研究生的数值分析课，打算整理一下平时的上机代码，做一个汇总。

数值分析大致八部分内容——

解线性方程组的直接法：Gauss消去法与矩阵三角分解法（Doolittle分解法相比Crout分解法更常用）及其选择列主元的改进方法、Doolittle分解法的延伸（实对称正定矩阵利用Cholesky分解得到的平方根法、三对角矩阵作为线性方程组系数矩阵的追赶法）

解线性方程组的迭代法：Jacobi迭代法、Gauss-Seidel迭代法（利用前者每次迭代已得到的最新分量加速）、逐次超松弛（SOR，Successive Over-Relaxation）方法（取适当松弛因子ω可以比Gauss-Seidel方法收敛得更快）以及共轭梯度法（教材未讲，前提条件：系数矩阵A为对称正定矩阵）

函数拟合的插值法：拉格朗日（Lagrange）插值法与牛顿（Newton）插值法（根据插值多项式的存在唯一性定理，对同一插值问题两者展开的结果相同。牛顿插值法的优势在于：每次增加或删减一个节点，只需插值多项式增加或减少一项，相当于对不同插值节点数的情形给出了迭代公式）、分段插值法（克服插值节点过多时逼近效果反而很差的Runge现象，常见的是分段线性插值，也包括分段抛物线插值法）、Hermite插值法（插值条件包括所有已知节点的函数值和导数值，即两种条件的等式数量相等）、样条插值法（最基本的是三次样条插值法，也属于分段插值法的一种，整体具有二阶连续微商）。需要指出的是，除了分段插值法和样条插值法，其他几种方法经常是交叉混合使用的，因为实际插值条件中节点函数值和节点导数值条件的数目往往不同，可以灵活选择部分条件先构造Lagrange或Hermite插值多项式，然后再添加一项，使得满足附加插值条件的同时，已有插值条件不受影响。两项相加并解出待定系数就是完整的一般形式的插值多项式。此时插值余项也需自行推导，但一般是反复运用罗尔定理得出的带f(ξ)的高阶导数的表达式。

函数近似替代（复杂化简单、高次化低次）的数值逼近方法：数值逼近中引入了函数范数和函数内积的概念。前者用来度量逼近函数与原函数在一个区间内的整体误差，后者广泛用于各种数值逼近方法的计算过程中。函数的∞-范数对应最佳一致逼近；函数的2-范数（Euclid-范数）对应最佳平方逼近。这两种是最基本的数值逼近方法。（上课只选讲了最佳平方逼近，按老师原话，对工科学生来说最佳一致逼近理论性太强，且实际中用得不是很多）。将函数逼近中的线性无关函数族（类比线性代数中n维向量空间的向量线性表示选择的一组基）进行正交化（类似于线性代数中的施密特正交化），取适当的积分区间和权函数，可以衍生出最佳一致逼近中的切比雪夫（Chebyshev）多项式和最佳平方逼近中的勒让德（Legendre）多项式。以引入勒让德多项式为例，此时解法方程组，系数矩阵G为对角阵，可以大大降低函数逼近计算中解线性方程组的工作量。使用的教材中，将统计学中数据拟合的最小二乘法也归到数值逼近一章，个人理解是因为它的理论证明过程和计算套路都和最佳平方逼近及其相似，可以类比学习。

数值积分算法与数值微分（未学）：数值积分一章从机械求积公式出发，引入插值法并采用等距插值节点得到Newton-Cote公式。积分区间等分数n取不同整数有一系列不同数值求积公式，常见的为梯形公式（n=1)、Simpson公式（n=2）及Cotes公式（n=4）。当n≥8时，Cotes系数出现负数，数值不稳定，又有与Newton-Cotes公式对应的一系列复化求积公式（为保证数值积分精度，需积分步长相对固定，维持在一个较小的值。可以等分积分区间，在每个子区间中用n较小(＜8)的Newton-Cotes公式，这是复化求积公式的基本思想。）除此之外数值积分还有变步长的Romberg方法，适用于无法事先选择恰当的步长，在计算过程需要改变步长的情形。

解非线性方程及方程组的数值方法：对于非线性方程，最简单的方法是二分法，这个在高中就涉及过。常用的迭代法包括不动点迭代法、牛顿迭代法（Newton-Raphson迭代法或N-R法）和割线法。对于不动点迭代法，需将f(x)=0适当变形，取合适的φ（x）满足不动点迭代法的收敛条件（①"压缩镜像“；②迭代初值选取区间上导数值的绝对值恒不大于一个小于1的正数）；N-R法可以看做不动点迭代法φ（x）取一特殊函数的特例，也具有鲜明的几何意义，又称切线法；割线法是考虑到部分函数计算导数困难，用差商近似导数所做的改进。Newton法可以推广到非线性方程组的情形，在此基础上有一类拟Newton法，是针对前者的Jacobian矩阵中精确求解数目为方阵阶数n的平方的偏导数值的困难而做的改进。

矩阵特征值的数值解法：乘幂法与反幂法、HouseHolder方法、QR分解法（上课只讲了乘幂法，其他的以后用到了再学习。）

常微分方程的数值解法：欧拉（Euler）方法(前进与后退的Euler公式，综合两者的梯形公式，以及梯形公式中u_k+1项再代入前进的Euler公式，消除隐式参数，得到改进的欧拉公式)、龙格-库塔（Runge-Kutta）法、线性多步法。（P.S. 龙格-库塔法虽然是选学内容，上课没讲，但是对于电气、计算机、土木等工科研究生来说可能是一个重要的数学工具，给他个triple star 吧，今后用到了会细究。印象中本科的电器学教材中也有一节提到了这个）

数值分析课程中还有一些进阶模块，例如积分方程或积分-微分方程的数值解法，各种最优化方法等，也许超出了研究生教学大纲要求掌握的基本内容，今后如有用到再做补充。

下面逐一上代码，有些数值方法可能理论性较强，教材没有给出示例MATLAB程序或还没有MATLAB运行验证过，日后再更新。

1、选取列主元的高斯消去法

function x= pivot_Gauss(A,b) %功能：用Gauss列主元消去法解n阶线性方程组Ax=b n=length(b); for k=1:n-1 %步骤2：消元 %步骤2.1：选列主元 [max_value,max_index]=max(abs(A(k:n,k))); rk = k+max_index-1; %加上该列中前n-1个元素，还原主元在矩阵A中的行标 if max_value == 0 %步骤2.2：若A(rk,k)=0,则输出信息“A为奇异阵”，算法终止；否则执行步骤2.3 warning('系数矩阵奇异！'); return; end if rk~=k %步骤2.3：若rk!=k，则交换增广矩阵的第k行和第rk行 t=A(k,:);A(k,:)=A(rk,:);A(rk,:)=t; %交换矩阵的第k行与第rk行（从第k列开始即可，因为前k-1列下三角元素均已消零） t=b(k);b(k)=b(rk);b(rk)=t; end for i=k+1:n %步骤2.4：对i=k+1,...,n做以下两步：计算乘子l_ik;计算a_ij和b_i L(i,k)=A(i,k)/A(k,k); A(i,k+1:n)=A(i,k+1:n)-L(i,k)*A(k,k+1:n); b(i)=b(i)-L(i,k)*b(k); end end %步骤3：若a_nn=0，则输出信息“A为奇异阵！”，算法终止 if A(n,n)==0 warning('系数矩阵奇异！'); return; end %步骤4：回代求解 for k=n:-1:1 if k==n x(n)=b(n)/A(n,n); else x(k)=(b(k)-sum(A(k,k+1:n).*x(k+1:n)))/A(k,k); end end x

2、选取列主元的Doolittle分解法

% Doolittle.m function x=Doolittle2_3(A,b) %功能：用按列选主元的Doolittle方法求解n阶线性方程组Ax=b %假设A的各阶顺序主子式均不为零 n=length(A);L=eye(n);U=zeros(n); %length(A)返回矩阵A最大维数； %I = eye(n) 返回一个主对角线元素为 1 且其他位置元素为 0 的 n×n 单位矩阵。 %步骤1：输入系数矩阵A和及右端向量b %步骤2：对k=1,2,...,n-1执行步骤2.1~2.4 for k=1:n-1 %省略了中间步长为1 for i=k:n %步骤2.1 计算下三角矩阵L第k列元素的一部分 s(i)=A(i,k)-L(i,1:k-1)*U(1:k-1,k); end %步骤2.2 选列主元 [s_q,q]=max(abs(s(k:n))); %[M,I] = max(___) 使用先前语法中的任何输入参数，查找 A 的最大值的索引，并在输出矢量 I 中返回这些索引。 %如果最大值出现多次，则 max 返回对应于第一次出现位置的索引。 q=q+k-1;%还原s(i)的索引下标 %步骤2.3 若q!=k if q>k %绝对值最大的列主元不在第q行即当前第一行 %步骤2.3.1：交换A的第k行与第q行 t_A=A(k,:);A(k,:)=A(q,:);A(q,:)=t_A; %步骤2.3.2：交换下三角矩阵L第k列元素s_k与s_q（l_kk与l_qk） t_s=s(k);s(k)=s(q);s(q)=t_s; %步骤2.3.3：将下三角矩阵L的第k行前(k-1)个元素与第q行的前(k-1)个元素交换 t_L=L(k,1:k-1);L(k,1:k-1)=L(q,1:k-1);L(q,1:k-1)=t_L; %步骤2.3.4：交换b_k与b_q t_b=b(k);b(k)=b(q);b(q)=t_b; end %步骤2.4：计算矩阵U的第k行元素及矩阵L的第k列元素 U(k,k)=s(k); for i=k+1:n U(k,i)=A(k,i)-L(k,1:k-1)*U(1:k-1,i); L(i,k)=s(i)/U(k,k) end end %步骤3：求解上三角矩阵第n行的唯一待定元素u_nn U(n,n)=A(n,n)-L(n,1:n-1)*U(1:n-1,n); %步骤4：解单位下三角方程组Ly=b y=zeros(n,1); for k=1:n if k==1 y(1)=b(1); else y(k)=b(k)-L(k,1:k-1)*y(1:k-1); end end %步骤5：解上三角方程组Ux=y; x=zeros(n,1); for k=n:-1:1 if k==n x(n)=y(n)/U(n,n); else x(k)=(y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k); end end %步骤6：输出解x x=x';%把列向量转置成行向量，尽量不占用命令行的显示空间。

3、解线性方程组系数矩阵为三对角矩阵的追赶法

%ForwardBackward.m 追赶法解满足行对角占优条件的三对角矩阵的MATLAB程序 %行对角占优条件：(1)比较绝对值，第一行的主对角元大于第一行的另一个非零元，且绝对值均不为零； %(2)比较绝对值，最后一行的主对角元大于最后一行的另一个非零元，且绝对值均不为零； %(3)第一行和最后一行之外的其余各行，主对角元的绝对值大于另外两个非零元的绝对值之和且均不为0。 %说明：矩阵A的主对角线元素为d_1~d_n;次下对角线元素为a_2~a_n;次上对角线元素为c_1~c_n-1。 function x=ForwardBackward2_3(a,b,c,d) n=length(d); %步骤1：输入4个一维数组a,c,d(存放A的三条对角线上的元素),b(Ax=b，列向量)，由算例程序给出 %步骤2：若d_1==0，则输出信息“方法失败！”，算法终止。 if d(1)==0 warning('方法失败！') return end %步骤3：计算p1=d1,q1=c1/d1 p(1)=d(1);q(1)=c(1)/d(1); %步骤4：对i=2,3,...,n-1，执行步骤4.1~4.3 for i=2:n-1 %步骤4.1：计算p_i p(i)=d(i)-a(i)*q(i-1); %步骤4.2：若p_i=0,则输出信息“方法失败！"，算法终止(可以证明系数矩阵A满足对角占优条件，则p_1~p_n均不为0) if p(i)==0 warning('方法失败！') return end %步骤4.3：计算q_i q(i)=c(i)/p(i); end %步骤5：计算p_n=d_n-a_n*q_n-1。 p(n)=d(n)-a(n)*q(n-1); %步骤6：若p_n=0，则输出信息“方法失败！”，算法终止。 if p(n)==0 warning('方法失败！') return end %步骤7：解下三角方程组Ly=b y(1)=b(1)/p(1); for i=2:n y(i)=(b(i)-a(i)*y(i-1))/p(i); end %步骤8：解上三角方程组Ux=y x(n)=y(n); for i=n-1:-1:1 x(i)=y(i)-q(i)*x(i+1); end x=x'

4、Jacobi迭代法与 Gauss-Seidel迭代法

function x = Jacobi3_2(A,b,x0,eps,N) %功能：用Jacobi迭代法解n阶线性方程组Ax=b n=length(b);x=ones(n,1);k=0; %步骤1：输入系数矩阵A，右端向量b,以及初始向量x0，精读eps,以及最大迭代次数N；令k=1 %步骤2：当k≤N时，执行步骤2.1~2.3。 while k

【本文地址】

公司简介

联系我们