數(shù)值計算課程設(shè)計--典型數(shù)值算法的c++語言程序設(shè)計_第1頁
已閱讀1頁,還剩32頁未讀 繼續(xù)免費閱讀

下載本文檔

版權(quán)說明:本文檔由用戶提供并上傳,收益歸屬內(nèi)容提供方,若內(nèi)容存在侵權(quán),請進行舉報或認領(lǐng)

文檔簡介

1、<p>  數(shù)值計算課程設(shè)計說明書</p><p>  題目: 典型數(shù)值算法的C++語言程序設(shè)計 </p><p>  院 (系): 理學院 </p><p>  專業(yè)班級: 數(shù)學1xx </p><p>  學 號: </p><p>

2、  學生姓名: x </p><p>  指導教師: xxxx </p><p>  2014年7月 6日</p><p><b>  目 錄</b></p><p>  1高斯列主元法解線性方程組1</p><p><b>

3、;  1.1算法說明1</b></p><p><b>  1.2算例1</b></p><p><b>  1.3程序代碼1</b></p><p>  2牛頓法解非線性方程組2</p><p><b>  2.1算法說明2</b></p>

4、<p><b>  2.1算例:3</b></p><p><b>  2.2程序代碼5</b></p><p>  3經(jīng)典四階龍格庫塔法解一階微分方程組9</p><p><b>  3.1算法說明9</b></p><p><b>  3.2算

5、例10</b></p><p>  3.3程序代碼10</p><p>  4三次樣條插值算法(壓緊樣條)12</p><p>  4.1算法說明12</p><p><b>  4.2算例12</b></p><p>  4.3 C++程序代碼13</p>

6、<p>  4.4 matlab程序代碼16</p><p>  5龍貝格求積分算法17</p><p>  5.1算法說明17</p><p><b>  5.2算例17</b></p><p>  5.3 程序代碼17</p><p>  6 M次多項式曲線擬合18<

7、;/p><p>  6.1算法說明18</p><p><b>  6.2算例19</b></p><p>  6.3程序代碼19</p><p>  7拉格朗日插值解多項式22</p><p>  7.1算法說明22</p><p><b>  7.2算例

8、22</b></p><p>  7.3程序代碼23</p><p>  8二分法求解非線性方程24</p><p>  8.1算法說明24</p><p><b>  8.1算例24</b></p><p>  8.2程序代碼24</p><p>

9、  9 不動點迭代26</p><p>  9.1算法說明26</p><p><b>  9.2算例26</b></p><p>  9.3程序代碼26</p><p>  10復化梯形求積分公式27</p><p>  10.1算法說明27</p><p>

10、<b>  10.2算例28</b></p><p>  10.3程序代碼28</p><p><b>  11設(shè)計體會29</b></p><p>  參 考 文 獻30</p><p>  1高斯列主元法解線性方程組</p><p><b>  1.1算

11、法說明</b></p><p>  將線性方程組做成增廣矩陣,對增廣矩陣進行行變換。對第1列元素,在第i行及以下的元素中選取絕對值最大的元素,將該元素最大的行與第i行交換,然后采用高斯消元法將新得到的消去第i行以下的元素。一次進行直到。從而得到上三角矩陣。再對得到的上三角矩陣進行回代操作,即可以得到方程組的解。</p><p><b>  1.2算例</b>

12、;</p><p>  課本99頁例3.16,求解方程組的解</p><p>  運行結(jié)果如圖1-1所示。</p><p><b>  圖1-1</b></p><p><b>  1.3程序代碼</b></p><p>  #include<iostream>&

13、lt;/p><p>  #include<cmath></p><p>  using namespace std;</p><p>  const N=20;</p><p>  float a[N][N];</p><p><b>  int m;</b></p><

