数值计算课程设计(陕西科技大学)

1、数值计算课程设计1、 经典四阶龙格库塔法解一阶微分方程1.1、 算法说明龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。4阶龙格-库塔方法(RK4)可模拟N=4的泰勒方法的精度。这种算法可以描述为，自初始点开始，利用下面的计算方法生成近似序列(1-1)1.2、 经典四阶龙格库塔法解一阶微分方程算法流程图图1-1 经典四阶龙格库塔法解一阶微分方程算法流程图1.3、 经典四阶龙格库塔法解一阶微分方程程序调试图1-2 经典四阶龙格库塔法解一阶微分方程程序调试1.4、

2、经典四阶龙格库塔法解一阶微分方程代码#include #include using namespace std;/f为函数的入口地址，x0、y0为初值，xn为所求点，step为计算次数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

3、+h*k1/2); 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()do

4、uble f(double 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

5、r; r=(x-y)/2; return(r); 2、 高斯列主元法解线性方程组2.1、 算法说明首先将线性方程组做成增光矩阵，对增广矩阵进行行变换。对第元素，在第i列中，第i行及以下的元素选取绝对值最大的元素，将该元素最大的行与第i行交换，然后采用高斯消元法将新得到的消去第i行以下的元素。一次进行直到。从而得到上三角矩阵。再对得到的上三角矩阵进行回代操作，即可以得到方程组的解。2.2、 高斯列主元算法代码#include#include#include#define m 10using namespace std;double Amm,bm,Augmm+1;void equ(double A

6、mm,double bm,double xm,int n)int i,i1,j,k; double Augmm+1,maxele,Temp,l,s; for (i=0;in;i+) /构建增广矩阵Aug for (j=0;jn;j+) Augij=Aij; Augin=bi; for (i=0;in;i+) /输出增广矩阵Aug for (j=0;jn;j+) coutAugij ; coutAuginendl;for (i1=0;i1n-1;i1+) /求主元maxele=fabs(Augi1i1);k=i1;for(i=i1;in;i+)if (maxelefabs(Augii1)maxe

7、le=fabs(Augii1); k=i; couti1=i1 k=k maxele=maxele Augii1 endl;for (j=i1;jn+1;j+) /换行Temp=Augi1j;Augi1j=Augkj;Augkj=Temp;for (i=0;in;i+) /换行后输出增广矩阵Aug for (j=0;jn;j+) coutAugij ; coutAuginendl; for (k=i1+1;kn;k+)l=-Augki1/Augi1i1; for (j=i1;jn+1;j+)Augkj=Augkj+l*Augi1j;cout每次消元结束输出增广矩阵Augendl;for (i=

8、0;in;i+) /每次消元结束输出增广矩阵Aug for (j=0;jn;j+) coutAugij ; coutAugin=0;i=i-1)s=0; for (j=i+1;jn;j+) s=s+Augij*xj; couts=sendl; xi=(Augin-s)/Augii; coutxi=xiendl;int main()void equ(double Amm,double bm,double xm,int n); double Amm,bm,xm;int i,j,n;cout输入未知量的个数n; if (nm) cout问题规模太大，需更改原程序中符号常量mendl; return

