有限基础分析 4

NUMERICALCOMPUTATIONMETHODS第2章数值计算方法概论从误差分析到方程求解——构建有限元分析的数值基础工学院·有限元方法课程本章内容CONTENTS01引言数值计算的意义、误差来源与控制原则02插值与拟合拉格朗日插值、Runge现象、分段插值与最小二乘拟合03数值微分与数值积分三种差商近似、牛顿-柯特斯公式、梯形与Simpson法04线性方程组的求解高斯消去法、三角分解、雅可比与高斯-赛德尔迭代05非线性问题求解二分法、一般迭代法、牛顿迭代法及收敛性对比2.1引言为什么需要数值计算？误差从何而来？如何控制？为什么需要数值计算？多数实际物理问题难以通过理论推导获得解析解，需采用数值计算的三类典型场景01无解析解的模型数学模型本身可积，但无初等函数表达的解析解。典型示例：不定积分，如∫e^(-x²)dx无法用初等函数表示。02计算量爆炸的问题虽有解析解，但随问题规模增大，计算量激增，累计误差不可忽视。典型示例：n阶线性方程组求解，阶次增大时计算量呈指数增长。03离散数据驱动实验采集数据本身为离散点，无法直接通过解析法处理。典型示例：载荷-时间实验数据，需插值、拟合等方法处理。核心思想：将不可直接求解的数学模型转化为计算机可执行的算术与逻辑运算序列算法差异显著影响计算效率与精度效率对比：秦九韶算法多项式f(x)=a₀xⁿ+a₁xⁿ⁻¹+…+aₙ原始算法：需n(n+1)/2次乘法秦九韶算法：仅需n次乘法递推形式：f(x)=(…((a₀x+a₁)x+a₂)x+…)x+aₙ精度对比：等价多项式的计算差异P(x)=(x-2)⁹在[1.92,2.08]上计算展开为Q(x)=x⁹-18x⁸+144x⁷-…-512尽管P(x)≡Q(x)数学上完全等价，但计算机浮点计算显示两者曲线存在明显误差。a)P(x)=(x-2)⁹b)Q(x)展开式图2-1计算机绘图对比—算法差异导致的精度差异误差来源与数值算法设计原则三类误差来源模型误差与观测误差—源于问题设定与数据获取，难以通过数值方法弱化截断误差—数值方法本身的近似误差，可通过算法设计优化舍入误差—计算机存储限制造成，需选用稳定性好的方法控制五条算法设计原则①尽量减少计算步骤（如秦九韶算法）②防止大数"吃"小数（如10⁸+0.01中0.01被忽视）③避免相近数相减（有效数字急剧减少）④避免绝对值很小的数作分母（浮点溢出）⑤选用稳定性更好的数值方法构建优质数值算法的核心：从截断误差与舍入误差两个可优化环节入手，提升计算精度与可靠性2.2插值与拟合离散数据如何转化为连续函数？插值经过所有点，拟合逼近所有点插值的概念与几何意义问题：已知离散点，如何求未知点的函数值？设y=f(x)在[a,b]上连续，在n+1个互不相同的点x₀,x₁,…,xₙ上取值y₀,y₁,…,yₙ。若存在简单函数P(x)，使得P(xᵢ)=yᵢ(i=0,1,…,n)，则称P(x)为f(x)的插值函数。[a,b]为插值区间，x₀…xₙ为插值节点，f(x)为被插函数。常见插值类型：多项式插值·三角插值·分段插值本节重点：拉格朗日插值多项式（线性插值、抛物线插值）几何意义：插值函数经过所有插值点图2-2插值的几何意义拉格朗日插值—线性插值已知两点，构造一次多项式已知两点(x₀,f₀)和(x₁,f₁)，构造一次多项式P₁(x)=a+bx满足：拉格朗日形式：截断误差：R₁(x)=(f″(ξ)/2!)·(x-x₀)(x-x₁),ξ∈(x₀,x₁)数据点越密集，误差越小Matlab实现：xi=interp1(t,x,ti,'linear');yi=interp1(t,y,ti,'linear');被插函数与插值函数均通过2个插值点图2-3线性插值例2-1：线性插值结果%二维插值示例clearall;closeall;clc;%定义二维数据点x=[1.2,2.1,3.5,4.8,5.2,...6.7,7.3,8.9,9.1,10.0];y=[2.5,3.7,1.8,4.2,5.9,...6.1,7.8,8.3,9.5,10.2];t=1:length(x);%点的顺序ti=linspace(min(t),max(t),100);%线性插值xi=interp1(t,x,ti,'linear');yi=interp1(t,y,ti,'linear');%可视化plot(x,y,'ro','MarkerSize',8);holdon;plot(xi,yi,'b-','LineWidth',2);10个原始数据点的线性插值路径图2-4线性插值结果—红色圆点为原始数据，蓝色曲线为插值结果抛物线插值抛物线插值：已知三点，构造二次多项式已知三点(x₀,f₀),(x₁,f₁),(x₂,f₂)，构造P₂(x)=a+bx+cx²，满足拉格朗日形式：截断误差：R₂(x)=(f‴(ξ)/3!)·(x-x₀)(x-x₁)(x-x₂)ξ∈(x₀,x₂)，数据点越密集误差越小图2-5抛物线插值三次函数插值与分段插值例2-2：三次函数插值(cubic)图2-7三次函数插值结果Matlab实现：xi=interp1(t,x,ti,‘cubic');yi=interp1(t,y,ti,‘cubic');%二维插值示例clear