14、;p>  int main()</p><p>  {int i,j;</p><p>  int c,k,n,p,r;</p><p>  float x[N],l[N][N],s,d;</p><p>  cout<<"下面請輸入未知數(shù)的個數(shù)m=";</p><p><

15、b>  cin>>m;</b></p><p>  cout<<endl;</p><p>  cout<<"請按順序輸入增廣矩陣a:"<<endl;</p><p>  for(i=0;i<m;i++)</p><p><b>  {<

16、;/b></p><p>  for(j=0;j<m+1;j++)</p><p>  cin>>a[i][j];</p><p><b>  }</b></p><p>  for(i=0;i<m;i++)</p><p>  {for(j=i;j<m;j+

17、+) </p><p>  c=(fabs(a[j][i])>fabs(a[i][i]))?j:i; /*找列最大元素*/</p><p>  for(n=0;n<m+1;n++) </p><p>  {s=a[i][n]; a[i][n]=a[c][n]; a[c][n]=s;} /*將列最大數(shù)防在對角線上*/</p>

18、<p>  for(p=0;p<m+1;p++)</p><p>  cout<<a[i][p]<<"\t";</p><p>  cout<<endl;</p><p>  for(k=i+1;k<m;k++) </p><p>  {l[k][i]=a[k

19、][i]/a[i][i];</p><p>  for(r=i;r<m+1;r++) /*化成三角陣*/</p><p>  a[k][r]=a[k][r]-l[k][i]*a[i][r]; }</p><p><b>  }</b></p><p>  x[m-1]=a[

20、m-1][m]/a[m-1][m-1];</p><p>  for(i=m-2;i>=0;i--)</p><p><b>  {d=0;</b></p><p>  for(j=i+1;j<m;j++)</p><p>  d=d+a[i][j]*x[j];</p><p> 

21、 x[i]=(a[i][m]-d)/a[i][i]; /*求解*/</p><p><b>  }</b></p><p>  cout<<"該方程組的解為:"<<endl;</p><p>  for(i=0;i<m;i++)</p><p>  

22、cout<<"x["<<i<<"]="<<x[i]<<"\t"; //system("pause"); </p><p>  return 0; </p><p><b>  }</b></p><p> 

23、 2牛頓法解非線性方程組</p><p><b>  2.1算法說明</b></p><p>  設(shè)已知,第一步計算函數(shù)</p><p>  第二步計算雅可比矩陣</p><p>  第三步求線性方程組的解。第四步計算下一點重復以上過程。</p><p><b>  2.1算例:<

24、/b></p><p>  設(shè)有非線性方程組(課本137頁例3.32)</p><p>  設(shè)初始值,用牛頓法計算。</p><p>  解:迭代3次后的近似解向量為[1.9068,0.311219];具體運行結(jié)果如圖2-1,2-2,2-3所示。</p><p><b>  圖2-1</b></p>

25、<p><b>  圖2-2</b></p><p><b>  圖2-1</b></p><p><b>  2.2程序代碼</b></p><p>  #include<iostream></p><p>  #include<cmath>

26、</p><p>  #define N 2 //非線性方程組中方程的個數(shù)、未知量的個數(shù)</p><p>  #define Epsilon 0.0001 //差向量的上限</p><p>  #define Max 100 //最大迭代次數(shù)</p><p>  using namespace std;</p>

27、<p>  const int N2=2*N;</p><p>  int main()</p><p><b>  {</b></p><p>  void ff(float xx[N],float yy[N]); //計算向量函數(shù)的因變量向量yy[N]</p><p>  void ffyakebi(fl

28、oat xx[N],float yy[N][N]);//計算雅可比矩陣yy[N][N]</p><p>  void inv_yakebi(float yy[N][N],float inv[N][N]);//計算雅可比矩陣的逆矩陣inv</p><p>  void newton(float x0[N],float inv[N][N],float y0[N],float x1[N]);//

29、由近似解向量x0求近似解向量x1</p><p>  float x0[N]={2.0,0.25},y0[N],yakebi[N][N],invyakebi[N][N],x1[N],errornorm;</p><p>  int i,iter=0;</p><p>  cout<<"初始近似解向量:"<<endl;<

30、;/p><p>  for(i=0;i<N;i++)</p><p>  cin>>x0[i];</p><p>  for(i=0;i<N;i++)</p><p>  cout<<x0[i]<<" ";</p><p>  cout<<en

31、dl;cout<<endl;</p><p><b>  do</b></p><p><b>  {</b></p><p>  iter=iter+1;</p><p>  cout<<"第"<<iter<<"次迭代開

32、始"<<endl;</p><p>  ff(x0,y0);</p><p>  ffyakebi(x0,yakebi);</p><p>  inv_yakebi(yakebi,invyakebi);</p><p>  newton(x0,invyakebi,y0,x1);</p><p> 

33、 errornorm=0;</p><p>  for(i=0;i<N;i++)</p><p>  errornorm=errornorm+fabs(x1[i]);</p><p>  if(errornorm<Epsilon) break;</p><p>  for(i=0;i<N;i++)</p>&l

34、t;p>  x0[i]=x1[i];</p><p>  }while(iter<Max);</p><p><b>  return 0;</b></p><p><b>  }</b></p><p>  void ff(float xx[N],float yy[N])</p&

