計(jì)算方法實(shí)驗(yàn)報(bào)告

[[計(jì)算方法實(shí)驗(yàn)報(bào)告]班級(jí)：信息安全09-1班學(xué)號(hào)：20092988姓名：王義濤[合肥工業(yè)大學(xué)計(jì)算機(jī)與信息學(xué)院][2011/5/28]

目錄實(shí)驗(yàn)一牛頓下山法求解非線性方程的根 3一、實(shí)驗(yàn)?zāi)康?3二、實(shí)驗(yàn)內(nèi)容 3三、基本原理 3四、算法設(shè)計(jì)與實(shí)現(xiàn) 3五、輸入與輸出 4六、源代碼 4實(shí)驗(yàn)二高斯——賽德爾法求線性方程組 6一、實(shí)驗(yàn)?zāi)康?6二、實(shí)驗(yàn)內(nèi)容 6三、算法基本原理 6四、算法設(shè)計(jì)與實(shí)現(xiàn) 6五、計(jì)算用例的參考輸出 8六、源代碼 8實(shí)驗(yàn)三高斯消去法 10一、實(shí)驗(yàn)?zāi)康?10二、實(shí)驗(yàn)內(nèi)容 10三、算法基本原理 10四、算法設(shè)計(jì)與實(shí)現(xiàn) 10五、計(jì)算用例的參考輸出 12六、源代碼 12實(shí)驗(yàn)四Romberg算法 14一、 實(shí)驗(yàn)?zāi)康?14二、 實(shí)驗(yàn)內(nèi)容 14三、算法基本原理 14四、算法設(shè)計(jì)與實(shí)現(xiàn) 14五、計(jì)算用例的參考輸出 16六、源代碼 16實(shí)驗(yàn)五RungeKutta算法 18一、實(shí)驗(yàn)?zāi)康?18二、實(shí)驗(yàn)內(nèi)容 18三、算法基本原理 18四、算法設(shè)計(jì)與實(shí)現(xiàn) 19五、計(jì)算用例的參考輸出 19六、源代碼 20感悟與體會(huì) 21實(shí)驗(yàn)一牛頓下山法求解非線性方程的根一、實(shí)驗(yàn)?zāi)康?1) 熟悉非線性方程求根簡(jiǎn)單迭代法，牛頓迭代及牛頓下山法(2) 能編程實(shí)現(xiàn)簡(jiǎn)單迭代法，牛頓迭代及牛頓下山法(3) 認(rèn)識(shí)選擇迭代格式的重要性(4) 對(duì)迭代速度建立感性的認(rèn)識(shí)；分析實(shí)驗(yàn)結(jié)果體會(huì)初值對(duì)迭代的影響二、實(shí)驗(yàn)內(nèi)容用牛頓下山法解方程（初值為0.6）輸入：初值，誤差限，迭代最大次數(shù)，下山最大次數(shù)輸出：近似根各步下山因子三、基本原理牛頓下山公式：下山因子下山條件四、算法設(shè)計(jì)與實(shí)現(xiàn)流程圖關(guān)鍵點(diǎn)：while(fabs(f(x1))>fabs(f(x0))) { cout<<"λ="<<r<<""<<endl; x1=x0-r*f(x0)/df(x0); r=r/2; k++; if(k>n){cout<<"迭代失敗"<<endl;break;} }//求解下山因子while(fabs(x0-x1)>=eps) { if(df(x0)==0){cout<<"***無法迭代***"<<endl;break;} x0=x1; x1=x0-r*f(x0)/df(x0); c++; if(c>m){cout<<"迭代失敗"<<endl;break;} }//迭代過程五、輸入與輸出 x0=0.6;e=0.0000001;m=100;n=100 輸出六、源代碼#include<iostream>#include<math.h>usingnamespacestd;doublef(doublex){ return(x*x-1)*x-1;}doubledf(doublex){ return3*x*x-1;}doublenewton(doublex0,doubleeps,intn,intm){ if(df(x0)==0)cout<<"***無法迭代***"<<endl; else { intk,c;//下山次數(shù),迭代次數(shù) doublex1=x0-f(x0)/df(x0),r=1; while(fabs(f(x1))>fabs(f(x0))) { cout<<"λ="<<r<<""<<endl; x1=x0-r*f(x0)/df(x0); r=r/2; k++; if(k>n){cout<<"迭代失敗"<<endl;break;} }//求解下山因子 while(fabs(x0-x1)>=eps) { if(df(x0)==0){cout<<"***無法迭代***"<<endl;break;} x0=x1; x1=x0-r*f(x0)/df(x0); c++; if(c>m){cout<<"迭代失敗"<<endl;break;} }//迭代過程 returnx1; } }voidmain(){ doublex0,x1,e; x0=0.6;e=0.0000001; x1=newton(x0,e,100,100); cout<<"x="<<x1<<endl;}實(shí)驗(yàn)二高斯——賽德爾法求線性方程組一、實(shí)驗(yàn)?zāi)康?1) 熟悉求解線性方程組的有關(guān)理論和方法；(2) 能編程實(shí)現(xiàn)雅可比及高斯-塞德爾迭代法；(4) 根據(jù)不同類型的方程組，選擇合適的數(shù)值方法。二、實(shí)驗(yàn)內(nèi)容用Gauss-Seidel迭代法求解方程組輸入：系數(shù)矩陣A，最大迭代次數(shù)N，初始向量，誤差限e輸出：解向量三、算法基本原理線性方程組大致分迭代法和直接法。只有收斂條件滿足時(shí)，才可以進(jìn)行迭代。高斯-塞德爾是最基本的一類迭代方法，其迭代過程中引用新值進(jìn)行剩下的計(jì)算。高斯-塞德爾迭代：四、算法設(shè)計(jì)與實(shí)現(xiàn)流程圖關(guān)鍵算法：圖4圖4.1X[i]解的向量，a[][]方程組的增廣矩陣for(k=1;k<=max_k;k++){ cout<<"第"<<k<<"次迭代"; e2=0; for(i=0;i<n;i++) { old_x=x[i]; sum=0; for(j=0;j<m;j++) if(i!=j) sum=sum+a[i][j]*x[j]; x[i]=(b[i]-sum)/a[i][i]; if(e2<fabs(old_x-x[i])) e2=fabs(old_x-x[i]); cout<<x[i]<<""; }cout<<endl; if(e2<e1) break;}五、計(jì)算用例的參考輸出六、源代碼#include<iostream>#include<cmath>usingnamespacestd;#defineMAXSIZE50voidinput(doublea[MAXSIZE][MAXSIZE],doubleb[],intn,intm){ inti,j; cout<<"請(qǐng)輸入原方程的增廣矩陣"<<endl; for(i=0;i<n;i++) { for(j=0;j<m;j++) cin>>a[i][j]; cin>>b[i]; cout<<endl; }}voidoutput(doublex[],intn){ cout<<"原方程組的解的向量為："<<endl; for(inti=0;i<n;i++) cout<<x[i]<<"";}#include"G_S.h"intmain(){ doublea[MAXSIZE][MAXSIZE],b[MAXSIZE],x[MAXSIZE]; doublee1,e2,sum,old_x; intn,m,i,j,k,max_k; cout<<"輸入原方程的行數(shù)n"<<endl; cin>>n; cout<<"輸入原方程的列數(shù)m"<<endl; cin>>m; input(a,b,n,m); cout<<"輸入迭代初始向量"<<endl; for(i=0;i<n;i++) cin>>x[i]; cout<<"輸入做大迭代次數(shù)"<<endl; cin>>max_k; cout<<"輸入誤差上限"<<endl; cin>>e1;for(k=1;k<=max_k;k++){ cout<<"第"<<k<<"次迭代"; e2=0; for(i=0;i<n;i++) { old_x=x[i]; sum=0; for(j=0;j<m;j++) if(i!=j) sum=sum+a[i][j]*x[j]; x[i]=(b[i]-sum)/a[i][i]; if(e2<fabs(old_x-x[i])) e2=fabs(old_x-x[i]); cout<<x[i]<<""; }cout<<endl; if(e2<e1) break;}if(k<=max_k) output(x,n);elsecout<<"超過最大迭代次數(shù)，迭代失?。?<<endl;}實(shí)驗(yàn)三高斯消去法一、實(shí)驗(yàn)?zāi)康?1) 熟悉求解線性方程組的有關(guān)理論和方法；(2) 能編程實(shí)現(xiàn)列主元高斯消去法；(4) 根據(jù)不同類型的方程組，選擇合適的數(shù)值方法。二、實(shí)驗(yàn)內(nèi)容用選主元高斯消去求方程組三、算法基本原理A. ；B. 消元結(jié)果直接存儲(chǔ)在系數(shù)矩陣中；C. 當(dāng)消元過程發(fā)生兩行對(duì)調(diào)的情況為偶數(shù)次時(shí)，行列式值為對(duì)角線乘積，否則為對(duì)角線乘積的相反數(shù)。四、算法設(shè)計(jì)與實(shí)現(xiàn)列主元高斯消去法：列主元消元回代列主元的約當(dāng)消去約當(dāng)消去列主元的約當(dāng)消去關(guān)鍵點(diǎn)：//選主元素 for(k=0;k<n-1;k++) { max=a[k][k];max_i=k; for(i=k+1;i<n;i++) if(fabs(a[i][k])>fabs(max)) { max=a[i][k]; max_i=i; } if(max==0)break; if(max_i!=k)//交換兩行 for(j=k;j<n+1;j++) { t=a[k][j]; a[k][j]=a[max_i][j]; a[max_i][j]=t; } for(i=k+1;i<n;i++) { a[i][k]=a[i][k]/-a[k][k]; for(j=k+1;j<n+1;j++) a[i][j]=a[i][j]+a[i][k]*a[k][j]; }//消元 } for(k=n-1;k>=0;k--) { sum=0; for(j=k+1;j<n;j++) sum=sum+a[k][j]*x[j]; x[k]=(a[k][n]-sum)/a[k][k]; }//回代五、計(jì)算用例的參考輸出六、源代碼#include<iostream>#include<cmath>usingnamespacestd;#defineMAX50voidinput(doublea[MAX][MAX+1],intn){ cout<<"輸入原方程組的增廣矩陣"<<endl; for(inti=0;i<n;i++) for(intj=0;j<n+1;j++) cin>>a[i][j];}voidoutput(doublex[],intn){ cout<<"Gauss消去法得到的原方程組的解為"<<endl; for(intk=0;k<n;k++) cout<<x[k]<<"";}#include"Gauss.h"intmain(){ doublea[MAX][MAX+1],x[MAX],sum,max,t; intn,i,j,k,max_i; cout<<"輸入原方程組的階"<<endl; cin>>n; input(a,n); for(k=0;k<n-1;k++)//選主元素 { max=a[k][k];max_i=k; for(i=k+1;i<n;i++) if(fabs(a[i][k])>fabs(max)) { max=a[i][k]; max_i=i; } if(max==0)break; if(max_i!=k)//交換兩行 for(j=k;j<n+1;j++) { t=a[k][j]; a[k][j]=a[max_i][j]; a[max_i][j]=t; } for(i=k+1;i<n;i++) { a[i][k]=a[i][k]/-a[k][k]; for(j=k+1;j<n+1;j++) a[i][j]=a[i][j]+a[i][k]*a[k][j]; }//消元 } if(max==0)cout<<"原方程組無解"<<endl; else { for(k=n-1;k>=0;k--) { sum=0; for(j=k+1;j<n;j++) sum=sum+a[k][j]*x[j]; x[k]=(a[k][n]-sum)/a[k][k]; }//回代 output(x,n); cout<<endl; } return0;}實(shí)驗(yàn)四Romberg算法實(shí)驗(yàn)?zāi)康氖煜ぬ菪芜f推算法、龍貝格算法；能編程實(shí)現(xiàn)梯形遞推算法、龍貝格算法；分析實(shí)驗(yàn)結(jié)果體會(huì)各種方法的精確度，建立計(jì)算機(jī)求解定積分問題的感性認(rèn)識(shí)實(shí)驗(yàn)內(nèi)容用龍貝格算法計(jì)算輸入：積分區(qū)間，誤差限輸出：序列Tn，Sn,Cn,Rn及積分結(jié)果（參考書本P71的表2-5）三、算法基本原理由梯形遞推公式求得梯形序列，相鄰序列值作線性組合得Simpson序列,Simpson序列作線性組合得柯特斯序列,柯特斯序列作線性組合的龍貝格序列。若|R2-R1|<e,則輸出R2;否則…依此類推。梯形遞推公式加權(quán)平均公式：四、算法設(shè)計(jì)與實(shí)現(xiàn)kTSCR0R[0][0]………………1R[1][0]R[1][1]…………2R[2][0]R[2][1]R[2][2]……3R[3][0]R[3][1]R[3][2]R[3][3]4R[4][0]R[4][1]R[4][2]R[4][3]本實(shí)驗(yàn)，未采取書本所提供的框圖算法，而是直接使用T表，先通過梯形遞推公式求解和，緊接著就求解S，通過精度判斷，是否繼續(xù)，將所得到的T，S，C，R……都存儲(chǔ)在R[][]中，然后通過加權(quán)平均公式：推出r[i+1][j]=r[i+1][j-1]+(r[i+1][j-1]-r[i][j-1])/(m-1);每次循環(huán)將m=4*m,m初值為1。關(guān)鍵點(diǎn)：先用梯形遞推法求解出T存在r[][0]中，再通過一循環(huán)求解S、C、Rh=b-a;n=1; tn=h*(f(a)+f(b))/2; cout.precision(6); do { r[i][0]=tn; s=0; x=a+h/2; while(x<b) { s=s+f(x); x=x+h; } t2n=(tn/2+h*s/2); h=h/2; tn=t2n; i++; }while(i<MAXSIZE);加權(quán)平均公式 for(j=1;j<MAXSIZE-1;j++) { m=m*4; for(i=0;i<MAXSIZE;i++) r[i+1][j]=r[i+1][j-1]+(r[i+1][j-1]-r[i][j-1])/(m-1); }五、計(jì)算用例的參考輸出六、源代碼#include<iostream>#include<cmath>#include<iomanip>usingnamespacestd;#defineMAXSIZE5doublef(doublex);intmain(){ doublex,a,b,h,s,e,tn,t2n,r[MAXSIZE][MAXSIZE]; intn,i=0,j,m=1; cout<<"請(qǐng)輸入積分區(qū)間邊界a和b"<<endl; cin>>a>>b; cout<<"輸入精度"<<endl; cin>>e; h=b-a;n=1; tn=h*(f(a)+f(b))/2; cout.precision(6); do { r[i][0]=tn; s=0; x=a+h/2; while(x<b) { s=s+f(x); x=x+h; } t2n=(tn/2+h*s/2); h=h/2; tn=t2n; i++; }while(i<MAXSIZE); for(j=1;j<MAXSIZE-1;j++) { m=m*4; for(i=0;i<MAXSIZE;i++) { r[i+1][j]=r[i+1][j-1]+(r[i+1][j-1]-r[i][j-1])/(m-1); } } cout<<"Tn、Sn、Cn、Rn如下表所示"<<endl; for(i=0;i<4;i++) { for(j=0;j<i+1;j++) cout<<r[i][j]<<" "; cout<<endl; } for(i=0;i<4;i++)cout<<r[4][i]<<" "; cout<<endl; if(fabs(r[4][3]-r[3][3])<=e) { cout<<"積分="<<r[4][3]<<endl; } elsecout<<"超出最大迭代次數(shù)??！，異常退出！"<<endl;}doublef(doublex){ if(x==0)return1; elsereturnsin(x)/x;}實(shí)驗(yàn)五RungeKutta算法一、實(shí)驗(yàn)?zāi)康?1) 熟悉數(shù)值微分中Rung-Kutta方法；(2) 能編程實(shí)現(xiàn)Rung-Kutta方法；(4)明確步長(zhǎng)對(duì)算法的影響并理解變步長(zhǎng)的Rung-Kutta方法二、實(shí)驗(yàn)內(nèi)容（1）0<x<1取h=0.1時(shí)用Rung-Kutta方法求其數(shù)值解并與精確解進(jìn)行比較。輸入：求解區(qū)間，初值，數(shù)值解個(gè)數(shù)輸出：數(shù)值解三、算法基本原理通過龍格-庫(kù)塔法我們可以獲得更高精度。經(jīng)典龍格-庫(kù)塔法即在區(qū)間[xn,xn+1]取四點(diǎn)，并對(duì)這四點(diǎn)的斜率進(jìn)行加權(quán)平均作為平均斜率，通過泰勒公式尋找使局部截?cái)嗾`差為O(h5)（即4階精度）的參數(shù)滿足條件。改進(jìn)的歐拉公式：預(yù)測(cè)校正四階（經(jīng)典）龍格-庫(kù)塔公式四、算法設(shè)計(jì)與實(shí)現(xiàn)流程圖：經(jīng)典龍格庫(kù)塔算法經(jīng)典龍格庫(kù)塔算法五、計(jì)算用例的參考輸出六、源代碼#include<iostream>#include<cmath>#include<iomanip>usingnamespacestd;doublef(doublex,doubley){ return(y-2*x/y);}voidR_K(doublex0,doubley0,doubleh,int