《計算方法》實驗指導書

1、計算方法實驗指導書實驗1方程求根、實驗目的1. 通過對二分法、牛頓法、割線法作編程練習，進一步體會它們各自不同的特點;2. 了解二分法，切線法，割線法。3. 能熟練運用二分法，牛頓法進行方程求根4. 通過上機調試運行，對方程求根的幾種方法程序進行改進。、實驗要求1. 上機前作好充分準備，包括復習編程所需要的語言工具。2. 上機時要遵守實驗室的規(guī)章制度，愛護實驗設備。3. 記錄調試過程及結果，記錄并比較與手工運算結果的異同。4. 程序調試完后，須由實驗輔導教師在機器上檢查運行結果。5. 給出本章實驗單元的實驗報告。三、實驗環(huán)境、設備1.硬件設備：IBM PC以上計算 機，有硬盤和一個軟驅、單機

2、和網絡環(huán)境均可。2 .軟件環(huán)境：C語言運行環(huán)境。四、實驗原理、方法二分算法計算步驟：(1)輸入有根區(qū)間的端點a、b及預先給定的精度£；(2 )計算中點 x=(a+b)/2 ;(3) 若 f(x)f(b)<0 ,則 a=x,轉向下一 步；否則b=x，轉向下一步；(4) 若b-a< £ ,則輸出方程滿足精度 要求的根x,結束；否則轉向步驟(2)。 迭代法：圖2J二分法框團牛頓法：f(x)=O的求根問題歸結為計算一系牛頓迭代法是一種逐步線性化方法，即將非線性方程 列線性方程的根。設xk是方程f(x)=O的一個近似根，將f(x)在xk處作一階泰勒展開，即 f(x)f(+

3、f ' k(x<- x k)于是得到如下的近似方程(2.7)f(x k)+f ' k(X<- x k)=0設f ' kx工0則式(2.7)的解為X =Xkf (Xk)f'(Xk)取X作為原方程的新的近似根Xk+1，即令_ f(Xk)xk 1 xk 'f (Xk)k=0,1,2,-(2.8 )則稱式(2.8)為牛頓迭代公式。用牛頓迭代公式( 簡稱牛頓法，又稱切線法。2.8)求方程近似根的方法稱為牛頓迭代法,五、實驗內容1.以方程：x3-0.2x2-0.2x-1.2=0 為例,編寫程序求方程的根2 編寫二分法、迭代法、牛頓法程序，分析運行結果。3

4、 對用這兩種方法求解出的根進行對比分析六、實驗步驟1 根據實驗題目，給出題目的 C 程序。2 上機輸入和調試自己所編的程序。3 上機結束后，應整理出實驗報告。七、實驗報告要求及記錄、格式按金陵科技學院實驗報告（工科） 格式填寫附 1 ：牛頓法程序核心部分：for(i=0;i<N;i+)printf("x(%d)=%fn",i,x1); x1=x0-f(x0)/f1(x0); /* 牛頓迭代 */ if(fabs(x1-x0)<epsilon|fabs(f(x1)<epsilon)printf("n The root of the equation

5、 is x=%fn",x1);/* 滿足精度，輸出近似根 */ return;x0=x1;mk =aikk)/a：k出)_(k)aij=誦-mik akjb(f-mkbkk)計算n(3)回代過程i,j=k+1,k+2，(n) bnnn= (bi(i)高斯消去法框圖實驗2線性方程組數值解法、實驗目的1 掌握方程組的解法，迭代法及其收斂性。2 .能熟練掌握高斯消去法，列主元高斯消去法，三角分解法。3 掌握雅可比迭代法，高斯=賽德爾迭代求線性方程組的解。、實驗要求1上機前作好充分準備，比較不用的方法解決相同問題的不同。2 上機時要遵守實驗室的規(guī)章制度，愛護實驗設備。3 記錄調試過程及結果，

6、記錄并比較與手工 運算結果的異同。4 程序調試完后，須由實驗輔導教師在機器 上檢查運行結果。5 給出本章實驗單元的實驗報告。三、實驗設備、環(huán)境1 .硬件設備：IBM PC以上計算機，有硬盤和 一個軟驅、單機和網絡環(huán)境均可。2 .軟件環(huán)境：C語言運行環(huán)境。四、實驗原理、方法1、高斯消去法：1 )計算步驟(1)輸入方程組的階數n,系數矩陣A和右端常數矩陣b消元過程：設 ak=0 ,對k=1 , 2,，n-1na。Xj)/a(i),i = n -1,2,1j -i I(4)輸岀方程組的解2、列主元高斯消去法(1)、輸入方程組的階數n ,系數矩陣A和右端常數矩陣；.(kJ) (k J)9 :, a、列

7、主元素：對k=1 , 2,，n-1,選岀；akkk41,k:八，ank(k J)1中絕對值最大的元素Ek =a絆 /ad=akkL=kj)aij=a(k) mik akjk)K(k+)以k) 皿以k)i=k+1bi mik bkj從主程序來對k行和m行交換后，再作第 k步消元操作。(3)、消元過程：對 k=1 , 2 ,n-1計算(i,j=k+1, k+2，n)(4)、回代過程Xn(n) bnXin-aij )xj) /aii J (i - n 11,2,1)j ± 1(5)、輸岀方程組的解3、三角分解法:(1)根據方程組得到增廣矩陣(2)對 j=1 , 2,n計算u1 j - a1