35、gt;<p><b>  {</b></p><p>  float x,y;</p><p><b>  int i;</b></p><p><b>  x=xx[0];</b></p><p><b>  y=yy[1];</b><

36、;/p><p>  yy[0]=x*x-2*x-y+0.5;</p><p>  yy[1]=x*x+4*y*y-4;</p><p>  cout<<"向量函數(shù)的因變量向量是:"<<endl;</p><p>  for(i=0;i<N;i++)</p><p>  co

37、ut<<yy[i]<<" ";</p><p>  cout<<endl;</p><p>  cout<<endl;</p><p><b>  }</b></p><p>  void ffyakebi(float xx[N],float yy[N]

38、[N])</p><p><b>  {</b></p><p>  float x,y;</p><p><b>  int i,j;</b></p><p><b>  x=xx[0];</b></p><p><b>  y=xx[1];

39、</b></p><p>  yy[0][0]=2*x-2;</p><p>  yy[0][1]=-1;</p><p>  yy[1][0]=2*x;</p><p>  yy[1][1]=8*y;</p><p>  cout<<"雅可比矩陣是:"<<end

40、l;</p><p>  for(i=0;i<N;i++)</p><p>  {for(j=0;j<N;j++)</p><p>  cout<<yy[i][j]<<" ";</p><p>  cout<<endl;</p><p><b&g

41、t;  }</b></p><p>  cout<<endl;</p><p><b>  }</b></p><p>  void inv_yakebi(float yy[N][N],float inv[N][N])</p><p><b>  {</b></p>

42、;<p>  float aug[N][N2],L;</p><p>  int i,j,k;</p><p>  cout<<"開始計算雅可比矩陣的逆矩陣:"<<endl;</p><p>  for(i=0;i<N;i++)</p><p>  {for(j=0;j<N

43、;j++)</p><p>  aug[i][j]=yy[i][j];</p><p>  for(j=N;j<N2;j++)</p><p>  if(j==i+N)</p><p>  aug[i][j]=1;</p><p><b>  else</b></p><

44、p>  aug[i][j]=0;</p><p><b>  }</b></p><p>  for(i=0;i<N;i++)</p><p>  {for(j=0;j<N2;j++)</p><p>  cout<<aug[i][j]<<" ";</p

45、><p>  cout<<endl;</p><p><b>  }</b></p><p>  cout<<endl;</p><p>  for(i=0;i<N;i++)</p><p><b>  {</b></p><p&

46、gt;  for(k=i+1;k<N;k++)</p><p>  {L=-aug[k][i]/aug[i][i];</p><p>  for(j=1;j<N2;j++)</p><p>  aug[k][j]=aug[k][j]+L*aug[i][j];}</p><p><b>  }</b></

47、p><p>  for(i=0;i<N;i++)</p><p>  {for(j=0;j<N2;j++)</p><p>  cout<<aug[i][j]<<" ";</p><p>  cout<<endl;}</p><p>  cout<&

48、lt;endl;</p><p>  for(i=N-1;i>0;i--)</p><p>  {for(k=i-1;k>=0;k--)</p><p>  {L=-aug[k][i]/aug[i][i];</p><p>  for(j=N2-1;j>=0;j--)</p><p>  aug[k]

49、[j]=aug[k][j]+L*aug[i][j];}</p><p><b>  }</b></p><p>  for(i=0;i<N;i++)</p><p><b>  {</b></p><p>  for(j=0;j<N2;j++)</p><p> 

50、 cout<<aug[i][j]<<" ";</p><p>  cout<<endl;</p><p><b>  }</b></p><p>  cout<<endl;</p><p>  for(i=N-1;i>=0;i--)</p&g

51、t;<p>  for(j=N2-1;j>=0;j--)</p><p>  aug[i][j]=aug[i][j]/aug[i][i];</p><p>  for(i=0;i<N;i++)</p><p><b>  {</b></p><p>  for(j=0;j<N2;j++)&

52、lt;/p><p>  cout<<aug[i][j]<<" ";</p><p>  cout<<endl;</p><p>  for(j=N;j<N2;j++)</p><p>  inv[i][j-N]=aug[i][j];</p><p><b&

53、gt;  }</b></p><p>  cout<<endl;</p><p>  cout<<"雅可比矩陣的逆矩陣是:"<<endl;</p><p>  for(i=0;i<N;i++)</p><p>  {for(j=0;j<N;j++)</p&

54、gt;<p>  cout<<inv[i][j]<<" ";</p><p>  cout<<endl;</p><p><b>  }</b></p><p>  cout<<endl;</p><p><b>  }</

55、b></p><p>  void newton(float x0[N],float inv[N][N],float y0[N],float x1[N])</p><p><b>  {</b></p><p><b>  int i,j;</b></p><p>  float sum=0;

56、</p><p>  for(i=0;i<N;i++)</p><p><b>  {sum=0;</b></p><p>  for(j=0;j<N;j++)</p><p>  sum=sum+inv[i][j]*y0[j];</p><p>  x1[i]=x0[i]-sum;&

57、lt;/p><p><b>  }</b></p><p>  cout<<"近似解向量:"<<endl;</p><p>  for(i=0;i<N;i++)</p><p>  cout<<x1[i]<<" ";</p&g

58、t;<p>  cout<<endl;</p><p>  cout<<endl;</p><p><b>  }</b></p><p>  3經(jīng)典四階龍格庫塔法解一階微分方程組</p><p><b>  3.1算法說明</b></p><

