數(shù)值計算課程設(shè)計

1、1、 經(jīng)典四階龍格庫塔法解一階微分方程1.1、 算法說明龍格-庫塔(Runge-Kutta)方法是一種在工程上應(yīng)用廣泛的高精度單步算法。由于此算法精度高，采取措施對誤差進行抑制，所以其實現(xiàn)原理也較復(fù)雜。該算法是構(gòu)建在數(shù)學(xué)支持的基礎(chǔ)之上的。4階龍格-庫塔方法(RK4)可模擬N=4的泰勒方法的精度。這種算法可以描述為，自初始點開始，利用下面的計算方法生成近似序列(1-1)1.2、 經(jīng)典四階龍格庫塔法解一階微分方程算法流程圖圖1-1 經(jīng)典四階龍格庫塔法解一階微分方程算法流程圖1.3、 經(jīng)典四階龍格庫塔法解一階微分方程程序調(diào)試圖1-2 經(jīng)典四階龍格庫塔法解一階微分方程程序調(diào)試1.4、 經(jīng)典四階龍格庫塔

2、法解一階微分方程代碼#include #include using namespace std;/f為函數(shù)的入口地址，x0、y0為初值，xn為所求點，step為計算次數(shù)double Runge_Kuta( double (*f)(double x, double y), double x0, double y0, double xn, int step ) double k1,k2,k3,k4,result; double h=(xn-x0)/step; if(step=0) return(y0); if(step=1) k1=f(x0,y0); k2=f(x0+h/2, y0+h*k1/2)

3、; k3=f(x0+h/2, y0+h*k2/2); k4=f(x0+h, y0+h*k3); result=y0+h*(k1+2*k2+2*k3+k4)/6; else double x1,y1; x1=xn-h; y1=Runge_Kuta(f, x0, y0, xn-h,step-1); k1=f(x1,y1); k2=f(x1+h/2, y1+h*k1/2); k3=f(x1+h/2, y1+h*k2/2); k4=f(x1+h, y1+h*k3); result=y1+h*(k1+2*k2+2*k3+k4)/6; return(result);int main()double f(d

4、ouble x, double y); double x0,y0; double a,b;/ int step; coutx0y0; coutab; /double x0=0,y0=1; double x,y,step; int i; coutstep; /step=0.1; cout.precision(10); for(i=0;i=(b-a)/step;i+)x=x0+i*step; coutsetw(8)xsetw(18)Runge_Kuta(f,x0,y0,x,i)endl;return(0); double f(double x, double y) double r; r=(x-y

5、)/2; return(r);2、 高斯列主元法解線性方程組2.1、 算法說明首先將線性方程組做成增光矩陣，對增廣矩陣進行行變換。對第元素，在第i列中，第i行及以下的元素選取絕對值最大的元素，將該元素最大的行與第i行交換，然后采用高斯消元法將新得到的消去第i行以下的元素。一次進行直到。從而得到上三角矩陣。再對得到的上三角矩陣進行回代操作，即可以得到方程組的解。2.2、 高斯列主元算法流程圖i=0;j=0NY開始輸入未知數(shù)個數(shù)mim*(m+1)max j=i:max=|*(head+(m+1)*i+i)|;kmaxmax =|*(head+(m+1)*k+i)|;max i=k;NYYKmK=k

6、+1;YNmax i!=itemp=*(head*(m+1)*i+k);*(head*(m+1)*i+k)=*(head*(m+1)*max i+k)*(head*(m+1)*max i+k)=tempk=0Nk=m+1YNYi!=jk=0*(head*(m+1)*i+k)=*(head+(m+1)*j+k)*(head*(m+1)*j+i) *(head*(m+1)* i+k)/(*head+(m+1)*i+i)k=m+1YNj=j+1YjmYNi=0NY結(jié)束*(head+(m+1)*i+m)=*(head+(m+1)*i+m)/(*head+(m+1)*i+m)；i=i+1；圖2-1 算法