8、j對 i=2 , 3,n計算li1ai1U11(3)對 k=1 , 2,-a.對 j=k,k+1,n+1 計算 pk Ab.對 i=k+1 ,k+2 ,n 計算 l ik = (aik hquqk)/ ukk qmk Jl kqu qjq m(4)回代計算解Xn =yn/Unn,對 k=n-1 , n-2 ,2 , 1 計算nXk =(yk 一 ' UkqXq)/Ukkq »十(kJ) am,k否是d=aikL=i否i=i+1=n?是是.d=0?打印奇標志是L=k?否t=aij,aij=akj ,akj=t j=k,k+1,n t=bk,bk = bL,bL=t返回主程序列主

9、元框圖結束五、實驗內容1 求解方程組:40.24(2)0.0930.04-0.08-0.08 xi-0.15 X24 JLX389'10J2 編寫高斯消去法、三解分解法程序，分析運行結果。3 調試運行列主元高斯消去法、列主元三解分解法算法程序。4 并用上述幾種算法程序計算出上面兩個方程組的解。六、實驗步驟1. 根據實驗題目，給出解決問題的程序代碼。2. 上機輸入和調試自己所編的程序。3. 上機結束后，應整理出實驗報告。七、實驗報告要求及記錄、格式按金陵科技學院實驗報告（工科）格式填寫 附1 :列主元高斯消去法源程序:*/*/*列主元高斯消去法求線性方程組的解/*/#i nclude &

10、lt;stdio.h>#in elude <math.h>#define Max_N10/*方程組最大維數 */*列主元高斯消去法函數*/void ColPivot(float AMax_NMax_N,float B,int n) int i,j,k,m_i;float m_x,temp;for(i=0;i< n-1;i+) /*列主元*/j=i+1; m=i; m_x=fabs(Aii);for(;j< n;j+)if(fabs(Aji>m_x)/* 找主元素 */m_i=j； m_x=fabs(Aji);if(i<m)/*交換兩行*/ temp=B

11、i; Bi=Bm_i; Bm_i=temp; for(j=i;j< n;j+) temp=Aij; Aij=Amj; Amj=temp; /*消元*/for(j=i+1;j<n;j+) temp=-Aji/Aii;Bj+=Bi*temp; for(k=i;k<n;k+)Ajk+=Aik*temp;main() /* 主函數 */ int i,j,k,n;float aMax_NMax_N,bMax_N,xMax_N;printf("nPlease input n value(dim of Ax=b):"); /*輸入矩陣維數 */do scanf(&quo

12、t;%d",&n); if(n>Max_N) printf("nplease re-input n value:"); while(n>Max_N|n<=0); /* 輸入 Ax=b 的 A 矩陣 */ printf("Input the A(i,j):n"); for(i=0;i<n;i+) for(j=0;j<n;j+) scanf("%f",&aij);/*輸入 b 矩陣*/printf("Input b(i):n"); for(i=0;i<n;i

13、+) scanf("%f",&bi);ColPivot(a,b,n); /* 調用列主元消去法函數計算方程組的解 */ xn-1=bn-1/an-1n-1;/* 解方程 */for(i=n-2;i>=0;i-) xi=bi;for(j=i+1;j<n;j+) xi-=aij*xj;xi/=aii;printf("Solve is :"); /* 輸出方程組的解 */ for(i=0;i<n;i+)printf("x%d=%f ",i,xi); if(i%2=0)printf("n");/*

14、End of file*/*Please input n value(dim of Ax=b):3Input the A(i,j):2 1 11 3 21 2 2Input b(i):4 6 5Solve is:x0=1.000000 x1=1.000000 x2=1.000000 */附 2 ：三角分解法源程序： /*/* 直接三角分解法 (LU 分解法 )求線性方程組的解 */ /*/#include <stdio.h> #include <math.h>#define Max_N 10/* 最大維數 */* 直接三角分解法函數 */float *DirectLU(

15、float aMax_NMax_N,float b,int n)int i,j,k;float yMax_N,LMax_NMax_N,UMax_NMax_N,xMax_N;/*U 矩陣對角元素賦值為 1*/ for(i=0;i<n;i+) Uii=1; for(k=0;k<n;k+) for(i=k;i<n;i+) /* 計算 L 矩陣的第 k 列元素 */ Lik=aik; for(j=0;j<=k-1;j+) Lik-=(Lij*Ujk);for(j=k+1;j<n;j+) /* 計算 U 矩陣的第 k 行元素 */ Ukj=akj; for(i=0;i<

16、;=k-1;i+) Ukj-=(Lki*Uij);Ukj/=Lkk;for(i=0;i<n;i+) /* 計算 Ly=b 中的 y*/ yi=bi;for(j=0;j<=i-1;j+)yi-=(Lij*yj);yi/=Lii;for(i=n-1;i>=0;i-)/* 計算 Ux=y 中的 x*/xi=yi;for(j=i+1;j<n;j+)xi-=(Uij*xj);return(x);main() /* 主函數 */ int i,j,k,n;float temp;float aMax_NMax_N,bMax_N,*x;printf("nPlease input

17、 n value(dim of Ax=b):"); /*輸入矩陣維數 */do scanf("%d",&n); if(n>Max_N) printf("nplease re-input n value:");while(n>Max_N|n<=0);/* 輸入 Ax=b 的 A 矩陣 */ printf("Input the A(i,j):n");for(i=0;i<n;i+) for(j=0;j<n;j+)scanf("%f",&aij);/*輸入 b 矩陣*

18、/ printf("Input b(i):n");for(i=0;i<n;i+) scanf("%f",&bi);x=DirectLU(a,b,n); /* 調用直接三角分解法函數 */ printf("Solve is :");/* 輸出方程組的解 */for(i=0;i<n;i+) printf("x%d=%f",i,xi);if(i%2=0) printf("n");/* End of file */* 程序輸入輸出：Please input n value(dim o

19、f Ax=b):3Input the A(i,j):2 1 11 3 21 2 2Input b(i):4 6 5Solve is:x0=1.000000 x1=1.000000 x2=1.000000 */實驗3插值法、實驗目的1 掌握插值函數的概念，插值多項式的唯一性。2 掌握插值余項，差分及等距插值公式，高次插值的誤差分析。3 .掌握基本插值多項式，拉格朗日插值多項式，差商，牛頓插值多項式。、實驗要求1上機前作好充分準備，比較不用的方法解決相同問題的不同。2 上機時要遵守實驗室的規(guī)章制度，愛護實驗設備。3 記錄調試過程及結果，記錄并比較與手工運算結果的異 同。4 程序調試完后，須由實驗輔

20、導教師在機器上檢查運行結 果。5 給出本章實驗單元的實驗報告。三、實驗設備、環(huán)境1 .硬件設備：IBM PC以上計算機，有硬盤和一個軟驅、單 機和網絡環(huán)境均可。2 .軟件環(huán)境： C語言運行環(huán)境。四、實驗原理、方法1、拉格朗日插值算法步驟:(1 )、輸入 n, Xi , y (i=0, 1, 2,，n),給初值 Ln(x) =0(2)、對 i=0, 1, 2，n 計算nli(x)可丨j=o j wX _XjXi _XjLn(x) =Ln(x) Ti(x)yi(3)、輸出 Ln(x)2、牛頓插值法算法步驟(1 )、輸入 n, x'(要求其函數值),Xi , y (i=0,1,2,n)；(2

21、 )、對 k=1,2,3,，n ,i=1,2,，k計算各階均差拉格朗日插值框圖fX0,X1, ,Xk；(3)、利用下面的牛頓插值公式計算x'的函數值Nn(x) =f(X。) fX0,X1(X -X0) fX0,X1,X2(X -xj(x -X。)十+ f X。，知,Xn(X Xo)(X Xj(XX)(4)、輸出函數值。五、實驗內容X1.6151.6341.7021.8281.921Y=f(x)2.414502.464592.652713.030353.34066仁已知函數表(1 )用拉格朗日插值(二次、四次)計算 f(1.682)和f(1.813) 的近似值。(2)構造出均差表，并利用

22、牛頓(均差)插值多項式計算 f(1.682)和 f(1.813)的近似值。(3 )分析并比較兩種算法得到的近似值的精度。2 編寫拉格朗日和牛頓插值算法程序，分析運行結果。六、實驗步驟1. 根據實驗題目，給出解決問題的程序代碼。2. 上機輸入和調試自己所編的程序。3. 上機結束后，應整理出實驗報告。七、實驗報告要求及記錄、格式按金陵科技學院實驗報告(工科)格式填寫附1 ：拉格朗日插值法核心程序：for(i=0;i<N;i+)/*計算拉格朗日插值函數的值*/ fi=yi;for(j=0;j<N;j+)if(j!=i) fi*=(XX-Xj)/(Xi-Xj); yy+=fi;附2 :牛頓

23、插值法核心程序：for(l=1;l<=N;l+)/*計算牛頓插值函數的值*/ f0=yI;for(j=0;j<I;j+)/* 計算均差*/fj+1=(fj-yj)/(xI-xj);yI=fI;b=yN;牛頓插值法框圖for(I=N-1;I>=0;I-) b=b*(xx-xI)+yI; /* 計算函數值 */實驗4曲線擬合、實驗目的1 掌握最小二乘原理，正規(guī)方程組，超定方程組概念。2 掌握用最小二乘法擬合曲線，超定方程級的最小二乘解。3 掌握用最小二乘法擬合曲線。、實驗要求1 上機前作好充分準備，復習最小二乘擬合方法。2 上機時要遵守實驗室的規(guī)章制度，愛護實驗設備。3 記錄調試

24、過程及結果，記錄并比較與手工運算結果的異同。4 程序調試完后，須由實驗輔導教師在機器上檢查運行結果。5 給出本章實驗單元的實驗報告。、實驗設備、環(huán)境i .硬件設備：IBM PC以上計算機，有硬盤和一個軟驅、單機和網絡環(huán)境均可。2 .軟件環(huán)境：C語言運行環(huán)境。四、實驗原理、方法最小二乘法：對于給定的線性方程組Ax=b式中A=aii ai2 a2i a 22ain Ia2nb=b2_ami am2 a mn_bmX2x=_Xn當m>n時，稱為矛盾方程組，又稱超定方程組。對于這種方程組有m個方程，而只有比m小的n個變量，即方程的個數超過未知量的個數，這種方程組一般來說是沒有解的。我們 轉而尋求

25、在某種意義下的近似解。如果這組近似解對于矛盾方程組中的每個方程式的誤差的 平方和為最小，即mm n（Xi,X2，Xn）八汀 c ajXj - bi）2imi# j=!為最小值時，就可以認為該組近似值為矛盾方程組的近似解。這種近似解不是指對精確解的近似（因為精確解并不存在），而是指尋求各未知數的一組取值，使方程式（5.i）中各方程式近似相等，這就是最小二乘法的基本思想。五、實驗內容i 設有如下實驗數據x1.361.491.731.811.952.162.282.48y14.09415.06916.84417.37818.43519.94920.96322.495試用最小二乘法分別求一次及二次多項

26、式曲線擬合以上數據。2 編寫程序，分析運行結果。六、實驗步驟1. 根據實驗題目，給出解決問題的程序代碼。2. 上機輸入和調試自己所編的程序。3. 上機結束后，應整理出實驗報告。七、實驗報告要求及記錄、格式按金陵科技學院實驗報告（工科）格式填寫附：拋物函數擬合源程序/*/*/*最小二乘法-擬合拋物函數/*/#i nclude <stdio.h>#in elude <math.h>#define Max_N 25/*最大數據點的個數 */#define M 3/*正規(guī)方程組的階數*/*列主元高斯消去法求解線性方程組*/void ColPivot(float AMM,floa

27、t B,i nt n) int i,j,k,m_i;float m_x,temp;for(i=0;i< n-1;i+) /*列主元*/j=i+1; m_i=i; m_x=fabs(Aii);for(;j< n;j+) if(fabs(Aji>m_x) /* 找主元素 */m_i=j;m_x=fabs(Aji);if(i<m)/*交換兩行*/ temp=Bi; Bi=Bm_i; Bm_i=temp; for(j=i;j< n;j+) temp=Aij; Aij=Amj; Am口 j=temp; /*消元*/for(j=i+1;j<M;j+) temp=-Aji

28、/Aii;Bj+=Bi*temp; for(k=i;k<M;k+)Ajk+=Aik*temp;main()int i,j,k,n;float xMax_N,yMax_N,bM,aMM,cM; printf("nPlease input n value:"); /* 輸入點數 n*/ do scanf("%d",&n);if(n>Max_N) printf("nplease re-input n value:"); while(n>Max_N|n<=0); printf("Input xi,i=

29、0,.%d:n",n-1); for(i=0;i<n;i+) scanf("%f",&xi); printf("Input yi,i=0,.%d:n",n-1); for(i=0;i<n;i+) scanf("%f",&yi); for(i=0;i<M;i+) /* 構造正規(guī)方程組 */ for(j=0;j<M;j+)aij=0; bi=0; for(k=0;k<n;k+) aij=aij+pow(xk,i+j); bi=bi+pow(xk,i)*yk;/* 輸出正規(guī)方程組 fo

30、r(i=0;i<M;i+) for(j=0;j<M;j+) printf("%f",aij);printf(" %f",bi); printf("n");*/ColPivot(a,b,M); /* 調用列主元消去法函數計算方程組的解 */cM-1=bM-1/aM-1M-1; /* 解方程 */ for(i=M-2;i>=0;i-) ci=bi;for(j=i+1;j<M;j+)ci-=aij*cj;ci/=aii;printf("Solve is :n"); /* 輸出方程組的解 */for

31、(i=0;i<M;i+) printf("c%d=%fn",i,ci);printf("Result ： y=%f+(%f)x+(%f)x 2",c0,c1,c2); getch();/* End of file */* 程序輸入輸出Please input n value:16Input xi,i=0,.N-1:1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16Input yi,i=0,.N-1:4 6.4 8 8.8 9.22 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10

32、.58 10.60 Solve is : c0=4.387500 c1=1.065962 c2=-0.044466Result:y=4.387500+1.065962x+(-0.044466)x 2*/實驗5數值積分與數值微分、實驗目的1.掌握插值型求積公式，梯形公式、辛卜生公式、柯特斯公式。2 .掌握高斯求積公式，數值微分的中點公式。3 掌握龍貝格求積公式。二、實驗要求1 .上機前作好充分準備，理解求積公式，梯形公 式的概念。2 上機時要遵守實驗室的規(guī)章制度，愛護實驗設 備。3 記錄調試過程及結果，記錄并比較與手工運算 結果的異同。4 程序調試完后，須由實驗輔導教師在機器上檢 查運行結果。5

33、給出本章實驗單元的實驗報告。三、實驗設備、環(huán)境1 .硬件設備：IBM PC以上計算機，有硬盤和一 個軟驅、單機和網絡環(huán)境均可。2 .軟件環(huán)境：C語言運行環(huán)境。四、實驗原理、方法a) 梯形求積公式：復化梯形公式計算步驟、令h=(b-a)/N , T=0 (h為等分數，T為存 放積分值的變量)(2) 、對k=1,2,N計算T=T+f(a+kh)(3) 、T=hf(a)+2T+f(b)/2b) 辛卜生求積公式 復化辛卜生公式計算步驟(1) 、令 h=(b-a)/N , s仁f(a+h/2), s2=0(2) 、對 k=1,2,N-1 計算s仁s1+f(a+kh+h/2), s2=s2+f(a+kh)