59、;p>  4階龍格-庫塔方法(RK4)可模擬N=4的泰勒方法的精度。這種算法可以描述為,自初始點開始,利用下面的計算方法生成近似序列</p><p><b>  3.2算例:</b></p><p>  用4階龍格-庫塔方法計算區(qū)間[0.0,0.2]上式(3)的數(shù)值解,采用10個子區(qū)間,步長h=0.02。(課本400頁例9.15)求解運行結(jié)果如圖3-1所示。&l

60、t;/p><p><b>  初值</b></p><p><b>  圖3-1</b></p><p><b>  3.3程序代碼</b></p><p>  #include<iostream></p><p>  #include<i

61、omanip></p><p>  using namespace std;</p><p>  void LG(double (*f)(double t,double x,double y),</p><p>  double (*g)(double t,double x,double y),</p><p>  double chu

62、zhi[3],double resu[3],double h)</p><p><b>  {</b></p><p>  double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1;</p><p>  t0=chuzhi[0];x0=chuzhi[1];y0=chuzhi[2];</p><

63、;p>  f1=f(t0,x0,y0);</p><p>  g1=g(t0,x0,y0);</p><p>  f2=f(t0+h/2,x0+h*f1/2,y0+h*g1/2);</p><p>  g2=g(t0+h/2,x0+h*f1/2,y0+h*g1/2);</p><p>  f3=f(t0+h/2,x0+h*f2/2,y0

64、+h*g2/2);</p><p>  g3=g(t0+h/2,x0+h*f2/2,y0+h*g2/2);</p><p>  f4=f(t0+h,x0+h*f3,y0+h*g3);</p><p>  g4=f(t0+h,x0+h*f3,y0+h*g3);</p><p>  x1=x0+h*(f1+2*f2+2*f3+f4)/6;<

65、/p><p>  y1=y0+h*(g1+2*g2+2*g3+g4)/6;</p><p>  resu[0]=t0+h;resu[1]=x1;resu[2]=y1;</p><p><b>  }</b></p><p>  int main()</p><p><b>  {</b

66、></p><p>  double f(double t,double x,double y);</p><p>  double g(double t,double x,double y);</p><p>  double chuzhi[3],resu[3];</p><p>  double a,b,S;</p>

67、<p>  double t,step;</p><p><b>  int i;</b></p><p>  cout<<"輸入微分方程的初值:";</p><p>  cin>>chuzhi[0]>>chuzhi[1]>>chuzhi[2];</p>

68、<p>  cout<<"輸入所求微分方程組的微分區(qū)間[a,b]:";</p><p>  cin>>a>>b;</p><p>  cout<<"輸入所求微分方程組所分解子區(qū)間的個數(shù)step:";</p><p>  cin>>step;</p&g

69、t;<p>  S=(b-a)/step;</p><p>  cout<<chuzhi[0]<<setw(10)<<chuzhi[1]<<setw(10)<<chuzhi[2]<<endl;</p><p>  for(i=0;i<step;i++)</p><p><

70、;b>  {</b></p><p>  LG(f,g,chuzhi,resu,S);</p><p>  cout<<resu[0]<<setw(10)<<resu[1]<<setw(10)<<resu[2]<<endl;</p><p>  chuzhi[0]=resu[0

71、];chuzhi[1]=resu[1];chuzhi[2]=resu[2];</p><p><b>  }</b></p><p><b>  return 0;</b></p><p><b>  }</b></p><p>  double f(double t,doub

72、le x,double y)</p><p><b>  {</b></p><p>  double dx;</p><p><b>  dx=x+2*y;</b></p><p>  return(dx);</p><p><b>  }</b>&l

73、t;/p><p>  double g(double t,double x,double y)</p><p><b>  {</b></p><p>  double dy;</p><p>  dy=3*x+2*y;</p><p>  return(dy);</p><p&g

74、t;<b>  }</b></p><p>  4三次樣條插值算法(壓緊樣條)</p><p><b>  4.1算法說明</b></p><p>  表5-1 三次樣條插值算法說明表</p><p><b>  4.2算例</b></p><p>  

75、求三次緊壓樣條曲線,以經(jīng)過點(0,0),(1,0.5),(2,2.0)和(3,1.5),且一階導數(shù)的邊界條件為S’(0)=0.2和S’(3)=-1。(課本222頁例5.7)</p><p>  求解運行結(jié)果界面如下:</p><p>  圖4-1 三次樣條插值算法(壓緊樣條)</p><p>  Matlab作圖如圖4-2所示。</p><p&g

76、t;<b>  圖4-2</b></p><p>  4.3 C++程序代碼</p><p>  #include<iostream></p><p>  #include<iomanip></p><p>  using namespace std;</p><p>  