7、流程圖輸出*(head+(m+1)*i+m)jm2.3、 高斯列主元程序調(diào)試對所編寫的高斯列主元程序進行編譯和鏈接，然后執(zhí)行得如下所示的窗口，我們按命令輸入增廣矩陣的行數(shù)為3，輸入3行5列的增廣矩陣,運行界面為：圖2-2 高斯列主元程序調(diào)試2.4、 高斯列主元算法代碼#include#includeusing namespace std;void load();const N=20;float aNN;int m;int main()int i,j;int c,k,n,p,r;float xN,lNN,s,d;coutm;coutendl;cout請按順序輸入增廣矩陣a:endl;load()

8、;for(i=0;im;i+) for(j=i;jfabs(aii)?j:i; /*找列最大元素*/for(n=0;nm+1;n+) s=ain; ain=acn; acn=s; /*將列最大數(shù)防在對角線上*/for(p=0;pm+1;p+)coutaipt;coutendl;for(k=i+1;km;k+) lki=aki/aii; for(r=i;r=0;i-) d=0;for(j=i+1;jm;j+)d=d+aij*xj;xi=(aim-d)/aii;cout該方程組的解為:endl;for(i=0;im;i+)coutxi=xit; /system(pause); return 0;

9、void load()int i,j;for(i=0;im;i+)for(j=0;jaij;3、 牛頓法解非線性方程組3.1、 算法說明設(shè)已知。第1步：計算函數(shù) （3-1）第2步：計算雅可比矩陣（3-2）第3步：求線性方程組 的解。第4步：計算下一點 （3-3）重復(fù)上述過程。3.2、 牛頓法解非線性方程組算法流程圖圖3-1 算法流程圖3.3、 牛頓法解非線性方程組算法程序調(diào)試圖3-2牛頓法解非線性方程組算法程序調(diào)試應(yīng)用本程序解方程組， 初始近似值x0,y0分別為2.00和0.25，經(jīng)過3次迭代求出X(1)=1.900691和X(2)=0.311213。圖3-2牛頓法解非線性方程組算法程序運行結(jié)

10、果3.4、 牛頓法解非線性方程組算法程序代碼#include#include#define N 2 / 非線性方程組中方程個數(shù)、未知量個數(shù) #define epsilon 0.0001 / 差向量1范數(shù)的上限#define max 10 /最大迭代次數(shù)using namespace std;const int N2=2*N;int main()void ff(float xxN,float yyN);void ffjacobian(float xxN,float yyNN);void inv_jacobian(float yyNN,float invNN);void newdim(float

11、x0N, float invNN,float y0N,float x1N);float x0N=2.0,0.25,y0N,jacobianNN,invjacobianNN,x1N,errornorm;int i,iter=0;cout初始解向量：endl;for (i=0;iN;i+) coutx0i ; coutendl;do iter=iter+1;cout第 iter 次迭代開始：endl;/jis uan xiang liang han shu zhi-yin bian liang xiang liang y0ff(x0,y0);/jis uan jacobian ju zhen ja

12、cobianffjacobian(x0,jacobian);/jis uan jacobian ju zhen de ni juzhen invjacobianinv_jacobian(jacobian,invjacobian);/you jie xiang liang x0 ji suan jie xiang liang x1newdim(x0, invjacobian,y0,x1);/ji suan cha xiang liang de 1 fan shuerrornorm=0;for (i=0;iN;i+) errornorm=errornorm+fabs(x1i-x0i);if (er

13、rornormepsilon) break;for (i=0;iN;i+)x0i=x1i; while (itermax);return 0;void ff(float xxN,float yyN) float x,y; int i; x=xx0;y=xx1;/非線性方程組yy0=x*x-2*x-y+0.5;yy1=x*x+4*y*y-4; cout因變量向量：endl;for( i=0;iN;i+) coutyyi ; coutendl; coutendl;void ffjacobian(float xxN,float yyNN)float x,y;int i,j;x=xx0;y=xx1;y

14、y00=2*x-2;yy01=-1;yy10=2*x;yy11=8*y;cout雅克比矩陣： endl;for( i=0;iN;i+)for(j=0;jN;j+) coutyyij ; coutendl;coutendl;void inv_jacobian(float yyNN,float invNN)float augNN2,L; int i,j,k; cout計算雅克比矩陣的逆: endl; for (i=0;iN;i+)for(j=0;jN;j+) augij=yyij; for(j=N;jN2;j+) if(j=i+N) augij=1; else augij=0; for (i=0;

15、iN;i+)for(j=0;jN2;j+) coutaugij ; coutendl;coutendl;for (i=0;iN;i+) for (k=i+1;kN;k+)L=-augki/augii; for(j=i;jN2;j+) augkj=augkj+L*augij;for (i=0;iN;i+)for(j=0;jN2;j+) coutaugij ; coutendl;cout0;i-) for (k=i-1;k=0;k-)L=-augki /augii; for(j=N2-1;j=0;j-) augkj=augkj+L*augij;for (i=0;iN;i+)for(j=0;jN2;

16、j+) coutaugij ; coutendl;cout=0;i-)for(j=N2-1;j=0;j-)augij=augij/augii;for (i=0;iN;i+)for(j=0;jN2;j+) coutaugij ; coutendl; for(j=N;jN2;j+) invij-N=augij;coutendl;cout雅克比矩陣的逆： endl;for (i=0;iN;i+) for(j=0;jN;j+) coutinvij ; coutendl; coutendl;void newdim(float x0N, float invNN,float y0N,float x1N)in

