数值计算课程设计

数值计算课程设计1 经典四阶龙格库塔法解一阶微分方程组1.1算法说明：4阶龙格-库塔方法(RK4)可模拟N=4的泰勒方法的精度。

这种算法可以描述为，自初始点00(,)t y 开始，利用12341(22)6k k h f f f f y y ++++=+生成的近似序列，其中1213243(,)(,)22(,)22(,)k k k k k k k k f f t y h hf f t y f h hf f t y f f f t h y hf ==++=++=++ 1.2 用龙格库塔法求解求解微分方程⎩⎨⎧+=+=yx y y x x 232 满足初值条件：()()⎩⎨⎧==4060y x1.3 算法流程图：图1-1 四阶龙格库塔法解一阶微分方程组经典四阶龙格库塔法解一阶微分方程组1.4 程序调试：组建调试，确保程序可运行时输入初值，区间，步长。

1.5 程序运行运行界面及运行结果：图1-2 四阶龙格库塔法解一阶微分方程组1.6 源程序代码:#include<iostream>#include<iomanip>using namespace std;#define N 10float ff(float t,float xx,float yy){xx=xx+2*yy;return xx;}float gg(float t,float xx,float yy){yy=3*xx+2*yy;return yy;}void rks4(float xx[N],float yy[N],float tt[N],double a,double b) { float h,p1,p2,p3,p4,q1,q2,q3,q4;xx[0]=6;yy[0]=4;int i,p;h=(b-a)/N;for(p=0;p<N;p++)数值计算课程设计 tt[p]=a+h*p;for(i=0;i<N;i++){ p1=h*ff(tt[i],xx[i],yy[i]);q1=h*gg(tt[i],xx[i],yy[i]);p2=h*ff(tt[i]+h/2,xx[i]+p1/2,yy[i]+q1/2); q2=h*gg(tt[i]+h/2,xx[i]+p1/2,yy[i]+q1/2); p3=h*ff(tt[i]+h/2,xx[i]+p2/2,yy[i]+q2/2); q3=h*gg(tt[i]+h/2,xx[i]+p2/2,yy[i]+q2/2); p4=h*ff(tt[i]+h,xx[i]+p3,yy[i]+q3);q4=h*gg(tt[i]+h,xx[i]+p3,yy[i]+q3);xx[i+1]=xx[i]+(p1+2*p2+2*p3+p4)/6;yy[i+1]=yy[i]+(q1+2*q2+2*q3+q4)/6;}}int main(){float xx[N],yy[N],tt[N];rks4(xx,yy,tt,0,0.2);int i,j,k;for(i=0;i<N;i++)cout<<setw(5)<<tt[i]<<" ";cout<<endl;for(j=0;j<N;j++)cout<<setw(5)<<xx[j]<<" ";cout<<endl;for(k=0;k<N;k++)cout<<setw(6)<<yy[k]<<" ";cout<<endl;return 0;}高斯列主元法解线性方程组2 高斯列主元法解线性方程组2.1 算法说明：Gauss列主元序消元法主要根据线性方程组人以交换两个方程的次序，方程组的同解解性不变，且解的分量次序也不变。

于是，第k步在顺序学消元法进行之前，从Ak的第k’列元素a[k][k],a[k+1][k],……a[n][k]中选出绝对值最大者，并进行记录所在行，即|a[i][k]|=max|a[i][k]|如果l不等k，则交换矩阵的第k’行与第l列所对应的元素，然后，再进行第k步顺序消元法。

2.2 算法流程图：数值计算课程设计图2-1 高斯列主元法解线性方程组2.3高斯列主元程序调试：对所编写的高斯列主元程序进行编译和链接，然后执行得如下所示的窗口，我们按命令输入增广矩阵的行数为3，输入3行4列的增广矩阵运行。

2.4 程序运行运行界面及运行结果：图2-2 高斯列主元法解线性方程组高斯列主元法解线性方程组2.5 源程序代码:#include<iostream>#include<stdlib.h>#include<cmath >using namespace std;int FindMax(int p,int N,double *A){int i=0,j=0;double max=0.0;for(i=p;i<N;i++){if(fabs(A[i*(N+1)+p])>max){j=i;max=fabs(A[i*(N+1)+p]);}}return j;}void ExchangeRow(int p,int j,double *A,int N){int i=0;double C=0.0;for(i=0;i<N+1;i++){C=A[p*(N+1)+i];A[p*(N+1)+i]=A[j*(N+1)+i];A[j*(N+1)+i]=C;}}void uptrbk(double *A,int N){int p=0,k=0,q=0,j=0;double m=0.0;for(p=0;p<N-1;p++){j=FindMax(p,N,A);ExchangeRow(p,j,A,N);if(A[p*(N+1)+p]==0){cout<<"矩阵是一个奇异矩阵。

没有唯一解。

