工程數(shù)值分析報(bào)告實(shí)驗(yàn)龍格庫塔,最小二乘法

標(biāo)準(zhǔn)文案工程數(shù)值分析實(shí)驗(yàn)報(bào)告指導(dǎo)老師班級(jí)學(xué)號(hào)姓名大全標(biāo)準(zhǔn)文案實(shí)驗(yàn)一：最小二乘法擬合曲線實(shí)驗(yàn)一、實(shí)驗(yàn)名稱：最小二乘法擬合曲線實(shí)驗(yàn)實(shí)驗(yàn)時(shí)間：2015-5-14實(shí)驗(yàn)地點(diǎn)：主樓機(jī)房實(shí)驗(yàn)器材：計(jì)算機(jī)matlab二、實(shí)驗(yàn)?zāi)康模簩W(xué)會(huì)用最小二乘法求擬合數(shù)據(jù)的多項(xiàng)式，并應(yīng)用算法于實(shí)際問題。三、實(shí)驗(yàn)要求：（1）根據(jù)最小二乘法和加權(quán)最小二乘法的基本理論，編寫程序構(gòu)造擬合曲線的法方程，要求可以方便的調(diào)整擬合多項(xiàng)式的次數(shù)；（2）采用列主元法解（1）中構(gòu)造的法方程，給出所擬合的多項(xiàng)式表達(dá)式；（3）編寫程序計(jì)算所擬合多項(xiàng)式的均方誤差，并作出離散函數(shù) 和擬合函數(shù)的圖形；（4）用MATLAB的內(nèi)部函數(shù)polyfit求解上面最小二乘法曲線擬合多項(xiàng)式的系數(shù)及平方誤差，并用MATLAB的內(nèi)部函數(shù)plot作出其圖形，并與（1）的結(jié)果進(jìn)行比較。四、算法描述（實(shí)驗(yàn)原理與基礎(chǔ)理論）基本原理：從整體上考慮近似函數(shù)同所給數(shù)據(jù)點(diǎn)(xi,yi)(i=0,1,?,m)誤差rip(xi)yi(i=0,1,?,m)的大小，常用的方法有以下三種：一是誤差rip(xi)yi(i=0,1,?,m)絕maxri，即誤差向量r(r0,r1,rm)T對(duì)值的最大值的∞—范數(shù)；二是誤差絕對(duì)值的和0immrim2ri的算術(shù)平方根，即誤差向量r的2i0，即誤差向量r的1—范數(shù)；三是誤差平方和i0—范數(shù)；前兩種方法簡單、自然，但不便于微分運(yùn)算，后一種方法相當(dāng)于考慮2—范數(shù)的平m2方，因此在曲線擬合中常采用誤差平方和rii0來度量誤差ri(i=0，1，?，m)的整體大小。五、實(shí)驗(yàn)內(nèi)容：共有兩組給定數(shù)據(jù)，把給定的數(shù)據(jù)擬合成多項(xiàng)式。第一組給定數(shù)據(jù)點(diǎn)如表 1所示如下：xiyi

表1數(shù)據(jù)表0 0．5 0．6 0．7 0．8 0．9 1．01．751．962．192．442．713．00表2數(shù)據(jù)表大全標(biāo)準(zhǔn)文案xiyii六、程序流程圖

00．50．60．70．80．91．011．751．962．192．442．713．001236421開始輸入擬合次數(shù)N,X,Y的坐標(biāo)矩陣A 計(jì)算X,Y平均值獲得A的逆矩陣 B獲得矩陣 C計(jì)算相關(guān)系數(shù) r解方程組獲得系數(shù) a0,a1輸出結(jié)束大全標(biāo)準(zhǔn)文案七、實(shí)驗(yàn)結(jié)果zuixiaoerchenfaans=27-May-2015ans=7.3611e+05ans=1.0e+03*2.0150 0.0050 0.0270 0.0140 0.0010 0.0213>>最小二乘法3最小二乘法3y0y0yy2.52.522yy1.51.5110.50.10.20.30.40.50.60.70.80.910.50.10.20.30.40.50.60.70.80.9100xxFigure1Figure2大全標(biāo)準(zhǔn)文案最小二乘法3.5y0y32.5y 21.510.50 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1xFigure3八、實(shí)驗(yàn)結(jié)果分析實(shí)驗(yàn)程序quxiannihe.mclear alldate,now,clockx0=[0.00.50.60.70.80.91.0];y0=[11.751.962.192.442.713.00];w=ones(size(x0));x=0:0.01:1;%進(jìn)行五次曲線擬合N=5;fori=1:Na1=LSF(x0,y0,w,i);y=polyval(a1,x);figure(i)plot(x0,y0, 'ok',x,y, 'r' )