17、t i,j;/計算非線性方程組的近似解向量float sum=0;for(i=0;iN;i+) sum=0; for(j=0;jN;j+) sum=sum+invij*y0j;x1i=x0i-sum; cout近似解向量：endl;for (i=0;iN;i+) coutx1i endl; 4、 龍貝格求積分算法4.1、 算法說明生成的逼近表，并以為最終解來逼近積分（4-1）逼近存在于一個特別的下三角矩陣中，第0列元素用基于個a,b子區(qū)間的連續(xù)梯形方法計算，然后利用龍貝格公式計算。當(dāng)時，第行的元素為 （4-2）當(dāng)時，程序在第行結(jié)束。4.2、 龍貝格求積分算法流程圖圖4-1 算法流程圖圖4-1

18、算法流程圖4.3、 龍貝格求積分算法程序調(diào)試我們以求解積分，精度為0.0001，最高迭代10次為例，對所編寫的龍貝格求積分算法程序進行編譯和鏈接，經(jīng)執(zhí)行后得如下所示的窗口圖4-2龍貝格求積分算法程序調(diào)試說明：應(yīng)用Romberg算法求在區(qū)間上的精度為0.0001的積分為0.494508。4.4、 龍貝格求積分算法代碼#include#includeusing namespace std;#define f(x) sin(x*x) /舉例函數(shù)#define epsilon 0.0001 /精度#define MAXREPT 10 /迭代次數(shù),到最后仍達不到精度要求,則輸出T(m=10).doubl