77、const int max = 50;</p><p>  float x[max], y[max], h[max];</p><p>  float c[max], a[max], fxym[max];</p><p>  float f(int x1, int x2, int x3)</p><p><b>  {</b&

78、gt;</p><p>  float a=(y[x3]-y[x2])/(x[x3]-x[x2]);</p><p>  float b=(y[x2]-y[x1])/(x[x2]-x[x1]);</p><p>  return (a-b)/(x[x3]-x[x1]);</p><p>  } //求差分</p><

79、p>  void cal_m(int n)//用追趕法求解出彎矩向量M</p><p>  { </p><p>  float B[max];</p><p>  B[0] = c[0] / 2;</p><p>  for(int i = 1; i < n; i++)</p>&

80、lt;p>  B[i] = c[i] / (2 - a[i]*B[i-1]);</p><p>  fxym[0] = fxym[0] / 2;</p><p>  for(i = 1; i <= n; i++)</p><p>  fxym[i] = (fxym[i] - a[i]*fxym[i-1]) / (2 - a[i]*B[i-1]);<

81、/p><p>  for(i = n-1; i >= 0; i--)</p><p>  fxym[i] = fxym[i] - B[i]*fxym[i+1];</p><p><b>  }</b></p><p>  void printout(int n);</p><p>  int m

82、ain()</p><p><b>  {</b></p><p>  int n,i; char ch;</p><p>  do{ cout<<"輸入x的最大下標:";</p><p><b>  cin>>n;</b></p><

83、p>  for(i = 0; i <= n; i++)</p><p><b>  {</b></p><p>  cout<<"Please put in X"<<i<<':';</p><p>  cin>>x[i]; </p>

84、<p>  cout<<"Please put in Y"<<i<<':';</p><p>  cin>>y[i]; </p><p><b>  }</b></p><p>  for(i = 0; i < n; i++)

85、 //求步長</p><p>  h[i] = x[i+1] - x[i];</p><p>  cout<<"Please 輸入邊界條件\n 1 :已知兩端的一階導數(shù)\n 2 :兩端的二階導數(shù)已知\n";</p><p><b>  int t;</b></p><p> 

86、 float f0, f1;</p><p><b>  cin>>t;</b></p><p><b>  switch(t)</b></p><p><b>  {</b></p><p>  case 1:cout<<"輸入 Y0\

87、9; Y"<<n<<"\'\n";</p><p>  cin>>f0>>f1;</p><p>  c[0] = 1; a[n] = 1;</p><p>  fxym[0] = 6*((y[1] - y[0]) / (x[1] - x[0]) - f0) / h[0];<

88、/p><p>  fxym[n] = 6*(f1 - (y[n] - y[n-1]) / (x[n] - x[n-1])) / h[n-1];</p><p><b>  break;</b></p><p>  case 2:cout<<"輸入 Y0\" Y"<<n<<"

89、\"\n";</p><p>  cin>>f0>>f1;</p><p>  c[0] = a[n] = 0;</p><p>  fxym[0] = 2*f0; fxym[n] = 2*f1;</p><p><b>  break;</b></p><

90、p>  default:cout<<"不可用\n";//待定</p><p><b>  };</b></p><p>  for(i = 1; i < n; i++)</p><p>  fxym[i] = 6 * f(i-1, i, i+1);</p><p>  for(

91、i = 1; i < n; i++)</p><p><b>  {</b></p><p>  a[i] = h[i-1] / (h[i] + h[i-1]);</p><p>  c[i] = 1 - a[i];</p><p><b>  }</b></p><p&g

92、t;  a[n] = h[n-1] / (h[n-1] + h[n]);</p><p><b>  cal_m(n);</b></p><p>  cout<<"\n輸出三次樣條插值函數(shù):\n";</p><p>  printout(n); </p><p>  cout<<

93、;"Do you have another try ? y/n :";</p><p><b>  cin>>ch;</b></p><p>  }while(ch == 'y' || ch == 'Y');</p><p>  return 0;}</p><

