計算方法實驗報告

[[計算方法實驗報告]班級：信息安全09-1班學號：20092988姓名：王義濤[合肥工業(yè)大學計算機與信息學院][2011/5/28]

目錄實驗一牛頓下山法求解非線性方程的根 3一、實驗目的 3二、實驗內容 3三、基本原理 3四、算法設計與實現(xiàn) 3五、輸入與輸出 4六、源代碼 4實驗二高斯——賽德爾法求線性方程組 6一、實驗目的 6二、實驗內容 6三、算法基本原理 6四、算法設計與實現(xiàn) 6五、計算用例的參考輸出 8六、源代碼 8實驗三高斯消去法 10一、實驗目的 10二、實驗內容 10三、算法基本原理 10四、算法設計與實現(xiàn) 10五、計算用例的參考輸出 12六、源代碼 12實驗四Romberg算法 14一、 實驗目的 14二、 實驗內容 14三、算法基本原理 14四、算法設計與實現(xiàn) 14五、計算用例的參考輸出 16六、源代碼 16實驗五RungeKutta算法 18一、實驗目的 18二、實驗內容 18三、算法基本原理 18四、算法設計與實現(xiàn) 19五、計算用例的參考輸出 19六、源代碼 20感悟與體會 21實驗一牛頓下山法求解非線性方程的根一、實驗目的(1) 熟悉非線性方程求根簡單迭代法，牛頓迭代及牛頓下山法(2) 能編程實現(xiàn)簡單迭代法，牛頓迭代及牛頓下山法(3) 認識選擇迭代格式的重要性(4) 對迭代速度建立感性的認識；分析實驗結果體會初值對迭代的影響二、實驗內容用牛頓下山法解方程（初值為0.6）輸入：初值，誤差限，迭代最大次數(shù)，下山最大次數(shù)輸出：近似根各步下山因子三、基本原理牛頓下山公式：下山因子下山條件四、算法設計與實現(xià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;}實驗二高斯——賽德爾法求線性方程組一、實驗目的(1) 熟悉求解線性方程組的有關理論和方法；(2) 能編程實現(xiàn)雅可比及高斯-塞德爾迭代法；(4) 根據(jù)不同類型的方程組，選擇合適的數(shù)值方法。二、實驗內容用Gauss-Seidel迭代法求解方程組輸入：系數(shù)矩陣A，最大迭代次數(shù)N，初始向量，誤差限e輸出：解向量三、算法基本原理線性方程組大致分迭代法和直接法。只有收斂條件滿足時，才可以進行迭代。高斯-塞德爾是最基本的一類迭代方法，其迭代過程中引用新值進行剩下的計算。高斯-塞德爾迭代：四、算法設計與實現(xià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;}五、計算用例的參考輸出六、源代碼#include<iostream>#include<cmath>usingnamespacestd;#defineMAXSIZE50voidinput(doublea[MAXSIZE][MAXSIZE],doubleb[],intn,intm){ inti,j; cout<<"請輸入原方程的增廣矩陣"<<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;}實驗三高斯消去法一、實驗目的(1) 熟悉求解線性方程組的有關理論和方法；(2) 能編程實現(xiàn)列主元高斯消去法；(4) 根據(jù)不同類型的方程組，選擇合適的數(shù)值方法。二、實驗內容用選主元高斯消去求方程組三、算法基本原理A. ；B. 消元結果直接存儲在系數(shù)矩陣中；C. 當消元過程發(fā)生兩行對調的情況為偶數(shù)次時，行列式值為對角線乘積，否則為對角線乘積的相反數(shù)。四、算法設計與實現(xià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]; }//回代五、計算用例的參考輸出六、源代碼#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;}實驗四Romberg算法實驗目的熟悉梯形遞推算法、龍貝格算法；能編程實現(xiàn)梯形遞推算法、龍貝格算法；分析實驗結果體會各種方法的精確度，建立計算機求解定積分問題的感性認識實驗內容用龍貝格算法計算輸入：積分區(qū)間，誤差限輸出：序列Tn，Sn,Cn,Rn及積分結果（參考書本P71的表2-5）三、算法基本原理由梯形遞推公式求得梯形序列，相鄰序列值作線性組合得Simpson序列,Simpson序列作線性組合得柯特斯序列,柯特斯序列作線性組合的龍貝格序列。若|R2-R1|<e,則輸出R2;否則…依此類推。梯形遞推公式加權平均公式：四、算法設計與實現(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]本實驗，未采取書本所提供的框圖算法，而是直接使用T表，先通過梯形遞推公式求解和，緊接著就求解S，通過精度判斷，是否繼續(xù)，將所得到的T，S，C，R……都存儲在R[][]中，然后通過加權平均公式：推出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。關鍵點：先用梯形遞推法求解出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);加權平均公式 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); }五、計算用例的參考輸出六、源代碼#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ū)間邊界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;}實驗五RungeKutta算法一、實驗目的(1) 熟悉數(shù)值微分中Rung-Kutta方法；(2) 能編程實現(xiàn)Rung-Kutta方法；(4)明確步長對算法的影響并理解變步長的Rung-Kutta方法二、實驗內容（1）0<x<1取h=0.1時用Rung-Kutta方法求其數(shù)值解并與精確解進行比較。輸入：求解區(qū)間，初值，數(shù)值解個數(shù)輸出：數(shù)值解三、算法基本原理通過龍格-庫塔法我們可以獲得更高精度。經(jīng)典龍格-庫塔法即在區(qū)間[xn,xn+1]取四點，并對這四點的斜率進行加權平均作為平均斜率，通過泰勒公式尋找使局部截斷誤差為O(h5)（即4階精度）的參數(shù)滿足條件。改進的歐拉公式：預測校正四階（經(jīng)典）龍格-庫塔公式四、算法設計與實現(xiàn)流程圖：經(jīng)典龍格庫塔算法經(jīng)典龍格庫塔算法五、計算用例的參考輸出六、源代碼#include<iostream>#include<cmath>#include<iomanip>usingnamespacestd;doublef(doublex,doubley){ return(y-2*x/y);}voidR_K(doublex0,doubley0,doubleh,int