19、e Romberg(double aa, double bb) /aa,bb 積分上下限 int m, n;/m控制迭代次數(shù), 而n控制復(fù)化梯形積分的分點數(shù). n=2m double h, x; double s, q; double ep; /精度要求 double *y = new doubleMAXREPT;/為節(jié)省空間,只需一維數(shù)組 /每次循環(huán)依次存儲Romberg計算表的每行元素,以供計算下一行,算完后更新double p;/p總是指示待計算元素的前一個元素(同一行) /迭代初值 h = bb - aa; y0 = h*(f(aa) + f(bb)/2.0; m = 1; n = 1

20、; ep = epsilon + 1.0; /迭代計算 while (ep = epsilon) & (m MAXREPT) /復(fù)化積分公式求T2n(Romberg計算表中的第一列),n初始為1,以后倍增 p = 0.0; for (int i=0; in; i+)/求Hn x = aa + (i+0.5)*h; p = p + f(x); p = (y0 + h*p)/2.0;/求T2n = 1/2(Tn+Hn),用p指示 /求第m行元素,根據(jù)Romberg計算表本行的前一個元素(p指示), /和上一行左上角元素(yk-1指示)求得. s = 1.0; for (int k=1; k=m;

21、k+) s = 4.0*s; q = (s*p - yk-1)/(s - 1.0); yk-1 = p; p = q; p = fabs(q - ym-1); m = m + 1; ym-1 = q; n = n + n; h = h/2.0; return (q);int main() double a,b; coutRomberg積分,請輸入積分范圍a,b:ab; cout積分結(jié)果:Romberg(a, b)endl; return 0;5、 三次樣條插值算法5.1 三次樣條插值算法說明表表5-1三次樣條插值算法說明表策略描述包含和的方程(i)三次緊壓樣條，確定，（如果導(dǎo)數(shù)已知，這是“最佳

22、選擇”）(ii)natural三次樣條（一條“松弛曲線”），(iii)外掛到端點(iv)是靠近端點的常量，(v)在每個端點處指定，5.2、 三次樣條插值算法（壓緊樣條）程序調(diào)試我們將所編寫的程序三次樣條插值算法（壓緊樣條）程序進行調(diào)試圖5-1三次樣條插值算法（壓緊樣條）程序輸入界面、運行結(jié)果圖5-2三次樣條插值算法程序運行結(jié)果（a） 圖5-2三次樣條插值算法程序運行結(jié)果 （b）運行結(jié)果分析： （5-1）作圖程序（Matlab）：x1=0:0.01:1;y1=0.06*(x1 - 1).3 + 0.42*(x1 - 0).3 - 0.06*(x1 - 1) + 0.08*(x1 - 0);x2=

23、1:0.01:2;y2=-0.42*(x2 - 2).3 - 0.62*(x2 - 1).3 - 0.08*(x2 - 2) + 2.62*(x2 - 1);x3=2:0.01:3;y3=0.62*(x3 - 3).3 + 0.06*(x3 - 2).3 - 2.62*(x3 - 3) + 1.44*(x3 - 2);X=0 1 2 3;Y=0 0.5 2 1.5;plot(x1,y1,x2,y2,x3,y3,X,Y,*)gtext(S1)gtext(S2)gtext(S3)圖形為：圖5-3 三次樣條插值算法Matlab作圖分析5.3、 三次樣條插值算法（壓緊樣條）代碼#include#inc

24、ludeusing namespace std;const int max = 50;float xmax, ymax, hmax;float cmax, amax, fxymmax;float f(int x1, int x2, int x3) float a = (yx3 - yx2) / (xx3 - xx2); float b = (yx2 - yx1) / (xx2 - xx1); return (a - b)/(xx3 - xx1); /求差分void cal_m(int n) /用追趕法求解出彎矩向量M float Bmax; B0 = c0 / 2; for(int i = 1

25、; i n; i+) Bi = ci / (2 - ai*Bi-1); fxym0 = fxym0 / 2; for(i = 1; i = 0; i-) fxymi = fxymi - Bi*fxymi+1;void printout(int n);int main() int n,i; char ch; do coutn; for(i = 0; i = n; i+) coutPlease put in Xixi; coutPlease put in Yiyi; for(i = 0; i n; i+) /求步長 hi = xi+1 - xi; coutt; switch(t) case 1:c

26、out輸入 Y0 Ynf0f1; c0 = 1; an = 1; fxym0 = 6*(y1 - y0) / (x1 - x0) - f0) / h0; fxymn = 6*(f1 - (yn - yn-1) / (xn - xn-1) / hn-1; break; case 2:cout輸入 Y0 Ynf0f1; c0 = an = 0; fxym0 = 2*f0; fxymn = 2*f1; break; default:cout不可用n;/待定 ; for(i = 1; i n; i+) fxymi = 6 * f(i-1, i, i+1); for(i = 1; i n; i+) ai

27、 = hi-1 / (hi + hi-1); ci = 1 - ai; an = hn-1 / (hn-1 + hn); cal_m(n); coutn輸出三次樣條插值函數(shù)：n; printout(n); coutch; while(ch = y | ch = Y); return 0;void printout(int n)coutsetprecision(6); for(int i = 0; i n; i+)couti+1: xi , xi+1nt;coutSi+1 0) cout-t*(x - xi+1)3; else cout-t*(x - xi+1 0) cout + t*(x -

28、xi)3; else cout - t*(x - xi)3; cout 0) cout- t*(x - xi+1); else cout- -t*(x - xi+1 0) cout + t*(x - xi); else cout - -t*(x - xi); coutendl; coutendl;6、M次多項式曲線擬合6.1、算法說明設(shè)有N個點，橫坐標是確定的。最小二乘拋物線的系數(shù)表示為（6-1）求解A，B和C的線性方程組為 （6-2）6.2、 M次多項式曲線擬合算法流程圖圖6-1 算法流程圖6.3、 M次多項式曲線擬合算法程序調(diào)試我們按命令依次輸入命令如下命令后，得程序執(zhí)行結(jié)果如下圖6-2