9、0 ;for (i=0;in;i+)cout请输入A的第i+1行:;for(j=0;jAij; cout请输入b的第i+1bi; equ( A,b, x,n);for (i=0;in;i+) coutxi ; coutendl;return(0);2.3、 高斯列主元程序调试对所编写的高斯列主元程序进行编译和链接，然后执行得如下所示的窗口，我们按命令输入增广矩阵的行数为3，输入3行5列的增广矩阵,运行界面为：图2-2 高斯列主元程序调试3、 牛顿法解非线性方程组3.1、 算法说明设已知。第1步：计算函数 （3-1）第2步：计算雅可比矩阵 （3-2）第3步：求线性方程组 的解。第4步：计算下一点

10、 （3-3）重复上述过程。3.2、 牛顿法解非线性方程组算法程序代码#include#include#include#define N 2 / 非线性方程组中方程个数、未知量个数 #define Epsilon 0.0001 / 差向量1范数的上限#define Max 100 /最大迭代次数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 invN

11、N);void newdundiedai(float x0N, float invNN,float y0N,float x1N);float x0N=2.0,0.25,y0N,jacobianNN,invjacobianNN,x1N,errornorm;int i,j,iter=0;/for( i=0;ix0i;cout初始解向量是：endl;for (i=0;iN;i+) coutx0i ; coutendl;coutendl;do iter=iter+1;cout第 iter 次迭代开始endlendl; /计算向量函数值-因变量向量存放在一维数组y0中，y0做函数ff的实参，ff(x0,

12、y0)的调用将传出因变量向量ff(x0,y0); /计算雅克比矩阵，二维数组jacobian作为ffjacobian的实参;调用ffjacobian(x0,jacobian)， /按照地址传送方式，二维数组jacobian将传出雅克比矩阵的计算结果ffjacobian(x0,jacobian); /调用函数inv_jacobian计算雅克比矩阵的逆矩阵，计算结果存放在实参invjacobian中inv_jacobian(jacobian,invjacobian); /由近似解向量x0迭代计算近似解向量x1 cout第 iter 次迭代得到; newdundiedai(x0, invjacobi

13、an,y0,x1); /计算x1与x0的差向量的1范数errornorm=0;for (i=0;iN;i+)errornorm=errornorm+fabs(x1i-x0i);if (errornormEpsilon) break; /为下次循环计算准备新的x0for (i=0;iN;i+)x0i=x1i; while (iteriter;return 0;void ff(float xxN,float yyN)float x,y; int i; x=xx0;y=xx1;yy0=x*x-2*x+y+1;yy1=x*x-2*y*y+4; cout当前自变量向量为：endl;for( i=0;iN

14、;i+) coutsetw(10)setiosflags(ios:right)setprecision(5)setiosflags(ios:fixed)xxi ; coutendl;coutendl; cout当前因变量向量为：endl;for( i=0;iN;i+) coutsetw(10)setiosflags(ios:right)setprecision(5)setiosflags(ios:fixed)yyi ; coutendl;coutendl;void ffjacobian(float xxN,float yyNN) float x,y; int i,j;x=xx0;y=xx1;/

15、jacobian have n*n elementyy00=2*x-2;yy01=1;yy10=2*x;yy11=-4*y; cout当前的雅克比矩阵是： endl;for( i=0;iN;i+) for(j=0;jN;j+) coutsetw(10)setiosflags(ios:right)setprecision(5)setiosflags(ios:fixed)yyij ; coutendl; coutendl;void inv_jacobian(float yyNN,float invNN)float augNN2,L; int i,j,k; cout计算雅克比矩阵的逆矩阵endl;

16、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;iN;i+) for(j=0;jN2;j+) coutsetw(10)setiosflags(ios:right)setprecision(5)setiosflags(ios:fixed)augij ; coutendl;coutendl; for (i=0;iN;i+) for (k=i+1;kN;k+) L=-augki/augii;for(j=i;jN2;j+) augkj=augkj+L*a

17、ugij; for (i=0;iN;i+) for(j=0;jN2;j+) coutsetw(10)setiosflags(ios:right)setprecision(5)setiosflags(ios:fixed)augij ; 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;j+) coutsetw(10)setiosflags(ios:right)setprecision(5)setiosfla

18、gs(ios:fixed)augij ; coutendl;cout=0;i-) for(j=N2-1;j=0;j-)augij=augij/augii;for (i=0;iN;i+) for(j=0;jN2;j+) coutsetw(10)setiosflags(ios:right)setprecision(5)setiosflags(ios:fixed)augij ; coutendl; for(j=N;jN2;j+) invij-N=augij;coutendl;cout雅克比矩阵的逆矩阵是（Inverse of Jacobian is） endl;for (i=0;iN;i+) for

19、(j=0;jN;j+) coutsetw(10)setiosflags(ios:right)setprecision(5)setiosflags(ios:fixed)invij ; coutendl;coutendl;void newdundiedai(float x0N, float invNN,float y0N,float x1N)int 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 ; co

20、utendl;coutendl;#include#include#define N 2 / 非线性方程组中方程个数、未知量个数 #define epsilon 0.0001 / 差向量1范数的上限#define max 10 /最大迭代次数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 x0N, f

21、loat 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 jacobian

22、ffjacobian(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 (errornor

23、mepsilon) 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;yy00=2*

24、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;iN;i+)

25、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;j+) co

26、utaugij ; 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)int i,j;