94、p>  void printout(int n)</p><p>  {cout<<setprecision(6);</p><p>  for(int i = 0; i < n; i++)</p><p>  {cout<<i+1<<": ["<<x[i]<<"

95、, "<<x[i+1]<<"]\n"<<"\t";</p><p>  cout<<"S"<<i+1<<"=";</p><p>  float t = fxym[i]/(6*h[i]);</p><p>

96、  if(t > 0) cout<<-t<<"*(x - "<<x[i+1]<<")^3";</p><p>  else cout<<-t<<"*(x - "<<x[i+1]<<")^3";</p><p>

97、  t = fxym[i+1]/(6*h[i]);</p><p>  if(t > 0) cout<<" + "<<t<<"*(x - "<<x[i]<<")^3";</p><p>  else cout<<" - "<&l

98、t;t<<"*(x - "<<x[i]<<")^3";</p><p>  cout<<"\n\t";</p><p>  t = (y[i] - fxym[i]*h[i]*h[i]/6)/h[i];</p><p>  if(t > 0) cout&l

99、t;<"- "<<t<<"*(x - "<<x[i+1]<<")";</p><p>  else cout<<"- "<<-t<<"*(x - "<<x[i+1]<<")";<

100、;/p><p>  t = (y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i];</p><p>  if(t > 0) cout<<" + "<<t<<"*(x - "<<x[i]<<")";</p><p>  els

101、e cout<<" - "<<-t<<"*(x - "<<x[i]<<")";</p><p>  cout<<endl; }</p><p>  cout<<endl;</p><p><b>  }</

102、b></p><p>  4.4 matlab程序代碼</p><p>  x1=0:0.01:1;</p><p>  y1=0.06*(x1 - 1).^3 + 0.42*(x1 - 0).^3 - 0.06*(x1 - 1) + 0.08*(x1 - 0);</p><p>  x2=1:0.01:2;</p>&l

103、t;p>  y2=-0.42*(x2 - 2).^3 - 0.62*(x2 - 1).^3 - 0.08*(x2 - 2) + 2.62*(x2 - 1);</p><p>  x3=2:0.01:3;</p><p>  y3=0.62*(x3 - 3).^3 + 0.06*(x3 - 2).^3 - 2.62*(x3 - 3) + 1.44*(x3 - 2);</p>

104、<p>  X=[0 1 2 3];</p><p>  Y=[0 0.5 2 1.5];</p><p>  plot(x1,y1,x2,y2,x3,y3,X,Y,'*')</p><p>  gtext('S1')</p><p>  gtext('S2')</p>

105、<p>  gtext('S3')</p><p><b>  5龍貝格求積分算法</b></p><p><b>  5.1算法說明</b></p><p>  生成的逼近表,并以為最終解來逼近積分</p><p>  逼近存在于一個特別的下三角矩陣中,第0列元素用基

106、于個[a,b]子區(qū)間的連續(xù)梯形方法計算,然后利用龍貝格公式計算。當時,第行的元素為</p><p>  當時,程序在第行結(jié)束。</p><p><b>  5.2算例</b></p><p>  用龍貝格求積公式計算的值,運行結(jié)果如圖5-1所示。</p><p><b>  圖5-1</b><

107、/p><p><b>  5.3 程序代碼</b></p><p>  #include<iostream></p><p>  #include<stdio.h></p><p>  #include<math.h></p><p>  using namespa

108、ce std;</p><p>  #define f(x) (sin(x)) //列舉函數(shù)</p><p>  #define N 20 //區(qū)間等分數(shù)</p><p>  #define MAX 10 //最大迭代次數(shù)</p><p>  #define a 0 //所求積分的上下限</p><p>

109、;  #define b 1</p><p>  #define epsilon 0.0001 //精度</p><p>  double Romberg(double aa,double bb,long int n)</p><p><b>  {</b></p><p><b>  int i;</b

110、></p><p>  double sum,h=(bb-aa)/n;sum=0;</p><p>  for(i=1;i<n;i++)</p><p>  sum=sum+f(aa+i*h);</p><p>  sum=sum+(f(aa)+f(bb))/2;</p><p>  return (h*s

111、um);</p><p><b>  }</b></p><p>  void main()</p><p><b>  {</b></p><p><b>  int i;</b></p><p>  long int n=N,m=0;</p&g

112、t;<p>  double T[2][MAX+1];</p><p>  T[1][0]=Romberg(a,b,n);</p><p><b>  n=n*2;</b></p><p>  for(m=1;m<MAX;m++)</p><p><b>  {</b></

113、p><p>  for(i=0;i<=m;i++)</p><p>  {T[0][i]=T[1][i];}</p><p>  T[1][0]=Romberg(a,b,n);</p><p><b>  n=n*2;</b></p><p>  for(i=1;i<=m;i++)<

114、/p><p>  T[1][i]=T[1][i-1]+(T[1][i-1]-T[0][i-1])/(pow(2,2*m)-1);</p><p>  if((T[0][m-1]-T[1][m])<epsilon)</p><p>  cout<<"T="<<T[1][m]<<endl;</p>

115、<p><b>  return;</b></p><p><b>  }</b></p><p><b>  }</b></p><p>  6 M次多項式曲線擬合</p><p><b>  6.1算法說明</b></p>&l

116、t;p>  設(shè)有N個點,橫坐標是確定的。最小二乘拋物線的系數(shù)表示為(6-1)求解A,B和C的線性方程組為</p><p><b>  6.2算例</b></p><p>  根據(jù)4個點(-3,3),(0,1),(2,1)和(4,3),求解最小二乘拋物線。(課本211頁例5.6)求解結(jié)果如圖6-1所示。</p><p><b> 