29、算法調(diào)試圖作圖程序： 圖形為：X=-3 -1 1 3;Y=15 5 1 5;x=-3:0.01:3;y=0.875*x.2-1.7*x+2.145;plot(X,Y,go,x,y,r-)gtext(擬合曲線)圖6-3 圖形分析6.4、 M次多項式曲線擬合算法代碼#include#include#define MAX 20using namespace std;/求解任意可逆矩陣的逆，X為待求解矩陣，E為全零矩陣，非單位矩陣，也可以是單位矩陣void inv(double XMAXMAX,int n,double EMAXMAX)int i,j,k;double temp=0;for(i=0;i

30、MAX;i+)for(j=0;jMAX;j+)if(i=j)Eij=1;for(i=0;in-1;i+)temp=Xii;for(j=0;jn;j+)Xij=Xij/temp;Eij=Eij/temp;for(k=0;kn;k+)if(k=i)continue;temp=-Xii*Xki;for(j=0;jn;j+)Xkj=Xkj+temp*Xij;Ekj=Ekj+temp*Eij;int main()int n,M,i,j,k;double XMAX=0,YMAX=0,FMAXMAX=0,BMAX=0;double AMAXMAX=0,BFMAXMAX=0,CMAX=0,EMAXMAX=0;

31、coutn;coutn請輸入n個點的X坐標序列:n;for(i=0;iXi;coutn請輸入n個點的Y坐標序列:n;for(i=0;iYi;coutM;for(i=0;in;i+)for(k=1;k=M+1;k+)Fik-1=pow(Xi,k-1);/求F的轉(zhuǎn)置for(i=0;in;i+)for(j=0;jM+1;j+)BFji=Fij; /計算其轉(zhuǎn)置的BF與F的乘 for(i=0;iM+1;i+)for(j=0;jM+1;j+)for(k=0;kn;k+)Aij+=BFik*Fkj;/計算F的轉(zhuǎn)置BF與Y的乘for(i=0;iM+1;i+)for(j=0;jn;j+)Bi+=BFij*Yj;

32、/調(diào)用inv函數(shù)求解矩陣A的逆矩陣Einv(A,n,E);/計算A的逆BF與B的乘for(i=0;iM+1;i+)for(j=0;jn;j+)Ci+=Eij*Bj;coutn擬合后的M=0;i-)coutCit;coutn擬合后的M次多項式為:n;cout=0;i-)if(i=0)cout+Ci;elsecout+Ci*xi;coutendl;return 0;7、 牛頓-拉弗森迭代解非線性方程7.1、 牛頓-拉弗森迭代解非線性方程算法概要使用初始近似值，利用迭代 其中計算函數(shù)的根的近似值。7.2、 牛頓-拉弗森迭代解非線性方程算法程序調(diào)試圖7-1牛頓-拉弗森迭代解非線性方程算法程序調(diào)試7.3

33、、 牛頓-拉弗森迭代解非線性方程算法代碼#includeusing namespace std;#includedouble f(double x)return x*x-2*x-1double df(double x)return 2*x-2;int main() int k=0,max1=0; double p0,delta,epsilon,p1,err,relerr,y;cout牛頓拉弗森法解非線性方程f(x)=x2-2x-1endl;cout初始值為p0=4.8endl; p0=4.8,delta=0.0001,epsilon=0.0001,max1=100;for(k=1;k=max1

