清華大學(xué)數(shù)學(xué)實驗小結(jié)：科學(xué)計算中的基本概念課件

大學(xué)數(shù)學(xué)實驗MathematicalExperiments小結(jié)：科學(xué)計算中的基本概念大學(xué)數(shù)學(xué)實驗MathematicalExperimentsAJokeAJokeAnotherJokeAnotherJoke科學(xué)計算研究的核心問題算法的構(gòu)造算法的分析

了解、理解、選擇、使用分析、改進、創(chuàng)造、創(chuàng)新 構(gòu)造算法的基本手段:近似研究算法的核心問題:近似對計算結(jié)果的影響

科學(xué)計算研究的核心問題算法的構(gòu)造構(gòu)造算法的基本手段:近似科學(xué)計算中的基本概念收斂性(or復(fù)雜度)-----誤差估計和分析-----收斂速度病態(tài)性穩(wěn)定性研究的出發(fā)點:誤差!!科學(xué)計算中的基本概念收斂性(or復(fù)雜度)研究的出發(fā)點:誤差計算地球的表面積:A=4πr2模型誤差：地球被看成是一個球地球的簡單理想模型

測量誤差:測量儀器誤差和前面的計算誤差地球的半徑要經(jīng)過測量和計算得到

截斷誤差：公式中的π是無理數(shù)，取為3.14舍入誤差：浮點數(shù)的計算誤差計算地球的表面積:A=4πr2浮點數(shù)(L≤s≤U)一般非零β進制數(shù)(d1非0)(非零)浮點數(shù)f尾數(shù)階碼IEEE標準(754):2進制(Float)符號尾數(shù)階碼位移(Bias)single1238127double152111023

±(1+f)·2e,0≤f<1-126≤e<127-1022≤e<1023(階碼全0或全1保留)浮點數(shù)(L≤s≤U)一般非零β進制數(shù)(d1非0)(非零)浮浮點數(shù)取β=2,t=4，L=-3，U=3，(正)浮點數(shù)的集合為

特點：分布不均勻但各處相對誤差一樣浮點數(shù)取β=2,t=4，L=-3，U=3，(正)浮點數(shù)的浮點數(shù)β=2,如果取t=6，L=-3，U=4，這時采用對數(shù)坐標，則集合Ｆ（正數(shù)部分）為能夠精確表達的數(shù)總是有限的！MATLABeps=2-52=2.2204e-016realmin=2-1022=2.2251e-308realmax=(2-eps)*21023=1.7977e+308浮點數(shù)β=2,如果取t=6，L=-3，U=4，這時采用對浮點數(shù)運算---example(MATLAB)Results浮點數(shù)運算---example(MATLAB)Resul浮點數(shù)運算---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)浮點數(shù)運算---example(MATLAB)Resul浮點數(shù)運算---example(MATLAB)Results浮點數(shù)計算滿足（加法，乘法）交換率？結(jié)合律？分配律？NO!!浮點數(shù)運算---example(MATLAB)Resul復(fù)雜度—行列式，anexample回憶:2階問題,3階問題考慮一般矩陣的行列式計算需要的乘法次數(shù)復(fù)雜度—行列式，anexample回憶:2階問題,3復(fù)雜度指數(shù)型算法算法計算量是問題規(guī)模的指數(shù)函數(shù)只能夠處理規(guī)模很小的問題多項式型算法算法計算量是問題規(guī)模的多項式函數(shù)可以處理規(guī)模較大的問題復(fù)雜度指數(shù)型算法復(fù)雜度—examples

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

)

PlotaScatter-plotO(n)

CalculateMeans,Variances,KernelDensityEstimatesO(nlog(n))

CalculateFastFourierTransformsO(nc)

CalculateSingularValueDecompositionofan rxcMatrix;SolveaMultipleLinearRegressionO(n2)

SolvemostClusteringAlgorithmsO(an)

DetectMultivariateOutliersAlgorithmicComplexity復(fù)雜度—examplesO(n1/2) PComplexityComplexityComplexityComplexityComplexityComplexity復(fù)雜度

----對于直接方法的度量標準Ax=b的Gauss消去法線性規(guī)劃問題的Simplex方法組合優(yōu)化的問題和方法例：多項式計算f(x)=ax^5+bx^4+cx^3+dx^2+ex+ff(x)=((((ax+b)x+c)x+d)x+e)x+f復(fù)雜度