all;close

all;clc;

%定义二维数据点x=[1.2,2.1,3.5,4.8,5.2,6.7,7.3,8.9,9.1,10.0];%x坐标y=[2.5,3.7,1.8,4.2,5.9,6.1,7.8,8.3,9.5,10.2];%y坐标

%创建参数t表示点的顺序t=1:length(x);

%创建细粒度参数用于插值ti=linspace(min(t),max(t),100);%对x和y分别进行插值

%三次函数插值xi=interp1(t,x,ti,'cubic');yi=interp1(t,y,ti,'cubic');%可视化原始数据点和插值路径figure('Position',[100,100,800,600]);

%原始数据点和插值路径plot(x,y,'ro','MarkerSize',8,'MarkerFaceColor','r');holdon;plot(xi,yi,'b-','LineWidth',2);xlabel('X坐标');ylabel('Y坐标');gridon;legend('原始数据点','三次函数插值');holdoff;

分段插值图2-8分段线性插值问题：插值点增多→多项式阶数升高→计算量呈指数增长且可能出现龙格(Runge)现象：高阶多项式在区间端点附近产生剧烈振荡。图2-6Runge现象解决：分段插值Matlab实现：（三次函数分段插值）xi=interp1(t,x,ti,'pchip');yi=interp1(t,y,ti,'pchip');图2-9三次函数分段插值结果拟合：最小二乘法插值vs拟合：插值：函数严格经过所有数据点拟合：函数尽可能逼近所有点，不要求经过每个点最小二乘法原理给定数据集(xᵢ,yᵢ)，设拟合函数为φ(x)偏差平方和：当J取最小值时，φ(x)即为最优拟合函数。线性拟合：φ(x)=a+bx，通过求偏导为零解得a,b实验数据呈波动但整体遵循线性规律图2-10线性拟合%线性拟合p=polyfit(x,y,1);y1=polyval(p,x);plot(x,y,'ro');holdon;plot(x,y1,'b-');例2-4：线性拟合效果对比图2-11拟合效果对比已知某实验，获取数据点x=[1.2,2.1,3.5,4.8,5.2,6.7,7.3,8.9,9.1,10.0]，y=[2.5,3.7,1.8,4.2,5.9,6.1,7.8,8.3,9.5,10.2]%二维数据线性拟合示例clearall;closeall;clc;

%定义二维数据点（仅X和Y坐标）x=[1.2,2.1,3.5,4.8,5.2,6.7,7.3,8.9,9.1,10.0];%x坐标y=[2.5,3.7,1.8,4.2,5.9,6.1,7.8,8.3,9.5,10.2];%y坐标

%拟合p=polyfit(x,y,1);%一次多项式（线性）拟合，数字越大，阶数越高y1=polyval(p,x);plot(x,y,'ro','MarkerSize',8,'MarkerFaceColor','r');holdon;plot(x,y1,'b-','LineWidth',2)holdoff

xlabel('X坐标');ylabel('Y坐标');gridon;legend('原始数据点','线性拟合');

2.3数值微分与数值积分离散数据如何求导？曲线下面积如何近似计算？数值微分：三种差商近似导数定义：本质：数据点x0处的斜率近似。Matlab实现：向前差分：dy(i)=(y(i+1)-y(i))/h向后差分：dy(i)=(y(i)-y(i-1))/h中心差分：dy(i)=(y(i+1)-y(i-1))/(2*h)图2-12三种数值微分对比—a)向前b)向后c)中心向前差商向后差商中心差商例2-5：数值微分结果对比以正弦函数y=sin(x)为例，其解析导数为y'=cos(x)。对比三种数值微分方法与解析解的差异，中心差分精度最高。%中心差分(i=2到n-1)fori=2:length(y)-1dy_central(i)=(y(i+1)-y(i-1))/(2*h);enddy_central(1)=dy_central(2);%第一点使用后一点的值dy_central(end)=dy_central(end-1);%最后一点使用前一点的值%计算误差error_forward=abs(dy_forward-dy_exact);error_backward=abs(dy_backward-dy_exact);error_central=abs(dicentra-dy_exact);%可视化结果figure('Position',[100,100,1200,800]);

