数值计算课设.docx

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

2、果然后将这些点进行插值或者拟合后就可以得到微分方程的解。1.4 经典四阶龙格库塔法解一阶微分方程程序代码：#include #include using namespace std;double Runge_Kuta( double (*f)(double x, double y), double x0, double y0, double xn, long 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

3、*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() dou

4、ble f(double x, double y); double x0,y0; double a,b; coutx0y0; coutab; /double x0=0,y0=1; double x,y,step; long 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-3*y)

5、/5; return(r);2. 高斯列主元法解线性方程组2.1 算法说明首先将线性方程组做成增光矩阵，对增广矩阵进行行变换。对第元素，在第i列中，第i行及以下的元素选取绝对值最大的元素，将该元素最大的行与第i行交换，然后采用高斯消元法将新得到的消去第i行以下的元素，一次进行直到，从而得到上三角矩阵。再对得到的上三角矩阵进行回代操作，即可以得到方程组的解。2.2 高斯列主元算法流程图开 始输入（增广矩阵）是否交换中两行是否输出结 束图2-1 高斯列主元法解线性方程组算法流程图2.3 高斯列主元程序调试对所编写的高斯列主元程序进行编译和链接，然后执行得如下所示的窗口，我们按命令输入增广矩阵的行数

6、为3，输入3行4列的增广矩阵，输入时界面图如下，按回车键后，程序执行结果如下图所示：图2-2 运行结果图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(); for(i=0;im;i+) for(j=i;jfabs(aii)?j:i; for(n=0;nm+1;n+) s

7、=ain; ain=acn; acn=s; 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;void load() int i,j; for(i=0;im;i+) for(j=0;jaij;3.牛顿法解非线性方程组3.1 算法说明 设已知。 第1步：计算函数 第2步：计算雅可比矩阵 第3步：求

8、线性方程组 的解。 第4步：计算下一点 重复上述过程。3.2 牛顿法解非线性方程组算法流程图图3-1 算法流程图3.3 牛顿法解非线性方程组算法程序调试以方程组为例，初始近似值x0,y0分别为2.00和0.25，经过10次迭代求出X(1)=4.30017和X(2)=17.6736。运行程序，可以得到结果如下：图3-2 运行结果图3.4 牛顿法解非线性方程组算法程序代码#include#include#define N 2 #define epsilon 0.0001 #define max 10 using namespace std;const int N2=2*N;int main()vo

9、id ff(float xxN,float yyN);void ffjacobian(float xxN,float yyNN);void inv_jacobian(float yyNN,float invNN);void newdim(float 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=ite

10、r+1; cout第 iter 次迭代开始：endl; ff(x0,y0); ffjacobian(x0,jacobian); inv_jacobian(jacobian,invjacobian); newdim(x0, invjacobian,y0,x1); errornorm=0; for (i=0;iN;i+) errornorm=errornorm+fabs(x1i-x0i); if (errornormepsilon) break; for (i=0;iN;i+) x0i=x1i; while (itermax);return 0;void ff(float xxN,float yy

11、N) float x,y; int i; x=xx0; y=xx1; yy0=3*x*x+2*y-4; yy1=2*x+2*y-3; 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*x-2; yy01=-1; yy10=2*x; yy11=8*y;cout雅克比矩阵： endl; for( i=0;iN;i+) for(j=0;jN;j+) coutyyij ; cou

12、tendl; 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+) for(j=0;jN2;j+) coutaugij ; coutendl; coutendl; for (i=0;iN;i+) for (k=i+1;kN;k+) L=-a

13、ugki/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+) coutaugij ; coutendl; cout=0;i-) for(j=N2-1;j=0;j-) augij=augij/augii; for (i=0

14、;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; float sum=0; for(i=0;iN;i+) sum=0; for(j=0;jN;j+) sum=sum+invij*y0j; x

15、1i=x0i-sum; cout近似解向量：endl; for (i=0;iN;i+) coutx1i endl;4.龙贝格求积分算法4.1 算法说明 生成的逼近表，并以为最终解来逼近积分 逼近存在于一个特别的下三角矩阵中，第0列元素用基于个a,b子区间的连续梯形方法计算，然后利用龙贝格公式计算。当时，第行的元素为 当时，程序在第行结束。4.2 龙贝格求积分算法流程图图4-1 算法流程图4.3 龙贝格求积分法程序调试 我们以求解积分方程，精度为0.0001，最高迭代10次为例，对所编写的龙贝格求积分算法程序进行编译和链接，经执行后得如下所示的窗口：图4-2 运行结果4.4 龙贝格求积分算法程序

16、代码#include#include#define f(x) sin(x*x) #define epsilon 0.0001 #define MAXREPT 10 using namespace std;double Romberg(double aa, double bb) int m, n; double h, x; double s, q; double ep; double *y = new doubleMAXREPT;double p; h = bb - aa; y0 = h*(f(aa) + f(bb)/2.0; m = 1; n = 1; ep = epsilon + 1.0;

17、while (ep = epsilon) & (m MAXREPT) p = 0.0; for (int i=0; in; i+) x = aa + (i+0.5)*h; p = p + f(x); p = (y0 + h*p)/2.0; 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; ym-1 = q; n = n + n; h = h/2.0; return (q);int main() double

18、a,b; coutRomberg积分,请输入积分范围a,b:ab; cout积分结果:Romberg(a, b)endl; return 0;5.三次样条插值算法5.1 算法说明表5-1 算法说明策略描述包含和的方程(i)三次紧压样条，确定，（如果导数已知，这是“最佳选择”）(ii)natural三次样条（一条“松弛曲线”）， (iii)外挂到端点(iv)是靠近端点的常量， (v)在每个端点处指定， 5.2 三次样条插值算法（压紧样条）程序调试求解三次紧压样条曲线，以过点（0，0.0），（1，1.5），（2，2.5）和（3，3.5），而且一阶导数边界条件S(0)=0.5和S(2)=-5的三次压

19、紧样条曲线。将所编写的程序进行编译，链接和执行后得如下所示的结果：图5-1 运行结果图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);

20、 void cal_m(int n) 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; for(i = 0

21、; 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(i = 1; i n; i+) f

22、xymi = 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; coutSi+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个点，横坐标是确定的。最小二乘抛物线的系数表示为 求解A，B和C的线性方程组为 6.2 M