清华大学数学实验小结：科学计算中的基本概念课件

大学数学实验MathematicalExperiments小结：科学计算中的基本概念大学数学实验MathematicalExperimentsAJokeAJokeAnotherJokeAnotherJoke科学计算研究的核心问题算法的构造算法的分析

了解、理解、选择、使用分析、改进、创造、创新 构造算法的基本手段:近似研究算法的核心问题:近似对计算结果的影响

科学计算研究的核心问题算法的构造构造算法的基本手段:近似科学计算中的基本概念收敛性(or复杂度)-----误差估计和分析-----收敛速度病态性稳定性研究的出发点:误差!!科学计算中的基本概念收敛性(or复杂度)研究的出发点:误差计算地球的表面积:A=4πr2模型误差：地球被看成是一个球地球的简单理想模型

测量误差:测量仪器误差和前面的计算误差地球的半径要经过测量和计算得到

截断误差：公式中的π是无理数，取为3.14舍入误差：浮点数的计算误差计算地球的表面积:A=4πr2浮点数(L≤s≤U)一般非零β进制数(d1非0)(非零)浮点数f尾数阶码IEEE标准(754):2进制(Float)符号尾数阶码位移(Bias)single1238127double152111023

±(1+f)·2e,0≤f<1-126≤e<127-1022≤e<1023(阶码全0或全1保留)浮点数(L≤s≤U)一般非零β进制数(d1非0)(非零)浮浮点数取β=2,t=4，L=-3，U=3，(正)浮点数的集合为

特点：分布不均匀但各处相对误差一样浮点数取β=2,t=4，L=-3，U=3，(正)浮点数的浮点数β=2,如果取t=6，L=-3，U=4，这时采用对数坐标，则集合Ｆ（正数部分）为能够精确表达的数总是有限的！MATLABeps=2-52=2.2204e-016realmin=2-1022=2.2251e-308realmax=(2-eps)*21023=1.7977e+308浮点数β=2,如果取t=6，L=-3，U=4，这时采用对浮点数运算---example(MATLAB)Results浮点数运算---example(MATLAB)Resul浮点数运算---example(MATLAB)Resultsx=0.988:.0001:1.012;y=x.^7-7*x.^6+21*x.^5-35*x.^4+35*x.^3-21*x.^2+7*x-1;plot(x,y)浮点数运算---example(MATLAB)Resul浮点数运算---example(MATLAB)Results浮点数计算满足（加法，乘法）交换率？结合律？分配律？NO!!浮点数运算---example(MATLAB)Resul复杂度—行列式，anexample回忆:2阶问题,3阶问题考虑一般矩阵的行列式计算需要的乘法次数复杂度—行列式，anexample回忆:2阶问题,3复杂度指数型算法算法计算量是问题规模的指数函数只能够处理规模很小的问题多项式型算法算法计算量是问题规模的多项式函数可以处理规模较大的问题复杂度指数型算法复杂度—examples

DescriptorDataSetSizeinBytes StorageModeTiny 102 PieceofPaperSmall 104 AFewPiecesofPaperMedium 106 AFloppyDiskLarge 108 HardDiskHuge 1010 MultipleHardDisksMassive 1012 RoboticMagneticTape StorageSilosSuper-massive 1015 DistributedDataArchivesTheHuber-WegmanTaxonomyofDataSetSizes复杂度—examplesDescriptorDa复杂度—examplesO(n1/2

)

PlotaScatter-plotO(n)

CalculateMeans,Variances,KernelDensityEstimatesO(nlog(n))

CalculateFastFourierTransformsO(nc)

CalculateSingularValueDecompositionofan rxcMatrix;SolveaMultipleLinearRegressionO(n2)

SolvemostClusteringAlgorithmsO(an)

DetectMultivariateOutliersAlgorithmicComplexity复杂度—examplesO(n1/2) PComplexityComplexityComplexityComplexityComplexityComplexity复杂度

----对于直接方法的度量标准Ax=b的Gauss消去法线性规划问题的Simplex方法组合优化的问题和方法例：多项式计算f(x)=ax^5+bx^4+cx^3+dx^2+ex+ff(x)=((((ax+b)x+c)x+d)x+e)x+f复杂度

----对于直接方法的度量标准Ax=b的Gauss收敛性

----刻划算法的另外一个重要概念误差收敛性收敛性

----刻划算法的另外一个重要概念误差拉格朗日插值多项式的不收敛性Runge现象Example拉格朗日插值多项式的不收敛性Runge现象Example梯形公式和辛普森公式的收敛性（非零常数）则称In是p阶收敛的。梯形公式2阶收敛，辛普森公式4阶收敛。若对I某个数值积分有梯形公式和辛普森公式的收敛性（非零常数）则称In是p阶线性方程组迭代法的收敛性一般迭代形式线性方程组迭代法的收敛性一般迭代形式微分方程初值问题算法的收敛性与收敛速度

一个算法的局部截断误差为O(hp+1)该算法具有p阶精度

局部截断误差精度向前欧拉公式O(h2)1阶向后欧拉公式O(h2)1阶梯形公式O(h3)2阶改进欧拉公式O(h3)2阶经典龙格-库塔公式O(h5)4阶微分方程初值问题一个算法的该算法具有p阶精度局部截断误差精病态性