----對于直接方法的度量標準Ax=b的Gauss收斂性

----刻劃算法的另外一個重要概念誤差收斂性收斂性

----刻劃算法的另外一個重要概念誤差拉格朗日插值多項式的不收斂性Runge現(xiàn)象Example拉格朗日插值多項式的不收斂性Runge現(xiàn)象Example梯形公式和辛普森公式的收斂性（非零常數(shù)）則稱In是p階收斂的。梯形公式2階收斂，辛普森公式4階收斂。若對I某個數(shù)值積分有梯形公式和辛普森公式的收斂性（非零常數(shù)）則稱In是p階線性方程組迭代法的收斂性一般迭代形式線性方程組迭代法的收斂性一般迭代形式微分方程初值問題算法的收斂性與收斂速度

一個算法的局部截斷誤差為O(hp+1)該算法具有p階精度

局部截斷誤差精度向前歐拉公式O(h2)1階向后歐拉公式O(h2)1階梯形公式O(h3)2階改進歐拉公式O(h3)2階經(jīng)典龍格-庫塔公式O(h5)4階微分方程初值問題一個算法的該算法具有p階精度局部截斷誤差精病態(tài)性

-----刻劃模型的概念考慮如下的問題f(x)=(x-1)(x-2)…….(x-20)顯然方程f(x)=0的解是1234………1920請問:如下方程的解是什么？病態(tài)性

-----刻劃模型的概念考慮如下的問題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較小的根不敏感較大的根較敏感病態(tài)??！=10e-5=10e-8較小的根不敏感病態(tài)??！病態(tài)的線性方程組（1,1）01.201.121=+xxx22x120201.122121=+=+xxxxx對b的擾動敏感(2,0)病態(tài)的線性方程組（1,1）01.201.121=+xxx22穩(wěn)定性

-----刻劃算法的關(guān)鍵概念考慮如下的序列可以證明穩(wěn)定性

-----刻劃算法的關(guān)鍵概念考慮如下的序列兩個算法

----有什么差別，哪個可以用??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科學(xué)結(jié)論的取得，不能依靠感覺簡單的計算發(fā)現(xiàn)，可以使用的算法是--

Algorithm2！計算中誤差并不可怕，重要的是誤差在算法中的傳播。穩(wěn)定----算法中產(chǎn)生的任何誤差，對后續(xù)計算的影響是衰減或可以控制的。不穩(wěn)定的算法=不能用的垃圾！科學(xué)結(jié)論的取得，不能依靠感覺簡單的計算發(fā)現(xiàn)，可以使用的算法是微分方程初值問題數(shù)值算法的穩(wěn)定性穩(wěn)定性計算中舍入誤差不會隨步數(shù)的增加無限增大yn的誤差

n

>0

微分方程穩(wěn)定向前歐拉公式向后歐拉公式經(jīng)典龍格-庫塔公式算法穩(wěn)定(特征根-

)h任意微分方程初值問題數(shù)值算法的穩(wěn)定性穩(wěn)定性計算中舍入誤差不會隨思考與練習(xí)二次代數(shù)方程的求根，下面公式分別在何時更好？sin(x)用下面公式計算，如何控制精度？exp(x):x<0時，用下面公式計算好不好？思考與練習(xí)二次代數(shù)方程的求根，下面公式分別在何時更好？sin附：一道習(xí)題的求解附：一道習(xí)題的求解假設(shè)小船始終向?qū)Π赌繕饲斑M:如右圖,由于河水的流動,所以小船實際走的是一條曲線小船開始坐標為A(0,d),終點坐標為B(0,0).小船行至點P(x,y)，記向量BP與x正方向的夾角為α.小船x方向的速度為dx/dt，y方向的速度為dy/dt微分方程數(shù)值解小船過河a)建模xyAB0v1v2vP(x,y)假設(shè)小船始終向?qū)Π赌繕饲斑M:小船開始坐標為A(0,d),小船過河xyAB0v1v2vP(x,y)小船過河xyAB0v1v2vP(x,y)小船軌跡小船軌跡d=100(米),v1=1,v2=2(米/秒)，在[0，T]內(nèi)求數(shù)值解問題：在MATLAB中T=?(如用ode45求解)當T大于小船由A點到達B點所需時間時，出現(xiàn)“異?！碑攜由正變負時，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