117、 圖6-1</b></p><p><b>  6.3程序代碼</b></p><p>  #include<iostream></p><p>  #include<cmath></p><p>  #define MAX 20</p><p>  using

118、 namespace std;</p><p>  int main()</p><p><b>  {</b></p><p>  int n,m,i,j,k;</p><p>  void inv(double X[MAX][MAX],int n,double E[MAX][MAX]);</p><

119、;p>  double X[MAX]={0},Y[MAX]={0},Z[MAX][MAX]={0},B[MAX]={0},C[MAX]={0};</p><p>  double A[MAX][MAX]={0},BZ[MAX][MAX]={0},E[MAX][MAX]={0};</p><p>  cout<<"M次多項式曲線擬合\n請輸入待擬合的點的個數(shù):&q

120、uot;;</p><p><b>  cin>>n;</b></p><p>  cout<<endl;</p><p>  cout<<"請輸入"<<n<<"個點的X坐標序列:";</p><p>  for(i=0;

121、i<n;i++)</p><p>  cin>>X[i];</p><p>  cout<<endl;</p><p>  cout<<"請輸入"<<n<<"個點的Y坐標序列:";</p><p>  for(i=0;i<n;i++

122、)</p><p>  cin>>Y[i];</p><p>  cout<<endl;</p><p>  cout<<"請輸入需要擬合的次數(shù):";</p><p><b>  cin>>m;</b></p><p>  for

123、(i=0;i<n;i++)</p><p>  for(k=1;k<=m+1;k++)</p><p>  Z[i][k-1]=pow(X[i],k-1);</p><p><b>  //求Z的轉(zhuǎn)置</b></p><p>  for(i=0;i<n;i++)</p><p>

124、<b>  {</b></p><p>  for(j=0;j<m+1;j++)</p><p><b>  {</b></p><p>  BZ[j][i]=Z[i][j];</p><p><b>  }</b></p><p><b&g

125、t;  } </b></p><p>  //計算其轉(zhuǎn)置的BF與Z的乘 </p><p>  for(i=0;i<m+1;i++)</p><p>  for(j=0;j<m+1;j++)</p><p>  for(k=0;k<n;k++)</p><p>  A[i][j]+=BZ[i

126、][k]*Z[k][j];//計算Z的轉(zhuǎn)置BZ與Y的乘</p><p>  for(i=0;i<m+1;i++)</p><p>  for(j=0;j<n;j++)</p><p>  B[i]+=BZ[i][j]*Y[j];//調(diào)用inv函數(shù)求解矩陣A的逆矩陣E</p><p>  inv(A,n,E);//計算A的逆BZ與B

127、的乘</p><p>  for(i=0;i<m+1;i++)</p><p>  for(j=0;j<n;j++)</p><p>  C[i]+=E[i][j]*B[j];</p><p>  cout<<endl;</p><p>  cout<<"擬合后的"

128、;<<m<<"次多項式系數(shù)為,冪次由高到低:"<<endl;</p><p>  for(i=m;i>=0;i--)</p><p>  cout<<C[i]<<"\t";</p><p>  cout<<endl;</p><p

129、>  cout<<"擬合后的"<<m<<"次多項式為:"<<endl;</p><p>  cout<<"P(x)=";</p><p>  for(i=m;i>=0;i--)</p><p><b>  {</b>

130、;</p><p><b>  if(i==0)</b></p><p>  cout<<"+"<<C[i];</p><p><b>  else</b></p><p>  cout<<"+"<<C[i]&l

131、t;<"*x^"<<i;</p><p><b>  }</b></p><p>  cout<<endl;</p><p><b>  return 0;</b></p><p><b>  }</b></p>

132、<p>  void inv(double X[MAX][MAX],int n,double E[MAX][MAX])</p><p>  //求解任意可逆矩陣的逆,X為待求解矩陣,E為全零矩陣,非單位矩陣,也可以是單位矩陣</p><p><b>  {</b></p><p>  int i,j,k;</p><