34、(3) 、s=hf(a)+4s1+2s2+f(b)/2。復化梯形公式流程圖定義f(x)1/輸入 a,b,N /開始h=(b-a)/N,x=a+h/2,s1=f(x),s2=0s1=s1+f(a+kh+h/2)s2=s2+f(a+kh)(k=1,2,N-1)c)龍貝格求積公式算法流程見：龍貝格求積公式流程圖s2=hf(a)+4s1+2s2+f(b)/6結束復化辛卜生公式流程圖開始五、實驗內容1 .分別編寫出梯形公式、辛卜生求積公式、龍 貝格求積公式算法程序。輸入a,2 分別用上述算法程序求I01 xh=b-a, k=1T1=hf(a)+f(b)/2值，并分析每種算法程序求得的結果為什么不 同。六

35、、實驗步驟根據實驗題目，給出解決問題的程序代碼。 上機輸入和調試自己所編的程序。上機結束后，應整理出實驗報告。1.2.3.七、實驗報告要求及記錄、格式按金陵科技學院實驗報告（工科）格式填寫附：龍貝格求積分法源程序/*/*用龍貝格求積公式求定積分的解/*/*/#in elude <stdio.h>#in clude <math.h>#define epsilon0.0001 /* 精度要求,高，否則不能得岀結果*/*求積函數f(x)*/float f(float x)retur n(si n( x); /*梯形公式*/S=0，x=a+h/2S=S+f(x)， x=x+h是

