数值计算课程设计程序代码

1、这是本人做的数值计算课程设计的程序源代码。代码均为C+语言编写。前面十个程序代码分别为十种数值算法，最后一个程序代码头文件“的源代码。这十个程序有大局部要用到头文件“，具体请查看程序的头文件包含申明语句。这个头文件中包含矩阵常用运算的函数，可像matlab中一样方便的使用矩阵。1. 二分法解非线性方程#include#includeusing namespace std;/二分法解非线性方程double F(double x) /待解的f(x) return x*sin(x)-1;int main()double a,b,ya,yb,delta,t,yt,err;coutabdelta;int

2、 i,max1;ya=F(a); yb=F(b);if(ya*yb0)return 0;max1=1+int(log(b-a)-log(delta)/log(2);for(i=1;i0)b=t;yb=yt;elsea=t;ya=yt;if(b-adelta)break;t=(a+b)/2;err=fabs(a-b);cout方程的解为:tendl;cout误差: err = errendl;return 0;2. 高斯列主元消去法解线性方程组#include#includematrix.husing namespace std;/高斯列主元消去法解线性方程组int main()matrix A

3、,X; /A用来保存线性方程组的增广矩阵,X保存方程组的解int i,j,k,l,n;double temp;coutn;A(n,n+1); cout请输入线性方程组的增广矩阵:n A = endl;A.input();for(i=1;i=n;i+) /行变换为上三角形temp=A(i,i); k=0;for(j=i+1;jfabs(temp)temp=A(j,i);k=j;if(ki)A.M_h(i,k); /选主元后行交换for(j=i+1;j=n;j+)temp=-A(j,i)/A(i,i);A.M_h(temp,i,j);/cout=1;i-) /回代求解temp=0;for(j=1;

4、j=n-i;j+)temp+=A(i,i+j)*X(i+j);X(i)=(A(i,n+1)-temp)/A(i,i);cout线性方程组的解为(从x1到xn顺序排列): nX = endl;X.print(8);return 0;3. 龙贝格积分算法#include#includematrix.husing namespace std;/龙贝格积分算法double F(double x)/被积函数return (x*x+x+1)*cos(x);double power(double a,int n) /a的n次方double temp=1;for(int i=0;in;i+)temp*=a;r

5、eturn temp;int main()matrix R; /R为龙贝格积分表double a,b,tol,err=1,x,h,s;int i,n,m=1,j,k;coutabtoln;h=b-a;R=E(n);R(1,1)=(F(a)+F(b)*h/2;for(i=1;itol;i+)h=h/2;s=0;for(j=1;j=m;j+)x=a+h*(2*j-1);s+=F(x);R(i+1,1)=R(i,1)/2+h*s;m=2*m;for(k=1;k=i;k+)R(i+1,k+1)=R(i+1,k)+(R(i+1,k)-R(i,k)/(power(4,k)-1);err=fabs(R(i,

6、i)-R(i+1,k); /注意:结束循环时，k=i+1。cout误差: err = errendl;cout积分表为: R= n;R.print(8);return 0;4. 牛顿-拉夫森法解非线性方程组#include#includematrix.husing namespace std;/牛顿-拉夫森法解非线性方程组matrix F(matrix P)/非线性方程组matrix fp;double x,y;x=P(1,1);y=P(2,1);fp(1,1)=x*x-2*x-y+0.5;fp(2,1)=x*x+4*y*y-4;return fp;matrix JF(matrix Q)/雅可

7、比矩阵matrix jf;double x,y;x=Q(1,1);y=Q(2,1);jf(1,1)=2*x-2;jf(1,2)=-1;jf(2,1)=2*x;jf(2,2)=8*y;return jf;int main()matrix P,Q,Y,J;double delta,epsilon,err,relerr;int i,max;coutdeltaepsilonmax;P(1,1)=2.00;P(2,1)=0.25;for(i=1;i=max;i+)Y=F(P);J=JF(P);Q=inv(J)*Y;P=P-Q;coutPi=n;P.print(8);err=norm(Q);relerr=

8、err/norm(P);if(errdelta|relerrdelta|(fabs(Y(1,1)epsilon)&fabs(Y(2,1)epsilon)break;cout方程组的解为: P =n;P.print(7);cout迭代次数为: iter = iendl;cout误差: err = errendl;return 0;5. 牛顿-拉夫森法#include#includeusing namespace std;/牛顿-拉夫森法double F(double x) /待解的f(x) return x*sin(x)-1;double dF(double x) /f(x)的导数return

9、sin(x)+x*cos(x);int main()double p0,delta,epsilon,err,temp,relerr;int i,max1;coutp0max1deltaepsilon;for(i=1;i=max1;i+)temp=p0-F(p0)/dF(p0);err=fabs(temp-p0);relerr=2*err/(fabs(temp)+delta);p0=temp; if(errdelta|relerrdelta|fabs(F(p0)epsilon)break;cout迭代次数为: iendl;cout方程的解为: p0endl;cout误差: err = erren

10、dl;return 0;6. 欧拉法解一阶微分方程#include#includematrix.husing namespace std;/欧拉法解一阶微分方程double F(double t,double y)/待求解的函数return (t-y)/2;int main()double a,b,ya,h;int M,i;matrix T,Y,E;coutabyaM;h=(b-a)/M;for(i=1;a+(i-1)*h=b;i+)T(1,i)=a+(i-1)*h;Y(1,1)=ya;for(i=1;i=M;i+)Y(1,i+1)=Y(1,i)+h*F(T(1,i),Y(1,i);cout解

11、得结果如下:endl;E=AUG(T,Y);E.print(7);return 0;#include#includematrix.husing namespace std;/三次紧压样条matrix diff(matrix M_A) /求向量的差分matrix M_T;for(int i=1;iM_A.n;i+)M_T(1,i)=M_A(1,i+1)-M_A(1,i);return M_T;matrix dchu(matrix M_A,matrix M_B) /矩阵的点除matrix M_T;for(int i=1;i=M_A.m;i+)for(int j=1;j=M_A.n;j+)M_T(i

12、,j)=M_A(i,j)/M_B(i,j);return M_T;matrix qu(matrix M_A,int i,int j) /取出向量M_A的第i列至第j列matrix M_T;for(int k=i;k=j;k+)M_T(1,k-i+1)=M_A(1,k);return M_T;int main()matrix X,Y,S,H,D,A,B,C,U,M;double dx0,dxn,temp;int i,N;coutN;X(1,N);Y(1,N);cout请输入所有点的横坐标: X= ;X.input();cout请输入所有点的纵坐标: Y= ;Y.input();coutdx0dx

13、n;N=N-1;H=diff(X);D=dchu(diff(Y),H);A=qu(H,2,N-1);B=qu(H,1,N-1)+qu(H,2,N); B.M_h(2.00,1);C=qu(H,2,N);U=diff(D);U.M_h(6.00,1);B(1,1)-=H(1,1)/2;U(1,1)-=3*(D(1,1)-dx0);B(1,N-1)-=H(1,N)/2;U(1,N-1)-=3*(dxn-D(1,N);for(i=2;i=1;i-)M(1,i+1)=(U(1,i)-C(1,i)*M(1,i+2)/B(1,i);M(1,1)=3*(D(1,1)-dx0)/H(1,1)-M(1,2)/2

14、;M(1,N+1)=3*(dxn-D(1,N)/H(1,N)-M(1,N)/2;for(i=0;i=N-1;i+) /求系数S(i+1,1)=(M(1,i+2)-M(1,i+1)/(6*H(1,i+1);S(i+1,2)=M(1,i+1)/2;S(i+1,3)=D(1,i+1)-H(1,i+1)*(2*M(1,i+1)+M(1,i+2)/6;S(i+1,4)=Y(1,i+1);cout结果如下(系数从左向右按高次向低次排列):n S = n; S.print(8);return 0;8. 经典四阶龙格库塔法#include#includematrix.husing namespace std;

15、/经典四阶龙格库塔法double fun(double t,double y)/求微分的函数return (t-y)/2;int main()double a,b,ya,k1,k2,k3,k4;int M,i;matrix T,Y,R;coutabyaM;double h=(b-a)/M;for(i=1;a+(i-1)*h=b;i+)T(1,i)=a+(i-1)*h;Y(1,1)=ya;for(i=1;i=M;i+)k1=h*fun(T(1,i),Y(1,i);k2=h*fun(T(1,i)+h/2,Y(1,i)+k1/2);k3=h*fun(T(1,i)+h/2,Y(1,i)+k2/2);k

16、4=h*fun(T(1,i)+h,Y(1,i)+k3);Y(1,i+1)=Y(1,i)+(k1+2*k2+2*k3+k4)/6;cout解得结果如下:endl;R=AUG(T,Y);R.print(7);return 0;9. 雅可比迭代解线性方程组#include#includematrix.husing namespace std;/雅可比迭代解线性方程组int main()matrix A,X,P; /A用来保存线性方程组的增广矩阵,X保存方程组的解,P为迭代初始值int i,j,k,n,max1;double temp,delta,err,relerr;coutn;A(n,n+1);

17、cout请输入线性方程组的增广矩阵:n A = endl;A.input(); P(1,n);cout请输入迭代初始近似解: P =;P.input();coutdeltamax1;for(i=1;i=max1;i+)for(j=1;j=n;j+)temp=0;for(k=1;k=n;k+)if(k!=j)temp+=A(j,k)*P(1,k);X(1,j)=(A(j,n+1)-temp)/A(j,j);err=fabs(norm(X-P);P=X; P.M_copy();relerr=err/norm(X);if(errdelta|relerrdelta)break;cout线性方程组的解为

18、(从x1到xn顺序排列): nX = endl;X.print(8);cout误差: err = errendl;cout迭代次数: iendl;return 0;10. 最小二乘多项式拟合#include#includematrix.husing namespace std;/最小二乘多项式拟合double power(double a,int n) /a的n次方double temp=1;for(int i=0;in;i+)temp*=a;return temp;int main()matrix X,Y,A,B,C,F; /C用来保存解得的多项式的系数int i,k,n,M;coutn;X

19、(n);Y(n);cout请输入所有点的横坐标: X= ;X.input();cout请输入所有点的纵坐标: Y= ;Y.input();coutM;for(i=1;i=M+1;i+)for(k=1;k=n;k+)F(k,i)=power(X(k),i-1);/cout测试: F = n;F.print(8);A=(F)*F;B=F*Y;/cout测试: A = n;A.print(8);C=inv(A)*B;cout拟合多项式系数由低次至高次项排列如下: n;C.print(8);return 0;11. 头文件“源代码/新增数乘矩阵及向量定位，对常用运算符进行重载，使矩阵运算使用更加方便简

20、洁#include#include#includeusing namespace std;class m2d /二维链表节点类public:double num;m2d* next_h;m2d* next_l;class matrix /定义矩阵类public:int m; /矩阵行数int n; /矩阵列数m2d* mtr; /矩阵链表头指针matrix()m=0;n=0;mtr=0;double& operator()(int i,int j); /定位到矩阵i行j列的num的引用double& operator()(int i); /定位到向量i行(列)的num的引用m2d* M_p(i

21、nt i,int j); /定位到矩阵i行j列的指针void print(int M_n); /屏幕输出矩阵void input(); /从键盘输入矩阵void M_copy(); /新建动态链表，将原矩阵信息拷贝到其中，防止调用函数时矩阵被修改void M_h(double k,int i,int j); /行变换，将第i行的k倍加到第j行void M_h(double k,int i); /行变换，将第i行乘k倍void M_h(int i,int j); /行变换，交换第i行和第j行void M_delete(); /去除矩阵的动态链表，及时释放内存空间matrix operator +

22、(matrix M_B); /矩阵加法运算matrix operator -(matrix M_B); /矩阵减法运算matrix operator *(matrix M_B); /矩阵相乘private:void add_h(int h); /新增行void add_l(int l); /新增列;m2d* matrix:M_p(int i,int j) /定位到矩阵第i行j列，返回值为该处的指针m2d* temp;if(m=0&(i0)&(j0)mtr=new m2d;m=1;n=1;temp=mtr;if(im)add_h(i-m);if(jn) add_l(j-n);for(int p1

23、=1;p1next_h;for(int p2=1;p2next_l;/cout测试5:m n mtr0)&(j0)mtr=new m2d;m=1;n=1;temp=mtr;if(im)add_h(i-m);if(jn) add_l(j-n);for(int p1=1;p1next_h;for(int p2=1;p2next_l;double &_t=temp-num;return _t;double& matrix:operator()(int i) /定位到向量i行(列)的num的引用if(m=0|(m=1&n=1)return (*this)(1,i);if(n=1&m1)return (

24、*this)(i,n);void matrix:add_h(int h) /新增行int i1,i2;m2d* p;m2d* p1;m2d* p2;p=mtr;for(i1=1;i1next_h;p2=p;for(i2=0;i2next_h=new m2d;p1=p-next_h;p2=p1;for(i1=1;i1next_l=new m2d;p1=p1-next_l;p=p-next_l;p-next_h=p1;p=p2;m=m+1;void matrix:add_l(int l) /新增列int i1,i2;m2d* p;m2d* p1;m2d* p2;p=mtr;for(i1=1;i1n

25、ext_l;p2=p;for(i2=0;i2next_l=new m2d;p1=p-next_l;p2=p1;for(i1=1;i1next_h=new m2d;p1=p1-next_h;p=p-next_h;p-next_l=p1;p=p2;n=n+1;void matrix:print(int M_n) /屏幕输出矩阵,当n=0时按默认格式输出,否那么输出指定位数小数int i1,i2;m2d* p=mtr;m2d* p1;for(i1=1;i1=m;i1+)p1=p;for(i2=1;i2=n;i2+)if(M_n=0)coutsetw(10)num;else coutsetw(M_n+

26、4)setprecision(M_n)num;if(i2!=n)p1=p1-next_l;coutnext_h;void matrix:input() /从键盘输入矩阵int m_h,m_l,i1,i2;m2d* p;m2d* p1;if(m=0&n=0)coutm_hm_l;(*this)(m_h,m_l);cout请输入矩阵:endl;p=mtr;for(i1=1;i1=m;i1+)p1=p;coutt;for(i2=1;i2p1-num;if(i2!=n)p1=p1-next_l;if(i1!=m)p=p-next_h;void matrix:M_copy() /新建动态链表，将原矩阵信

27、息拷贝到其中，防止调用函数时矩阵被修改matrix _M,_T;_T.m=m;_T.n=n;_T.mtr=mtr;int i,j;for(i=1;i=m;i+)for(j=1;j0&j0&i=m&j=m)for(int I=1;I0&i=m)for(int I=1;I0&j0&i=m&j=m&i!=j)for(int I=1;I0;i-)for(int j=n;j0;j-)delete M_p(i,j);m=0;n=0;mtr=0;matrix matrix:operator +(matrix M_B) /矩阵加法运算matrix M_T;if(this-m=M_B.m&this-n=M_B.

28、n)for(int i=1;im;i+)for(int j=1;jn;j+)M_T(i,j)=(*this)(i,j)+M_B(i,j);return M_T;matrix matrix:operator -(matrix M_B) /矩阵减法运算matrix M_T;if(*this).m=M_B.m&(*this).n=M_B.n)for(int i=1;i=(*this).m;i+)for(int j=1;j=(*this).n;j+)M_T(i,j)=(*this)(i,j)-M_B(i,j);return M_T;matrix matrix:operator *(matrix M_B

29、) /矩阵相乘matrix TEMP;double temp;int i,j,k;if(*this).n!=M_B.m)cout错误!前一个矩阵的列数和后一个矩阵的行数不同，无法相乘。n;elsefor(i=1;i=(*this).m;i+)for(j=1;j=M_B.n;j+)temp=0;for(k=1;k=(*this).n;k+)temp=temp+(*this)(i,k)*(M_B(k,j);/if(fabs(temp)1e-15)temp=0; /将太小的数清零TEMP(i,j)=temp;return TEMP;/*以上为matrix类及成员函数定义*/matrix operat

30、or *(double k,matrix M_A) /数乘矩阵matrix M_T; int i,j;for(i=1;i=M_A.m;i+)for(j=1;j=M_A.n;j+)M_T(i,j)=k*M_A(i,j);return M_T;matrix operator (matrix MAT) /求矩阵的转置int I,J;matrix TEMP;for(I=1;I=MAT.m;I+)for(J=1;J=MAT.n;J+)TEMP(J,I)=MAT(I,J);return TEMP;matrix E(int N) /创立N阶单位矩阵int I,J;matrix TEMP;for(I=1;I=N;I+)for(J=1;J=N;J+)if(I=J)TEMP(I,J)=1;else TEMP(I,J)=0;return TEMP;matrix inv(matrix MAT) /求逆矩阵MAT.M_copy();matrix TEMP,T;int i,j,k;double temp,temp2=1;if(MAT.m!=MAT.n)cout错误!矩阵不是