27、/计算非线性方程组的近似解向量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; 3.3、 牛顿法解非线性方程组算法程序调试应用本程序解方程组初始近似值x0,y0分别为2.00和0.25，经过100次迭代求出X(1)=-1.58582和X(2)=0.58185。4、 龙贝格求积分算法4.1、 算法说明生成的逼近表，并以为最终解来逼近积分 （4-1）逼近存在于一个特别的下三角矩阵中，第0列元素用基于个

28、a,b子区间的连续梯形方法计算，然后利用龙贝格公式计算。当时，第行的元素为 （4-2）当时，程序在第行结束。4.2、 龙贝格求积分算法代码#include#includeusing namespace std;#define f(x) sin(x*x) /举例函数#define epsilon 0.0001 /精度#define MAXREPT 10 /迭代次数,到最后仍达不到精度要求,则输出T(m=10).double Romberg(double aa, double bb) /aa,bb 积分上下限 int m, n;/m控制迭代次数, 而n控制复化梯形积分的分点数. n=2m doub

29、le h, x; double s, q; double ep; /精度要求 double *y = new doubleMAXREPT;/为节省空间,只需一维数组 /每次循环依次存储Romberg计算表的每行元素,以供计算下一行,算完后更新double p;/p总是指示待计算元素的前一个元素(同一行) /迭代初值 h = bb - aa; y0 = h*(f(aa) + f(bb)/2.0; m = 1; n = 1; ep = epsilon + 1.0; /迭代计算 while (ep = epsilon) & (m MAXREPT) /复化积分公式求T2n(Romberg计算表中的第一

30、列),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行元素,根据Romberg计算表本行的前一个元素(p指示), /和上一行左上角元素(yk-1指示)求得. s = 1.0; for (int k=1; k=m; 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

31、; ym-1 = q; n = n + n; h = h/2.0; return (q);int main() double a,b; coutRomberg积分,请输入积分范围a,b:ab; cout积分结果:Romberg(a, b)endl; return 0;图4-1 算法流程图4.3、 龙贝格求积分算法程序调试我们以求解积分，精度为0.0001，最高迭代10次为例，对所编写的龙贝格求积分算法程序进行编译和链接，经执行后得如下所示的窗口图4-2龙贝格求积分算法程序调试说明：应用Romberg算法求在区间上的精度为0.0001的积分为0.494508。5、 三次样条插值算法5.1 三次样

32、条插值算法说明表表5-1 三次样条插值算法说明表策略描述包含和的方程(i)三次紧压样条，确定，（如果导数已知，这是“最佳选择”）(ii)natural三次样条（一条“松弛曲线”），(iii)外挂到端点(iv) 是靠近端点的常量，(v)在每个端点处指定，5.2、 三次样条插值算法（压紧样条）程序调试我们将所编写的程序三次样条插值算法（压紧样条）程序进行调试图5-1三次样条插值算法（压紧样条）程序输入界面、运行结果图5-2三次样条插值算法程序运行结果（a） 图5-2三次样条插值算法程序运行结果 （b）- 43 -运行结果分析： （5-1）作图程序（Matlab）：x1=0:0.01:1;y1=0.

33、06*(x1 - 1).3 + 0.42*(x1 - 0).3 - 0.06*(x1 - 1) + 0.08*(x1 - 0);x2=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)gt

34、ext(S3)图形为：图5-3 三次样条插值算法Matlab作图分析5.3、 三次样条插值算法（压紧样条）代码#include#includeusing 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); /求差分voi

35、d cal_m(int n) /用追赶法求解出弯矩向量M float Bmax; B0 = c0 / 2; for(int i = 1; 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;

36、 for(i = 0; i n; i+) /求步长 hi = xi+1 - xi; coutt; switch(t) case 1:cout输入 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

37、(i = 1; i n; i+) fxymi = 6 * f(i-1, i, i+1); for(i = 1; i n; i+) ai = hi-1 / (hi + hi-1); ci = 1 - ai; an = hn-1 / (hn-1 + hn); cal_m(n); coutn输出三次样条插值函数：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;cout

38、Si+1 0) cout-t*(x - xi+1)3; else cout-t*(x - xi+1 0) cout + t*(x - 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、算法说明设有N个点，横坐标是确定的。最小二乘抛物线的系数表示为 （6-1）求解A，B和C的线性方程组为 （6-2）6.2、 M次多项式曲线拟合算法流程图图6-1 算法流程图6.3、 M次多项式曲线拟合算法程序调试我们按命令依次输入命令如下命令后，得程序执行结果如