36、x<b?T2=(T1+h*S)/2S2=T2+(T2-T1)/3是k=1?S仁S2T1=T2h=h/2C2=S2+(S2-S1)/15d 1精度不能太JL是-+6= C2k=k+1float Romberg(float afloat b)int k=1;float S,x,T1,T2,S1,S2,C1,C2,R1,R2,h=b-a;T1=h*(f(a)+f(b)/2;while(1)S=0; x=a+h/2;doR2=C2+(C2-C1)/63是k=3?R仁R2|R2-R1|> £ ?S+=f(x);x+=h;while(x<b);T2=(T 1+h*S)/2.0;

37、if(fabs (T2-T1)<epsilo n) return (T2);S2=T2+(T2-T1)/3.0;輸出R2結束龍貝格求積公式流程圖if(k=1) T1=T2;S1=S2;h/=2;k+=1; continue; C2=S2+(S2-S1)/15.0;if(k=2) C1=C2;T1=T2;S1=S2;h/=2;k+=1; continue; R2=C2+(C2-C1)/63.0;if(k=3) R1=R2;C1=C2;T1=T2;S1=S2;h/=2;k+=1; continue; if(fabs(S2-S1)<epsilon) return(S2);R1=R2;C1

38、=C2;T1=T2;S1=S2;h/=2;k+=1;main() int I;float a,b,S;printf( “nInput the begin and end: ”);/* 輸入積分區(qū)間 */scanf( “ %f%f” ,&a,&b);S=Romberg(a,b); /* 調用龍貝格算法函數 */ printf( “ Solve is: %f” ,S);getch();/* End of file */*對函數 f(x)=sin(x) 在積分區(qū)間 1,2 內求定積分。 程序輸入輸出Input the begin and end: 1.0 2.0Solve is: 0

39、.956449*/充分理解各龍格-庫塔法流程圖改進歐拉法流程圖實習6常微分方程數值解法、實驗目的1、通過本實驗，充分理解常微分方程的初值問題的有關方法和理論理論；2、通過實際計算體會各種解法的功能、優(yōu)缺點及適用場合，會選取適當的求解方法。、實驗要求1、了解各種解決常微分方程初值問題的方法，并比較各種方法的不同之處, 種方法的特點及用途；2、針對后面的練習題目進行上機計算；3、在充分理解各種方法的基礎上，會編寫各種方法的程序；4、考慮其他龍格-庫塔公式的程序編寫。三、實驗設備、環(huán)境1 .硬件設備：IBM PC以上計算機，有硬盤和 個軟驅、單機和網絡環(huán)境均可。2 .軟件環(huán)境：C語言運行環(huán)境。四、實驗原理、方法1、改進歐拉公式算法框圖見右側流程圖2、龍格-庫塔公式（標準四階龍格 -庫塔公式）五、實驗內容1、用歐拉方式和改進歐拉法求初值問題y'=x*y; y（3.0）=1; 3.0<=x<=3.6 取 h=0.2（ 即 n=3）2、龍格 -庫塔公式求初值問題