%绘制原函数subplot(2,2,1);plot(x,y,'b-','LineWidth',2);title('原函数y=sin(x)');xlabel('x');ylabel('y');

gridon;

%绘制导数subplot(2,2,2);plot(x,dy_exact,'k-','LineWidth',2,'DisplayName','解析解');holdon;plot(x,dy_forward,‘r--','LineWidth',1.5,'DisplayName','向前差分');plot(x,dy_backward,'g--','LineWidth',1.5,'DisplayName','向后差分');plot(x,dy_central,'m--',

'LineWidth',1.5,'DisplayName','中心差分');title('导数计算结果对比');xlabel('x');ylabel('dy/dx');gridon;legend;%绘制误差subplot(2,2,3);semilogy(x,error_forward,'r-','LineWidth',1.5,'DisplayName','向前差分误差');holdon;semilogy(x,error_backward,'g-','LineWidth',1.5,'DisplayName','向后差分误差');title('数值微分误差对比(对数刻度)');semilogy(x,error_central,'m-','LineWidth',1.5,'DisplayName','中心差分误差');xlabel('x');ylabel('绝对误差');gridon;legend;

%误差统计subplot(2,2,4);bar([123],[mean(error_forward)mean(error_backward)mean(error_central)]);title('平均误差比较');xticklabels({'向前差分','向后差分','中心差分'});ylabel('平均绝对误差');gridon;

%数值微分示例clearall;closeall;clc;

%定义函数和采样点x=linspace(0,2*pi,100);%采样点y=sin(x);%原函数值

%计算理论导数（解析解）dy_exact=cos(x);%sin(x)的导数是cos(x)%计算步长h=x(2)-x(1);%采样间隔%使用不同方法计算数值导数dy_forward=zeros(size(y));%向前差分dy_backward=zeros(size(y));%向后差分dy_central=zeros(size(y));%中心差分%向前差分(i=1到n-1)fori=1:length(y)-1dy_forward(i)=(y(i+1)-y(i))/h;enddy_forward(end)=dy_forward(end-1);%最后一点使用前一点的值%向后差分(i=2到n)fori=2:length(y)dy_backward(i)=(y(i)-y(i-1))/h;enddy_backward(1)=dy_backward(2);%第一点使用后一点的值例2-5：数值微分结果对比图2-13数值微分结果数值积分：牛顿-柯特斯公式积分几何意义：曲线与坐标轴围成的面积核心思想：将区间[a,b]分割为若干小区域，通过求和近似面积。牛顿-柯特斯公式：其中Cᵢ⁽ⁿ⁾为柯特斯系数，仅与分割数n有关n=1梯形公式：用直线近似曲边梯形，O(h²)精度n=2Simpson公式：用抛物线近似，O(h⁴)精度Matlab实现：矩形法：I=h*sum(y(1:end-1))梯形法：I=trapz(x,y)Simpson法：自定义函数，要求n为奇数图2-14数值积分图2-15梯形公式与Simpson公式例2-6：三种积分方法比较以f(x)=sin(x)在[0,π]上的积分为例（解析解=2），对比矩形法、梯形法、Simpson法在不同采样点数下的误差与收敛率。收敛率：矩形法O(h)/梯形法O(h²)/Simpson法O(h⁴)—Simpson法精度最高%对每个数据点循环ifmod(n,2)==1I_simp=simpson(x,y);%调用自定义Simpson函数else

%如果n为偶数，去掉最后一个点使其变为奇数I_simp=simpson(x(1:end-1),y(1:end-1));end%自定义Simpson函数（处理离散数据）functionI=simpson(x,y)n=length(x);ifmod(n,2)==0error('Simpson法要求数据点数n为奇数（偶数个区间）');end

