科学计算与仿真应用（基于MWORKS）课件 第五章 数值计算方法

第5章数值计算方法本章目标理解并掌握在Syslab中Julia语言下的线性方程组的直接法与迭代法，以及矩阵特征值和特征向量的常用数值方法（如幂法、反幂法、Jacobi方法和QR方法）。掌握syslab中多项式插值与最小二乘曲线拟合的基本原理，能够应用这些方法进行数据拟合与函数逼近。理解数值积分与数值微分的基本方法，掌握常微分方程的Euler方法和Runge-Kutta方法，并了解相关函数的使用。同时，掌握使用julia语言对非线性方程求根的常用数值方法，如二分法、Newton迭代法和弦截法。知识点总览线性方程组的数值解法：掌握直接法（如LU分解、Cholesky分解、矩阵除法）和迭代法（Jacobi迭代、Gauss-Seidel迭代、SOR迭代）的基本原理与实现，能够针对不同类型（稠密/稀疏、正定/非正定）的线性方程组选择合适的数值解法。多项式插值与数据拟合技术：理解Lagrange插值、Hermite插值、样条插值等方法，掌握最小二乘曲线拟合的原理与实现，能够对离散数据进行函数逼近与预测。数值积分与数值微分的计算方法：掌握自适应Simpson积分、Romberg算法等数值积分方法，以及基于插值多项式的数值微分公式（如三点公式、差分法），能够对复杂函数或离散数据进行积分与导数的近似计算。矩阵特征值与特征向量的数值求解：理解幂法、反幂法、Jacobi方法、QR方法等经典算法的基本思想与适用场景，能够编写程序求解矩阵的主特征值、全部特征值及对应的特征向量。常微分方程与非线性方程的数值解法：掌握常微分方程初值问题的Euler方法、改进Euler方法、Runge-Kutta方法，以及非线性方程求根的二分法、Newton迭代法、弦截法，能够分析算法的精度、收敛性与稳定性。1.线性方程组的数值解法-直接法#将解向量中的元素转换为分数形式.x_fraction=map(x->Rational(x),x)直接法是经过有限步的算术运算求线性方程组精确解的方法（假设计算过程中没有舍入误差）。但由于在实际计算过程中存在舍入误差，因此这种方法只能求得线性方程组的近似解。示例：求解线性方程组运算规则：对于一个n阶线性方程组Ax=b，在Syslab中，只需要一个反斜杠运算符“/”或斜杠运算符“\”就可求解。#对矩阵A和b赋值.A=[5-41;-46-4;1-46];b=[2;-1;-1]x=A\b1.线性方程组的数值解法-迭代法#使用反斜杠运算符求解x.result=Jacobimethod(A,b,x0,Maxit);示例：使用Jacobi迭代法求解如下线性方程组，要求Jacobi迭代法：迭代的收敛性只和系数矩阵有关，和初值的选择没有关系，初值可任意选择。时停止：注意：需要先编辑Jacobi函数，再代入数值1.线性方程组的数值解法-迭代法#SOR更新公式.x1[i]=(1-Omega)*x0[i]+(Omega/A[i,i])*(b[i]-_sum1-_sum2);示例：求解如下线性方程组，要求SOR（SuccessiveOver-Relaxation，逐次超松弛）迭代法：是Gauss-Seidel迭代的一种加速方法，是解大型稀疏矩阵线性方程组的有效方法之一，其计算公式简单，但需要较好的加速因子（最佳松弛因子）。时停止：注意：SOR迭代法中最佳松弛因子ω的选择很困难，数值分析已经证明SOR迭代收敛的必要条件是ω在0和2之间，并且在系数矩阵是三角矩阵的情况下，给出了一个ω的计算公式。2.多项式插值与最小二乘曲线拟合-多项式插值Lagrange插值：已知个插值点及其对应的函数值若采用Lagrange插值方法，则对插值区间内任意x可以求得示例：已知数据如表5.1所示。试使用Lagrange插值方法求时，函数的近似值。