133、;p>  double temp=0;</p><p>  for(i=0;i<MAX;i++)</p><p><b>  {</b></p><p>  for(j=0;j<MAX;j++)</p><p><b>  if(i==j)</b></p><p

134、>  E[i][j]=1;</p><p><b>  }</b></p><p>  for(i=0;i<n-1;i++)</p><p><b>  {</b></p><p>  temp=X[i][i];</p><p>  for(j=0;j<n;

135、j++)</p><p><b>  {</b></p><p>  X[i][j]=X[i][j]/temp;</p><p>  E[i][j]=E[i][j]/temp;</p><p><b>  }</b></p><p>  for(k=0;k<n;k++)

136、</p><p><b>  {</b></p><p><b>  if(k==i)</b></p><p><b>  continue;</b></p><p>  temp=-X[i][i]*X[k][i];</p><p>  for(j=0;

137、j<n;j++)</p><p><b>  {</b></p><p>  X[k][j]=X[k][j]+temp*X[i][j];</p><p>  E[k][j]=E[k][j]+temp*E[i][j];</p><p><b>  }</b></p><p&g

138、t;<b>  }</b></p><p><b>  }</b></p><p><b>  }</b></p><p>  7拉格朗日插值解多項式</p><p><b>  7.1算法說明</b></p><p><b

139、>  7.2算例</b></p><p>  給定函數(shù)四個點坐標(1.1,3.887),(2.3,4.276),(3.9,4.651),(5.1,2.117),試用拉格朗日插值確定函數(shù)在x=2.101處的值。</p><p>  求解運行結(jié)果界面如圖7-1所示,可知函數(shù)在點2.101處的值為4.14569。</p><p><b>  圖

140、7-1</b></p><p><b>  7.3程序代碼</b></p><p>  #include<iostream></p><p>  #define N 4//插值節(jié)點的個數(shù)</p><p>  using namespace std;</p><p>  vo

141、id main()</p><p><b>  {</b></p><p>  float x[N],y[N],a,f=0,temp;</p><p><b>  int i,j;</b></p><p>  cout<<"輸入插值節(jié)點的坐標:"<<end

142、l;</p><p>  for(i=0;i<N;i++)</p><p>  {cin>>x[i];</p><p>  cin>>y[i];</p><p><b>  }</b></p><p>  cout<<"輸入所求點橫坐標:&qu

143、ot;;</p><p><b>  cin>>a;</b></p><p>  for(i=0;i<N;i++)</p><p><b>  {</b></p><p><b>  temp=1;</b></p><p><b

144、>  {</b></p><p>  for(j=0;j<N;j++)</p><p><b>  if(i!=j)</b></p><p>  temp=temp*(a-x[j])/(x[i]-x[j]);</p><p><b>  }</b></p>&l

145、t;p>  f=f+temp*y[i];</p><p><b>  }</b></p><p>  cout<<"所求值為:"<<f<<endl;</p><p><b>  }</b></p><p>  8二分法求解非線性方程<

146、;/p><p><b>  8.1算法說明</b></p><p>  若要求已知函數(shù)的根(x的解),則第一步先找出一個區(qū)間[a,b],使得與異號,因此這個區(qū)間內(nèi)一定包含著方程的根。第二步求該區(qū)間的中點 ,并找出的值。若與正負號相同則取[m,b]為新的區(qū)間, 否則取[a,m],重復第二和第三步至理想精確度為止。</p><p><b> 

147、 8.1算例</b></p><p>  利用二分法尋找函數(shù)在區(qū)間[0,2]內(nèi)的的零點。(課本44頁例2.7)運行結(jié)果如圖8-1所示。</p><p><b>  圖8-1</b></p><p><b>  8.2程序代碼</b></p><p>  #include <iost

148、ream.h></p><p>  #include <math.h></p><p>  #define eps 0.001</p><p>  double fun(double x)</p><p><b>  {</b></p><p>  return x*sin(

149、x)-1;</p><p><b>  }</b></p><p>  double dichotomy(double a,double b)</p><p><b>  {</b></p><p>  double c=0.0;</p><p>  if((fun(a)&l

150、t;0)&&(fun(b)>0))</p><p><b>  {</b></p><p>  while(true)</p><p><b>  {</b></p><p>  c=(a+b)/2;</p><p>  if(fun(c)<0)&

151、lt;/p><p><b>  {</b></p><p><b>  a=c;</b></p><p>  if(fabs((a-b))<eps)</p><p><b>  {</b></p><p>  return (a+b)/2;</p

溫馨提示

  • 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
  • 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯(lián)系上傳者。文件的所有權(quán)益歸上傳用戶所有。
  • 3. 本站RAR壓縮包中若帶圖紙,網(wǎng)頁內(nèi)容里面會有圖紙預覽,若沒有圖紙預覽就沒有圖紙。
  • 4. 未經(jīng)權(quán)益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
  • 5. 眾賞文庫僅提供信息存儲空間,僅對用戶上傳內(nèi)容的表現(xiàn)方式做保護處理,對用戶上傳分享的文檔內(nèi)容本身不做任何修改或編輯,并不能對任何下載內(nèi)容負責。
  • 6. 下載文件中如有侵權(quán)或不適當內(nèi)容,請與我們聯(lián)系,我們立即糾正。
  • 7. 本站不保證下載資源的準確性、安全性和完整性, 同時也不承擔用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。

評論

0/150

提交評論