h=(x(end)-x(1))/(n-1);I=h/3*(y(1)+y(end)+4*sum(y(2:2:end-1))+2*sum(y(3:2:end-2)));end图2-16三种积分方法的比较2.4线性方程组的求解有限元分析的核心：直接法vs迭代法，如何选择？直接求解法：高斯消去1.高斯消去法通过初等行变换将增广矩阵[A|b]化为上三角形式：步骤：①消元—逐列消去下方元素②回代—从最后一行向上求解计算量：O(n³)，适用于中小规模稠密矩阵增广矩阵回代求解矩阵形式对应的方程组变为直接求解法：三角分解法2.三角分解法(LU分解)则有则有以此类推，第n-1步时A(n)=Ln-1A(n-1)=Ln-1Ln-2A(n-2)=…=Ln-1Ln-2…L2L1A(1)A=A(1)=L1-1L2-1…Ln-1-1A(n)=LU。其中：L=L1-1L2-1…Ln-1-1,U=A(n)。三阶方程组五种方法对比例2-7：三阶方程组五种方法对比A=[4,-1,1;-1,4,-2;1,-2,4]b=[12;-1;5]线性方程组求解结果：----------------------------------------方法1（左除）：x=[3.000000,1.000000,1.000000],误差=0.0000000000方法2（求逆）：x=[3.000000,1.000000,1.000000],误差=0.0000000000方法3（LU分解）：x=[3.000000,1.000000,1.000000],误差=0.0000000000方法4（Cholesky分解）：x=[3.000000,1.000000,1.000000],误差=0.0000000000方法5（QR分解）：x=[3.000000,1.000000,1.000000],误差=0.0000000000----------------------------------------推荐使用左除运算符，效率最高%线性方程组求解示例clearall;closeall;clc;

%定义线性方程组Ax=bA=[4,-1,1;-1,4,-2;1,-2,4];%系数矩阵（对称正定）b=[12;-1;5];%常数向量

%方法1：左除运算符（推荐方法）x1=A\b;

%方法2：矩阵求逆（不推荐，计算效率低）x2=inv(A)*b;

%方法3：LU分解[L,U]=lu(A);%LU分解y=L\b;%解Ly=bx3=U\y;%解Ux=y

%方法4：Cholesky分解（适用于对称正定矩阵）R=chol(A);%Cholesky分解，A=R'*Ry=R'\b;%解R'y=bx4=R\y;%解Rx=y%方法5：QR分解[Q,R]=qr(A);%QR分解y=Q'*b;%解Qy=bx5=R\y;%解Rx=y%验证解的正确性err1=norm(A*x1-b);err2=norm(A*x2-b);err3=norm(A*x3-b);err4=norm(A*x4-b);err5=norm(A*x5-b);

%显示结果fprintf('线性方程组求解结果：\n');fprintf('----------------------------------------\n');fprintf('方法1（左除）：x=[%.6f,%.6f,%.6f],误差=%.10f\n',x1,err1);fprintf('\n');fprintf('方法2（求逆）：x=[%.6f,%.6f,%.6f],误差=%.10f\n',x2,err2);fprintf('\n');fprintf('方法3（LU分解）：x=[%.6f,%.6f,%.6f],误差=%.10f\n',x3,err3);fprintf('\n');fprintf('方法4（Cholesky分解）：x=[%.6f,%.6f,%.6f],误差=%.10f\n',x4,err4);fprintf('\n');fprintf('方法5（QR分解）：x=[%.6f,%.6f,%.6f],误差=%.10f\n',x5,err5);fprintf('----------------------------------------\n');迭代求解法：雅可比1.雅可比迭代法给定线性方程组可以表示为移项后可得即迭代关系式%线性方程组求解示例clearall;closeall;clc;

%定义线性方程组Ax=bA=[4,-1,1;-1,4,-2;1,-2,4];%系数矩阵（对称正定）b=[12;-1;5];%常数向量

%雅可比迭代法max_iter=1000;%最大迭代次数tol=1e-10;%收敛容差x6=zeros(size(b));%初始猜测值x_prev=x6;%上一次迭代值iter=0;%迭代计数器converged=false;%收敛标志

%分解矩阵A=D+L+UD=diag(diag(A));%对角矩阵L=tril(A,-1);%下三角部分（不含对角线）U=triu(A,1);%上三角部分（不含对角线）%雅可比迭代公式：x^(k+1)=D^(-1)*(b-(L+U)*x^(k))whileiter<max_iter&&~convergedx6=D\(b-(L+U)*x_prev);error=norm(x6-x_prev)/norm(x6);%相对误差iferror<tolconverged=true;endx_prev=x6;iter=iter+1;end

%验证解的正确性err6=norm(A*x6-b);