表5.1例5.6数据0.561600.562800.564010.565210.827410.826590.825770.814952.多项式插值与最小二乘曲线拟合-多项式插值functionlagrange(x,y,x0)n=length(x)m=length(x0)y0=zeros(m)k=0;i=0;j=0;forkin1:mz=x0[k];s=0.0;foriin1:np=1.0;forjin1:nifj!=ip*=(z-x[j])/(x[i]-x[j]);endends+=p*y[i];endy0[k]=s;endreturny0end使用Julia语言编写lagrange函数，之后赋值进行计算2.多项式插值与最小二乘曲线拟合-多项式插值代入数据后得到运行结果：pchip函数和spline函数的插值结果如图5.2所示。图5.22.多项式插值与最小二乘曲线拟合-最小二乘曲线拟合最小二乘曲线拟合：具体做法是，对于给定的一组数据functionlinearfit(x,y,m)nx=length(x)ny=length(y)ifnx!=nyerror("Theinputdataiswrong!");endA=zeros(m+1,m+1);b=zeros(m+1);foriin1:m+1forjin1:m+1A[i,j]=sum(x.^(i+j-2));endb[i]=sum(x.^(i-1).*y);enda=A\b;yy=polyval(reverse(a),x)，要求在函数类中找到一个函数，使误差的平方和最小，即温度/℃100110120130140150160170180190产品率/％45515461667074788589示例：已知某种化学反应过程中温度与产品率关系的数据如表5.2所示。请用最小二乘曲线拟合得出函数关系。2.多项式插值与最小二乘曲线拟合-最小二乘曲线拟合小二乘曲线拟合结果如图5.3所示图5.33.数值积分与数值微分-数值积分（1）quad函数。quad函数是Syslab自带的函数，其功能是使用自适应Simpson公式求给定函数的定积分。#q=quad(fun,a,b);q=quad(fun,a,b;Name=Value);语法规则：根据节点的选择方法，数值分析中的数值积分方法有Newton-Cotes公式法、Romberg算法、Gauss公式法等。示例：求积分注意：在运行脚本前需要安装QuadGK函数库，其语法格式如下：#julia>importPkg;julia>Pkg.add("QuadGK")usingTyMathusingQuadGKf(x)=sqrt(4-sin(x)^2)result,error=quadgk(f,0,pi/6)println("积分结果:",result)println("估计误差:",error)3.数值积分与数值微分-数值积分（2）Romberg算法和函数。usingTyMathfunctionromberg_integration(f,a,b,n=6)R=zeros(n,n);h=(b-a)./(2.^(0:n-1));R[1,1]=(b-a)*(f(a)+f(b))/2;foriin2:nsum_mid=0.0;step=2^(n-i);forjin1:2^(i-2)x=a+(2*j-1)*step*h[i];sum_mid+=f(x);endR[i,1]=R[i-1,1]/2+sum_mid*h[i];forkin2:iR[i,k]=R[i,k-1]+(R[i,k-1]-R[i-1,k-1])/(4^(k-1)-1);endendintegral=R[n,n];error=abs(R[n,n]-R[n,n-1]);returnintegral,errorendclear()clc()f(x)=sqrt(4-sin(x)^2)result,error=romberg_integration(f,0,pi/6,6)println("Lundberg积分结果:",result)println("估计误差:",error)3.数值积分与数值微分-数值微分本节讨论数值微分的方法。插值型求导公式用插值多项式的导数近似表示，即#functionnpointdiff(x,y,x0)A=copy(y);N=length(x);forjin2:NforkinN:-1:j#Calculatethej-thorderdifferenceatpointk.A[k]=(A[k]-A[k-1])/(x[k]-x[k-j+1]);endendX=x[1];df=A[2];prod=1;n1=length(A)–1;forkin2:n1prod*=(x0-x[k])df+=prod*A[k+1];endreturnA,dfend基于

个点的差分求解程序如下。示例：已知函数的数据如表5.4所示，计算处的一阶导数。x2.52.62.72.82.9y12.18513.463714.879716.444618.1741usingTyMathusingTyBasex=[2.5,2.6,2.7,2.8,2.9];y=[12.185,13.4637,14.8797,16.4446,18.1741];A,df=npointdiff(x,y,2.7);println("A是:",A)returnAprintln("df是:",df)returndf4.矩阵的特征值和特征向量