34、;k+)p1=p0-f(p0)/df(p0);err=fabs(p1-p0);relerr=2*err/(fabs(p1)+delta);p0=p1;y=f(p0);if(errdelta)|(relerrdelta)|(fabs(y)epsilon)break;cout方程的根的近似值為:p0endl;cout方程的根的誤差估計為:errendl;cout迭代次數(shù)為:kendl;cout方程在p0點的函數(shù)值f(p0):yendl;return 0;8、 不動點法解非線性方程8.1、 算法說明先將改寫成然后對進行迭代，即 其中然后判斷是否成立，成立則返回，不成立就重復(fù)以上步驟8.2、 不動點法

35、解非線性方程算法程序調(diào)試我們將編寫好的不動點迭代法解非線性方程算法程序進行編譯，鏈接和執(zhí)行后得如下所示結(jié)果。圖8-1 程序調(diào)試圖8.3、 不動點法解非線性方程算法代碼#include#include#include#define MAX 20#define eps 1e-10using namespace std;double g(double x)return 1-x*x/4+x; main()double PMAX=0,err=0.0,relerr=0.0,tol=0.0,p=0.0,p0=0.0;int k=0,max1=0,i=0;cout不動點法解非線性方程f(x)=1-x2/2en

36、dl;cout方程在0，1上有解，初始值為p0=0endl;/初始化P0=p0=0;max1=100;tol=0.001;for(k=2;k=max1;k+)Pk-1=g(Pk-2);err=fabs(Pk-1-Pk-2);relerr=err/(fabs(Pk-1+eps);p=Pk-1;if(errtol)|(relerrtol)break;if(k=max1)cout迭代次數(shù)超過允許的最大迭代次數(shù)！endl;cout不動點的近似值為:pendl; cout程序迭代次數(shù)為:kendl;cout近似值的誤差為:errendl;cout求解不動點近似值的序列:endl;for(i=0;ik;i

37、+)coutsetw(10)Pi;coutendl;return 0;9. 拉格朗日插值9.1、 算法說明對每個固定的k,拉格朗日系數(shù)多項式具有性質(zhì)為：. （9-1）9.2、 拉格朗日插值算法程序調(diào)試首先編寫好程序，然后編譯鏈接，運行程序，按照程序提示依次輸入點的個數(shù)、點數(shù)對應(yīng)的x值、y值、縱坐標、再輸入x值，輸入結(jié)果如下：圖9-1拉格朗日插值算法程序調(diào)試圖9-2拉格朗日插值算法程序計算結(jié)果9.3、 拉格朗日插值算法程序代碼：#includeusing namespace std;double func(double X,int k,double x,int n);int main()doub

38、le Sn=0; int n;coutn;double*x=(double*)malloc(n*sizeof(double); double*y=(double*)malloc(n*sizeof(double);double X;int i;for(i=0;in;i+)cout請輸入xi+1,yi+1:xiyi;coutX;for(i=0;in;i+) Sn=Sn+func(X,i,x,n)*yi;cout通過拉格朗日插值公式所得的結(jié)果：當(dāng)x=X時，y=Snendl; return 0;double func(double X,int k,double x,int n)int i;double

39、 Pn=1;double p;for(i=0;in;i+)if(i=k) continue;else p=(X-xi)/(xk-xi);Pn=Pn*p;return Pn;10、 雅克比迭代10.1、 雅克比迭代的基本算法說明（1），建立判定條件來判斷雅可比迭代是否收斂，因此我們定義矩陣的嚴格對角優(yōu)勢：，其中，可知可以收斂。（2）問題重點方程可以表示成下面的形式： （10-1）這樣就可以提出下列雅可比迭代過程： （10-2）代入初始點（x0,y0,z0）進行迭代。10.2、 雅克比迭代的程序運行圖10-1雅克比迭代的程序運行10.3、 雅克比迭代的程序代碼#include #include using namespace std;#define delta 0.0000001#define n 3#define