%显示结果fprintf('线性方程组求解结果：\n');fprintf('----------------------------------------\n');fprintf('雅可比迭代：x=[%.6f,%.6f,%.6f],误差=%.10f\n',x6,err6);fprintf('雅可比迭代：迭代次数=%d,收敛=%f\n',iter,converged);fprintf('----------------------------------------\n');

迭代求解法：高斯-赛德尔2.高斯-赛德尔迭代法%线性方程组求解示例clearall;closeall;clc;%定义线性方程组Ax=bA=[4,-1,1;-1,4,-2;1,-2,4];%系数矩阵（对称正定）b=[12;-1;5];%常数向量

%高斯-赛德尔迭代法max_iter=1000;%最大迭代次数tol=1e-10;%收敛容差x7=zeros(size(b));%初始猜测值x_prev=x7+inf;%上一次迭代值iter=0;%重置迭代计数器converged_gs=false;%高斯-赛德尔迭代收敛标志%高斯-赛德尔迭代公式：(D+L)x^(k+1)=b-Ux^(k)whileiter<max_iter&&~converged_gs

%逐分量更新fori=1:length(x7)x7(i)=(b(i)-A(i,[1:i-1,i+1:end])*x7([1:i-1,i+1:end]))/A(i,i);end

error=norm(x7-x_prev)/norm(x7);%相对误差iferror<tolconverged_gs=true;endx_prev=x7;iter=iter+1;enditer_gs=iter;%保存高斯-赛德尔迭代次数

%验证解的正确性err7=norm(A*x7-b);

%显示结果fprintf('线性方程组求解结果：\n');fprintf('----------------------------------------\n');fprintf('高斯-赛德尔迭代：x=[%.6f,%.6f,%.6f],误差=%.10f\n',x7);fprintf('高斯-赛德尔迭代：迭代次数=%d,收敛=%f\n',iter_gs,converged_gs);fprintf('----------------------------------------\n');

1、雅可比迭代法求解结果：----------------------------------------雅可比迭代：x=[3.000000,1.000000,1.000000],误差=0.0000000007雅可比迭代：迭代次数=61,收敛=1.000000----------------------------------------2、高斯-赛德尔迭代法求解结果：----------------------------------------高斯-赛德尔迭代：x=[3.000000,1.000000,1.000000],误差=高斯-赛德尔迭代：迭代次数=15,收敛=1.000000----------------------------------------2.5非线性问题求解f(x)=0无统一解析解法，二分法、迭代法、牛顿法如何选择？有根区间的确定介值定理：确定有根区间设f(x)在[a,b]内连续、严格单调，且f(a)·f(b)<0→在[a,b]内f(x)=0有且仅有一个实根作图法：对于多项式函数，可直接寻找其与横坐标的交点；对于超越方程，则可先将其转化为两个独立的函数，此时两函数曲线的交点便是原方程的解。图2-17有根区间的确定如何求解？二分法二分法原理①取区间中点c=(a+b)/2②若f(c)·f(a)<0→根在[a,c]，令b=c否则→根在[c,b]，令a=c③重复直到区间长度满足精度要求终止条件：|bₖ-aₖ|<ε或|f(c)|<ε优缺点：简单易编程，收敛稳定；但收敛速度慢(线性)，无法求重根、复根图2-18二分法图解例2-10：二分法求解f(x)=x³-2x-5=0二分法求解结果初始区间：[2,3]，容差：1e-6根的近似值：x=2.0945529938函数值：f(x)=1.688e-05迭代次数：19次区间长度：1.907e-06二分法收敛特性收敛阶：线性收敛（每次区间长度减半）误差：|xₖ-x*|≤(b-a)/2ᵏ优点：全局收敛，不依赖初值，稳定可靠缺点：收敛慢，仅求实根，无法求重根适用：快速确定根的大致范围，为其他方法提供初值Matlab实现：whileiter<max_iter&&~convergedc=(a+b)/2;%计算中点fc=f(c);%计算中点的函数值%存储迭代信息iterations(iter+1,:)=[iter+1,a,b,c];

%输出当前迭代信息fprintf('%d\t\t%.8f\t%.8f\t%.8f\t%.8e\n',iter+1,a,b,c,fc);%检查收敛条件ifabs(fc)<tol||abs(c-c_prev)/abs(c)<tolconverged=true;end%更新区间iffc*f(a)<0b=c;%根在左半区间elsea=c;%根在右半区间endc_prev=c;iter=iter+1;end一般迭代法不动点迭代原理将f(x)=0等价变形为x=g(x)建立迭代格式：xₖ₊₁=g(xₖ)选取初值x₀，递推计算直至收敛几何意义：y=g(x)