Syslab中求解特征值和特征向量的函数如下：eig函数：求矩阵的特征值和特征向量。eigvals函数：求矩阵的特征值。eigvecs函数：求矩阵的特征向量。svd函数：计算矩阵的奇异值分解。schur函数：Schur分解。示例：求矩阵的特征值和特征向量。clear()clc()usingTyMathusingTyBaseA=[2-10;-12-1;0-12];V,D=eig(A);println("矩阵A的特征向量为：",V)returnVprintln("矩阵A的特征值为：",D)returnD4.矩阵的特征值和特征向量-幂法和反幂法functionpowermethod(A,x0,Eps=1e-6,maxtimes=500)Lambda=0;cnt=0;err=1;state=1;while(cnt<=maxtimes)&&(state==1)Y=A*x0;m,j=findmax(abs.(Y));c1=mdc=abs(Lambda-c1);Y=(1/c1)*Y;dv=norm(x0-Y);err=max(dc,dv);x0=Y;Lambda=c1;state=0;幂法是用来求矩阵的模最大的特征值及对应的特征向量的迭代法。迭代公式如下：使用幂法计算矩阵的模最大的特征值及对应的特征向量的自编函数如下：iferr>Epsstate=1endcnt+=1end#返回结果.returnLambda,x0end示例：使用幂法求矩阵的模最大的特征值及对应的特征向量。clear()clc()usingTyMathusingTyBaseA=[011-5;-217-7;-426-10];x0=[1;1;1]；Lambda,V=powermethod(A,x0,1e-6,500)println("矩阵A的模最大的特征值为：",Lambda)returnLambdaprintln("矩阵A的模最大的特征值对应的特征向量为：",V)returnV4.矩阵的特征值和特征向量-幂法和反幂法反幂法是用来求矩阵的模最小的特征值及对应的特征向量的迭代法。用替代执行上述的幂法就是所谓的反幂法。使用反幂法计算矩阵的模最小的特征值及对应的特征向量的自编函数如下：functioninvpowermethod(A,x0,alpha,Eps,maxtimes)

n,_=size(A);A=A-alpha*I;Lambda=0;cnt=0;err=1;state=1;while(cnt<=maxtimes)&&(state==1)Y=A\x0;m,j=findmax(abs.(Y));c1=m;dc=abs(Lambda-c1);Y=(1/c1)*Y;dv=norm(x0-Y);err=max(dc,dv);x0=Y;

Lambda=c1;state=0;iferr>Epsstate=1;endcnt+=1;endLambda=alpha+1/c1;V=x0;returnLambda,Vend