最小二乘法3y0y2.52y1.510.50 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1xFigure4title( '最小二乘法');大全標(biāo)準(zhǔn)文案legend( 'y0','y');xlabel( 'x');ylabel( 'y');end實(shí)驗(yàn)二：4階經(jīng)典龍格庫塔法解常微分方程一、實(shí)驗(yàn)名稱：4階經(jīng)典龍格庫塔法解常微分方程實(shí)驗(yàn)時(shí)間：2015-5-14實(shí)驗(yàn)地點(diǎn)： 主樓機(jī)房實(shí)驗(yàn)器材： 計(jì)算機(jī)matlab二、實(shí)驗(yàn)?zāi)康模簩W(xué)習(xí)掌握4階經(jīng)典R-K方法，體會(huì)參數(shù)和步長對(duì)問題的影響。三、實(shí)驗(yàn)要求：1）用4階經(jīng)典R-K法編寫計(jì)算程序，要求用法與ode45一致。并將計(jì)算結(jié)果畫圖比較，并分析步長變化對(duì)解的影響。2）當(dāng)激勵(lì)力幅值F分別按0.3，0.33，0.4，0.43，0.54，0.58，0.75，0.84，11.21，13.34進(jìn)行計(jì)算。每一個(gè)數(shù)據(jù)畫出三幅圖，分別為時(shí)間位移曲線，時(shí)間速度曲線和相圖?？疾旒?lì)力幅值F變化引起的系統(tǒng)響應(yīng)的變化。（3）請(qǐng)采用MATLAB中的內(nèi)部庫函數(shù)ode45求解此常微分方程初值問題的解，并與（1）中的結(jié)果進(jìn)行比較。四、算法描述（實(shí)驗(yàn)原理與基礎(chǔ)理論）系統(tǒng)方程和表述如下：則系統(tǒng)的輸出按如下求解：其中：大全標(biāo)準(zhǔn)文案這樣，下一個(gè)值(yn+1)由現(xiàn)在的值(yn)加上時(shí)間間隔(h)和一個(gè)估算的斜率的乘積決定。該斜率是以下斜率的加權(quán)平均：k1是時(shí)間段開始時(shí)的斜率；k2是時(shí)間段中點(diǎn)的斜率，通過歐拉法采用斜率k1來決定y在點(diǎn)tn+h/2的值；k3也是中點(diǎn)的斜率，但是這次采用斜率k2決定y值；k4是時(shí)間段終點(diǎn)的斜率，其y值用k3決定。五、實(shí)驗(yàn)內(nèi)容：求解常微分方程初值問題，考慮著名的Duffing方程。G．Duffing在1918年引入了一個(gè)帶有立方項(xiàng)的非線性振子來描述出現(xiàn)在許多力學(xué)問題中的質(zhì)量、彈簧、阻尼系統(tǒng)。從那時(shí)起，Duffing方程在非線性動(dòng)力學(xué)系統(tǒng)的研究中占有重要的地位。Duffing方程的標(biāo)準(zhǔn)形式是d2xcdxf(x)g(t)dt2dt其中：f(x)是一個(gè)含有三次項(xiàng)的非線性函數(shù)，g(t)是一個(gè)周期函數(shù)。把f(x)xx3，g(t)Fcos(t)代入上式，可得d2xcdxxx3Fcos(t)（1）dt2dt式中：c0.301,1.201,步長h0.01;初值向量為:x0=(0,0.1)。要求考察激勵(lì)力幅值F變化引起的系統(tǒng)響應(yīng)的變化。積分時(shí)間區(qū)間為：[0,50]。六、程序流程圖開始輸入a,b,n,x a,y y0大全標(biāo)準(zhǔn)文案h(ba)nK1f(x,y)hhK2f(x2,y2K1)kh,y1K3f(xhK2)22K4f(xh,yhK3)2yyh(K12K22K3K4)6xxh輸出k,x,yk=n k k 1=結(jié)束七、實(shí)驗(yàn)結(jié)果ix(i)y(i)10.00001.000020.10000.990130.20000.961540.30000.917450.40000.862160.50000.800070.60000.735380.70000.671190.80000.6098100.90000.5525大全標(biāo)準(zhǔn)文案111.00000.5000121.10000.4525131.20000.4098>>八、實(shí)驗(yàn)結(jié)果分析實(shí)驗(yàn)程序function varargout=saxplaxliu(varargin)clc,clearx0=0;xn=1.2;y0=1;h=0.1;[y,x]=lgkt4j(x0,xn,y0,h);n=length(x);fprintf( 'ix(i)y(i)\n' );fori=1:nfprintf( '%2d%4.4f%4.4f\n' ,i,x(i),y(i));endfunction z=f(x,y)z=-2*x*y^2;function [y,x]=lgkt4j(x0,xn,y0,h)x=x0:h:xn;n=