-----刻划模型的概念考虑如下的问题f(x)=(x-1)(x-2)…….(x-20)显然方程f(x)=0的解是1234………1920请问:如下方程的解是什么？病态性

-----刻划模型的概念考虑如下的问题p=poly(1:20);ep=zeros(1,21);ep(3)=1.0e-5;re=roots(p+ep)plot(re,'b+');holdonplot(1:20,0,'r*');holdoffMatlabprogramp=poly(1:20);Matlabprogram=10e-5=10e-8较小的根不敏感较大的根较敏感病态！！=10e-5=10e-8较小的根不敏感病态！！病态的线性方程组（1,1）01.201.121=+xxx22x120201.122121=+=+xxxxx对b的扰动敏感(2,0)病态的线性方程组（1,1）01.201.121=+xxx22稳定性

-----刻划算法的关键概念考虑如下的序列可以证明稳定性

-----刻划算法的关键概念考虑如下的序列两个算法

----有什么差别，哪个可以用??Algorithm1Algorithm2两个算法

----有什么差别，哪个可以用??Algorithclearep(1)=1forn=2:100ep(n)=exp(1.0)-n*ep(n-1)endplot(ep,'b*');ProgramofAlgorithm1clearProgramofAlgorithm1Algorithm1withn=15Algorithm1withn=15Algorithm1withn=100Algorithm1withn=100clearep(100)=0forn=100:-1:2ep(n-1)=(exp(1.0)-ep(n))/n;endplot(ep,'b*');ProgramofAlgorithm2clearProgramofAlgorithm2Algorithm2withn=100Algorithm2withn=100科学结论的取得，不能依靠感觉简单的计算发现，可以使用的算法是--

Algorithm2！计算中误差并不可怕，重要的是误差在算法中的传播。稳定----算法中产生的任何误差，对后续计算的影响是衰减或可以控制的。不稳定的算法=不能用的垃圾！科学结论的取得，不能依靠感觉简单的计算发现，可以使用的算法是微分方程初值问题数值算法的稳定性稳定性计算中舍入误差不会随步数的增加无限增大yn的误差

n

>0

微分方程稳定向前欧拉公式向后欧拉公式经典龙格-库塔公式算法稳定(特征根-

)h任意微分方程初值问题数值算法的稳定性稳定性计算中舍入误差不会随思考与练习二次代数方程的求根，下面公式分别在何时更好？sin(x)用下面公式计算，如何控制精度？exp(x):x<0时，用下面公式计算好不好？思考与练习二次代数方程的求根，下面公式分别在何时更好？sin附：一道习题的求解附：一道习题的求解假设小船始终向对岸目标前进:如右图,由于河水的流动,所以小船实际走的是一条曲线小船开始坐标为A(0,d),终点坐标为B(0,0).小船行至点P(x,y)，记向量BP与x正方向的夹角为α.小船x方向的速度为dx/dt，y方向的速度为dy/dt微分方程数值解小船过河a)建模xyAB0v1v2vP(x,y)假设小船始终向对岸目标前进:小船开始坐标为A(0,d),小船过河xyAB0v1v2vP(x,y)小船过河xyAB0v1v2vP(x,y)小船轨迹小船轨迹d=100(米),v1=1,v2=2(米/秒)，在[0，T]内求数值解问题：在MATLAB中T=?(如用ode45求解)当T大于小船由A点到达B点所需时间时，出现“异常”当y由正变负时，dy/dt由负变正，引起的“震荡”.b)求解d=100(米),v1=1,v2=2(米/秒)，在[解决方法1)进入B(0,0)的小邻域后立即终止程序functionxdot=chuan(t,x)v1=1;v2=2;xdot=[v1-v2*x(1)/sqrt(x(1)^2+x(2)^2);-v2*x(2)/sqrt(x(1)^2+x(2)^2)];解决方法1)进入B(0,0)的小邻域后立即终止程序functd=100;ts=0:0.1:50;x0=[0,d];[t,x]=ode45('chuan',ts,x0);n=length(x);T=50;t=ts';whilemin(abs(x(n,1)),abs(x(n,2)))>=0.0005T=T+0.1;t1=[T-0.1,T];x1=x(n,:);[t2,x2]=ode45('chuan',t1,x1);n1=length(t2);n=n+1;x=[x;x2(n1,:)];t=[t;T];endn1=length(t);[t(1:10),x(1:10,:)],[t((n1-10):n1),x((n1-10):n1,:)],Tpauseplot(t,x),grid,xlabel('t'),ylabel('x1,x2'),gtext('x1(t)'),gtext('x2(t)'),pauseplot(x(:,1),x(:,2)),grid,xlabel('x1'),ylabel('x2')pausev1=1;v2=2;k=v1/v2;x1=d/2*((x(:,2)/d).^(1-k)-(x(:,2)/d).^(1+k));[x(1:10,:),x1(1:10)],[x((n1-10):n1,:),x1((n1-10):n1)],d=100;ts=0:0.1:50;x0=[0,d];[t,xiaochuan.mtxy00100.00000.10000.099999.80000.20000.199699.60000.30000.299199.40000.40000.398499.20000.50000.497599.00000.60000.596498.80000.70000.695198.60000.80000.793698.40000.90000.891998.2000……………..txy……………..65.60001.06600.045565.70000.96610.037465.80000.86630.030065.90000.76640.023566.00000.66650.017866.10000.56650.012866.20000.46660.008766.3