4.矩阵的特征值和特征向量-Jacobi方法Jacobi方法是用来计算实对称矩阵的全部特征值和特征向量的一种迭代法。使用Jacobi方法计算矩阵的特征值和特征向量的自编函数如下：functionjacobimethod(A,tol=1e-8,maxiter=1000)n=size(A,1);V=Matrix{Float64}(I,n,n);iters=0;whiletruemaxval=0.0;p,q=0,0;foriin1:n-1forjini+1:nifabs(A[i,j])>maxvalmaxval=abs(A[i,j]);p,q=i,j;endendendifmaxval<tol||iters>=maxiterbreakendθ=0.5*atan(2*A[p,q],A[q,q]-A[p,p]);R=Matrix{Float64}(I,n,n);c,s=cos(θ),sin(θ);R[p,p]=c;R[q,q]=c;R[p,q]=s;R[q,p]=-s;A=R'*A*R;V*=R;iters+=1;endidx=sortperm(diag(A));eigenvalues=diag(A)[idx];eigenvectors=V[:,idx];println("矩阵的特征值为:",eigenvalues)println("矩阵的特征向量为:",eigenvectors)end4.矩阵的特征值和特征向量-QR方法QR方法是目前求一般方阵的全部特征值最有效且被广泛应用的方法之一。通过一系列正交相似变换得到方阵序列usingLinearAlgebraXfunctionqrmethod(A::Matrix{<:Number};max_iter=1000,tol=1e-10)Ak=copy(A)n=size(Ak,1);foriterin1:max_iterQ,R=qr(Ak);Ak=R*Q;sub_diag_norm=norm([Ak[i+1,i]foriin1:n-1]);ifsub_diag_norm<tolbreakendendreturndiag(Ak)end{Ak}，在一定条件下，可以证明当，Ak时本质收敛于一个块对角矩阵，其对角子块是1阶或2阶方阵，因此容易求出它们的特征值。使用QR方法计算矩阵的特征值的自编函数如下。5.常微分方程的数值解法-Euler方法对于简单的一阶方程的初值问题，即functioneulermethod(f,x0,y0,h,n_steps)x_values=zeros(n_steps+1);y_values=zeros(n_steps+1);x_values[1]=x0;y_values[1]=y0;foriin1:n_stepsx=x_values[i];y=y_values[i];y_next=y+h*f(x,y);x_next=x+h;x_values[i+1]=x_next;y_values[i+1]=y_next;endreturnx_values,y_valuesend，由Euler方法可得使用Euler方法求解一阶方程的初值问题的自编函数如下。5.常微分方程的数值解法-Euler方法示例：使用Euler方法求解：并与精确解进行比较。根据自编程序代入数据后得到运行结果如图5.4所示：图5.45.常微分方程的数值解法-Runge-Kutta方法为了提高常微分方程的数值解法的精度，可间接使用Taylor展开的思想构造高阶的数值解法，四阶Runge-Kutta方法的经典公式如下：5.常微分方程的数值解法-Runge-Kutta方法使用Runge-Kutta方法求解常微分方程的自编函数如下。functionrungekuttamethod(f,a,b,y0,h)x=a:h:b;n=length(x);y=zeros(n);y[1]=y0foriin1:n-1k1=f(x[i],y[i]);k2=f(x[i]+h/2,y[i]+h/2*k1);k3=f(x[i]+h/2,y[i]+h/2*k2);k4=f(x[i+1],y[i]+h*k3);y[i+1]=y[i]+(k1+2*k2+2*k3+k4)*h/6;endreturnx,yendfunctionmyfun(x,y)return-1.0/xend5.常微分方程的数值解法表5.5例5.17的三种数值解法的对比i精确值Euler方法的误差改进的Euler方法的误差Runge-Kutta方法的误差01-100011.1-0.90910.0091-0.0009-0.000002405721.2-0.83330.0134-0.0013-0.000003416631.3-0.76920.0152-0.0015-0.000003723241.4-0.71430.0156-0.0015-0.000003676951.5-0.66670.0153-0.0015-0.00000346095.常微分方程的数值解法-Syslab中的相关函数函数算法适用条件精度特点ode454阶、5阶Runge-Kutta方法非刚性方程中最常用的单步算法。ode113可变阶算法非刚性方程低—高多步算法，在相同的精度下，它可能比ode45函数或ode23函数运算速度快一些，但不适用于不连续的系统。ode15s可变阶算法刚性方程低—中多步算法，方程刚性，用于ode45函数运算速度很慢时。ode232阶、3阶Runge-Kutta方法非刚性方程低有时对于轻度的刚性方程来说，它比ode45函数更有效。在相同的精度下，它需要比ode45函数更小的步长。ode23s基于改进的Rosenbrock公式刚性方程低单步算法，用于不适合使用ode15s函数的刚性方程。ode23t自由内插实现的梯形规则一般刚性低给出的解无数值衰减。ode23tbTR-BDF2算法刚性方程低当允许误差范围较大时，它比ode15s函数好。表5.6常微分方程的几种数值求解函数的对比6.非线性方程求根的数值计算方法-二分法

重复上述步骤，如果通过有限步求得根，则计算中止，否则每次有根区间长度减半，并得到一个闭区间：使用二分法求解非线性方程的自编函数如下。functionbisectionmethod(f,a,b;tol=1e-6,max_iter=100)iff(a)==0returna,0,trueendiff(b)==0returnb,0,trueendiff(a)*f(b)>0error("区间[a,b]两端点函数值必须异号")endconverged=false;iter=0;c=0.0;foriterin1:max_iterc=(a+b)/2;fc=f(c);ifabs(fc)<tol||(b-a)/2<tolconverged=true;break;endiffc*f(a)<0b=c;elsea=c;endendreturnc,iter,convergedend6.非线性方程求根的数值计算方法-Newton迭代法解非线性方程的Newton迭代法是一种将非线性函数线性化的方法。若令，则线性方程的根就是非线性方程的近似根。构造Newton迭代公式：当Newton迭代法收敛时，迭代序列的极限就是非线性方程的根，即6.非线性方程求根的数值计算方法-Newton迭代法functionnewtonmethod(f,df,x0;tol=1e-6,max_iter=100,verbose=false)x=x0;converged=false;foriterin1:max_iterfx=f(x);dfx=df(x);ifabs(dfx)<1e-12println("警告:在x=$(x)处导数接近零")breakendx_next=x-fx/dfx;ifabs(x_next-x)<tol||abs(fx)<tolconverged=true;x=x_next;breakendifverboseprintln("迭代$(iter):x=$(x_next),f(x)=$(f(x_next))")endx=x_next;endreturnx,max_iter,convergedend