";break;}for(k=p+1;k<N;k++){m=A[k*(N+1)+p]/A[p*(N+1)+p];for(q=p;q<N+1;q++)A[k*(N+1)+q]=A[k*(N+1)+q]-m*A[p*(N+1)+q];数值计算课程设计}}cout<<endl;cout<<"增广矩阵高斯列主元消去后的矩阵为:"<<endl;for(j=0;j<N*(N+1);j++){if(j%(N+1)==0)cout<<endl;cout<<A[j];}}double* backsub(double *A,int N){double* X=NULL,temp=0.0;int k=0,i=0;X=(double*)malloc(N*sizeof(double));X[N-1]=A[(N-1)*(N+1)+N]/A[(N-1)*(N+1)+N-1];for(k=N-2;k>=0;k--){temp=0.0;for(i=k+1;i<N;i++)temp=temp+A[k*(N+1)+i]*X[i];X[k]=(A[k*(N+1)+N]-temp)/A[k*(N+1)+k];}return X;}int main(){int N=0,i=0;double *A=NULL,*X=NULL;cout<<endl;cout<<"请输入待求解方程组的增广矩阵的行数:";cin>>N;A=(double*)calloc(N*(N+1),sizeof(double));cout<<"请输入待求解方程组的增广矩阵:"<<endl;for(i=0;i<N*(N+1);i++)cin>>A[i];system("cls");cout<<"方程的增广矩阵为:"<<endl;for(i=0;i<N*(N+1);i++){if(i%(N+1)==0)高斯列主元法解线性方程组cout<<endl;cout<<A[i]<<" ";}uptrbk(A,N);X=backsub(A,N);cout<<endl;cout<<"方程组的解为:"<<endl;for(i=0;i<N;i++)cout<<X[i]<<" ";free(A);free(X);cout<<endl;exit(0);}数值计算课程设计3 牛顿法解非线性方程组3.1 算法说明： 设k P 已知。

第1步：计算函数12(,)()(,)k k k k k f p q F P f p q ⎡⎤=⎢⎥⎣⎦第2步：计算雅可比矩阵1122(,)(,)()(,)(,)k k k k k k k k k f p q f p q x yJ P f p q f p q x y ∂∂⎡⎤⎢⎥∂∂⎢⎥=∂∂⎢⎥⎢⎥∂∂⎣⎦第3步：求线性方程组()()k k J P P F P ∆=-的解P ∆。

第4步：计算下一点1k k P P P +=+∆ 重复上述过程。

3.2 用牛顿法解方程组22220.50440x x y x y ⎧--+=⎨+-=⎩ 3.3 算法流程图：图3-1 牛顿法解非线性方程组牛顿法解非线性方程3.4 程序调试：编译组建并运行程序。

3.5 程序运行运行界面及运行结果：图3-2 牛顿法解非线性方程组3.6 源程序代码:#include<iostream>#include<stdlib.h>#include<cmath >#define N 2#define eps 2.2204e-16using namespace std;double* MatrixMultiply(double* J,double Y[]);double*Inv(double*J); double norm(double Q[]);double* F(double X[]);double* JF(double X[]);int method(double* Y,double epsilon);int newdim(double P[],double delta,double epsilon,int max1,double *err) {double *Y=NULL,*J=NULL,Q[2]={0},*Z=NULL,*temp=NULL;double relerr=0.0;int k=0,i=0,iter=0;Y=F(P);for(k=1;k<max1;k++){J=JF(P);temp=MatrixMultiply(Inv(J),Y);for(i=0;i<N;i++)Q[i]=P[i]-temp[i];Z=F(Q);for(i=0;i<N;i++)temp[i]=Q[i]-P[i];*err=norm(temp);relerr=*err/(norm(Q)+eps);for(i=0;i<N;i++)P[i]=Q[i];数值计算课程设计for(i=0;i<N;i++)Y[i]=Z[i];iter=k;if((*err<delta)||(relerr<delta)||method(Y,epsilon)) break;}return iter;}int method(double* Y,double epsilon){if(fabs(Y[0])<epsilon&&fabs(Y[1])<epsilon)return 1;elsereturn 0;}double *MatrixMultiply(double* J,double Y[]){double *X=NULL;int i=0,j=0;X=(double*)malloc(N*sizeof(double));for(i=0;i<N;i++)X[i]=0;for(i=0;i<N;i++)for(j=0;j<N;j++)X[i]+=J[i*N+j]*Y[j];return X;}double *Inv(double *J){double X[4]={0},temp=0.0;int i=0;temp=1/(J[0]*J[3]-J[1]*J[2]);X[0]=J[3];X[1]=-J[1];X[2]=-J[2];X[3]=J[0];for(i=0;i<4;i++)J[i]=temp*X[i];return J;}double norm(double Q[]){double max=0.0;int i=0;for(i=0;i<N;i++){if(Q[i]>max)牛顿法解非线性方程max=Q[i];}return max;}double* F(double X[]){ double x=X[0];double y=X[1];double *Z=NULL;Z=(double*)malloc(2*sizeof(double));Z[0]=x*x-2*x-y+0.5;Z[1]=x*x+4*y*y-4;return Z;}double* JF(double X[]){ double x=X[0];double y=X[1];double *W=NULL;W=(double*)malloc(4*sizeof(double));W[0]=2*x-2;W[1]=-1;W[2]=2*x;W[3]=8*y;return W;}int main(){double P[2]={0};double delta=0.0,epsilon=0.0,err=0.0;int max1=0,iter=0,i=0;cout<<"牛顿法解非线性方程组:\nx^2-2*x-y+0.5=0\nx^2+4*y^2-4=0\n";cout<<"输入的初始近似值x0,y0"<<'\n';for(i=0;i<2;i++)cin>>P[i];cout<<endl;cout<<"请依次输入P的误差限，F(P)的误差限，最大迭代次数"<<'\n';cin>>delta>>epsilon>>max1;iter=newdim(P,delta,epsilon,max1,&err);cout<<"收敛到P的解为:"<<'\n';for(i=0;i<2;i++)cout<<P[i]<<" ";cout<<endl;cout<<"迭代次数为:"<<iter;cout<<endl;cout<<"误差为:"<<err<<endl;return 0;}数值计算课程设计4 龙贝格求积分算法4.1算法说明生成J K ≥的逼近表(,)R J K ，并以(1,1)R J J ++为最终解来逼近积分()(,)ba f x dx R J J ≈⎰逼近(,)R J K 存在于一个特别的下三角矩阵中，第0列元素(,0)R J 用基于2J 个[a,b]子区间的连续梯形方法计算，然后利用龙贝格公式计算(,)R J K 。