数值计算课程设计.docx
《数值计算课程设计.docx》由会员分享，可在线阅读，更多相关《数值计算课程设计.docx（44页珍藏版）》请在冰点文库上搜索。
数值计算课程设计
1经典四阶龙格库塔法解一阶微分方程组
1.1算法说明：
4阶龙格-库塔方法（RK4）可模拟N=4的泰勒方法的精度。
这种算法可以描述为，自初始点
开始，利用
生成的近似序列，其中
1.2用龙格库塔法求解求解微分方程
满足初值条件：
1.3算法流程图：
图1-1四阶龙格库塔法解一阶微分方程组
1.4程序调试：
组建调试，确保程序可运行时输入初值，区间，步长。
1.5程序运行运行界面及运行结果：
图1-2四阶龙格库塔法解一阶微分方程组
1.6源程序代码:
#include
#include
usingnamespacestd;
#defineN10
floatff（floatt,floatxx,floatyy）
{
xx=xx+2*yy;
returnxx;
}
floatgg（floatt,floatxx,floatyy）
{
yy=3*xx+2*yy;
returnyy;
}
voidrks4（floatxx[N],floatyy[N],floattt[N],doublea,doubleb）
{floath,p1,p2,p3,p4,q1,q2,q3,q4;
xx[0]=6;
yy[0]=4;
inti,p;
h=（b-a）/N;
for（p=0;ptt[p]=a+h*p;
for（i=0;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;
}
}
intmain（）
{
floatxx[N],yy[N],tt[N];
rks4（xx,yy,tt,0,0.2）;
inti,j,k;
for（i=0;icout<cout<for（j=0;jcout<cout<for（k=0;kcout<cout<return0;
}
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
#include
#include
usingnamespacestd;
intFindMax（intp,intN,double*A）
{inti=0,j=0;
doublemax=0.0;
for（i=p;i{if（fabs（A[i*（N+1）+p]）>max）
{
j=i;
max=fabs（A[i*（N+1）+p]）;
}
}
returnj;
}
voidExchangeRow（intp,intj,double*A,intN）
{
inti=0;
doubleC=0.0;
for（i=0;i{
C=A[p*（N+1）+i];
A[p*（N+1）+i]=A[j*（N+1）+i];
A[j*（N+1）+i]=C;
}
}
voiduptrbk（double*A,intN）
{
intp=0,k=0,q=0,j=0;
doublem=0.0;
for（p=0;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{
m=A[k*（N+1）+p]/A[p*（N+1）+p];
for（q=p;qA[k*（N+1）+q]=A[k*（N+1）+q]-m*A[p*（N+1）+q];
}
}
cout<cout<<"增广矩阵高斯列主元消去后的矩阵为:
"<for（j=0;j{
if（j%（N+1）==0）
cout<cout<}
}
double*backsub（double*A,intN）
{
double*X=NULL,temp=0.0;
intk=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;itemp=temp+A[k*（N+1）+i]*X[i];
X[k]=（A[k*（N+1）+N]-temp）/A[k*（N+1）+k];
}
returnX;
}
intmain（）
{
intN=0,i=0;
double*A=NULL,*X=NULL;
cout<cout<<"请输入待求解方程组的增广矩阵的行数:
";
cin>>N;
A=（double*）calloc（N*（N+1）,sizeof（double））;
cout<<"请输入待求解方程组的增广矩阵:
"<for（i=0;icin>>A[i];
system（"cls"）;
cout<<"方程的增广矩阵为:
"<for（i=0;i{
if（i%（N+1）==0）
cout<cout<}
uptrbk（A,N）;
X=backsub（A,N）;
cout<cout<<"方程组的解为:
"<for（i=0;icout<free（A）;
free（X）;
cout<exit（0）;
}
3牛顿法解非线性方程组
3.1算法说明：
设
已知。
第1步：
计算函数
第2步：
计算雅可比矩阵
第3步：
求线性方程组
的解
。
第4步：
计算下一点
重复上述过程。
3.2用牛顿法解方程组
3.3算法流程图：
图3-1牛顿法解非线性方程组
3.4程序调试：
编译组建并运行程序。
3.5程序运行运行界面及运行结果：
图3-2牛顿法解非线性方程组
3.6源程序代码:
#include
#include
#include
#defineN2
#defineeps2.2204e-16
usingnamespacestd;
double*MatrixMultiply（double*J,doubleY[]）;
double*Inv（double*J）;doublenorm（doubleQ[]）;
double*F（doubleX[]）;
double*JF（doubleX[]）;
intmethod（double*Y,doubleepsilon）;
intnewdim（doubleP[],doubledelta,doubleepsilon,intmax1,double*err）
{double*Y=NULL,*J=NULL,Q[2]={0},*Z=NULL,*temp=NULL;
doublerelerr=0.0;
intk=0,i=0,iter=0;
Y=F（P）;
for（k=1;k{J=JF（P）;
temp=MatrixMultiply（Inv（J）,Y）;
for（i=0;iQ[i]=P[i]-temp[i];
Z=F（Q）;
for（i=0;itemp[i]=Q[i]-P[i];
*err=norm（temp）;
relerr=*err/（norm（Q）+eps）;
for（i=0;iP[i]=Q[i];
for（i=0;iY[i]=Z[i];
iter=k;
if（（*errbreak;
}
returniter;
}
intmethod（double*Y,doubleepsilon）
{
if（fabs（Y[0]）return1;
else
return0;
}
double*MatrixMultiply（double*J,doubleY[]）
{
double*X=NULL;
inti=0,j=0;
X=（double*）malloc（N*sizeof（double））;
for（i=0;iX[i]=0;
for（i=0;ifor（j=0;jX[i]+=J[i*N+j]*Y[j];
returnX;
}
double*Inv（double*J）
{
doubleX[4]={0},temp=0.0;
inti=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];
returnJ;
}
doublenorm（doubleQ[]）
{doublemax=0.0;
inti=0;
for（i=0;i{if（Q[i]>max）
max=Q[i];
}
returnmax;
}
double*F（doubleX[]）
{doublex=X[0];
doubley=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;
returnZ;
}
double*JF（doubleX[]）
{doublex=X[0];
doubley=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;
returnW;
}
intmain（）
{doubleP[2]={0};
doubledelta=0.0,epsilon=0.0,err=0.0;
intmax1=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<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<cout<cout<<"迭代次数为:
"<cout<cout<<"误差为:
"<return0;
}
4龙贝格求积分算法
4.1算法说明
生成
的逼近表
，并以
为最终解来逼近积分
逼近
存在于一个特别的下三角矩阵中，第0列元素
用基于
个[a,b]子区间的连续梯形方法计算，然后利用龙贝格公式计算
。
当
时，第
行的元素为
当
时，程序在第
行结束。
4.2求积分方程
4.3算法流程图：
图4-1龙贝格求积分算法
4.4程序调试
编译组建并运行程序。
4.5程序运行运行界面及运行结果
图4-2龙贝格求积分算法
4.6源程序代码
#include
#include
usingnamespacestd;
doublef（doublex）
{returnx*x;
}
intmain（）
{intM=1,n=0,p=0,K=0,i=0,j=0,J=0;
doubleh=0.0,a=0.0,b=0.0,err=1.0,quad=0.0,s=0.0,x=0.0,tol=0.0;
doubleR[30][30]={0};
a=0;
b=1;
h=b-a;
n=4;
tol=0.000001;
cout<<"求解函数y=x*x在（0,1）上的龙贝格矩阵"<<'\n';
cout<<"龙贝格矩阵最大行数为:
"<"<R[0][0]=h*（f（a）+f（b））/2;
while（（（err>tol）&&（J{
J=J+1;
h=h/2;
s=0;
for（p=1;p<=M;p++）
{
x=a+h*（2*p-1）;
s=s+f（x）;
}
R[J][0]=R[J-1][0]/2+h*s;
M=2*M;
for（K=1;K<=J;K++）
{
R[J][K]=R[J][K-1]+（R[J][K-1]-R[J-1][K-1]）/（pow（4,K）-1）;
}
err=fabs（R[J-1][J-1]-R[J][K]）;
}
quad=R[J][J];
cout<<'\n';
cout<<"龙贝格矩阵为:
"<for（i=0;i<（J+1）;i++）
{
for（j=0;j<（J+1）;j++）
{
cout<}
cout<}
cout<cout<<"积分值为:
"<cout<<"误差估计为"<cout<<"使用过的最小步长:
"<return0;
}
5三次样条插值算法（压紧样条）用C++语言进行编程计算依据计算结果，用Matlab画图并观察三次样条插值效果
5.1算法说明
（i）如果导数已知，这是“最佳选择”）三次紧压样条，确定
，
（ii）natural三次样条（一条“松弛曲线”）
，
（iii）外挂
到端点
（iv）
是靠近端点的常量
，
（v）在每个端点处指定
，
5.2程序调试：
编译组建并运行程序。
5.3程序运行运行界面及运行结果：
图5-1样条插值算法
借助Matlab绘制出以上三次压紧样条的函数图像如下所示
表5-1三次压紧样条的函数图像
5.4源程序代码:
#include
#include
#defineMAX4
usingnamespacestd;
double*diff（doubleX[],intn）
{inti=0;
double*H=NULL;
H=（double*）malloc（（n-1）*sizeof（double））;
for（i=1;i<=n-1;i++）
{H[i-1]=X[i]-X[i-1];
}
returnH;
}
double*divide（doubleY[],intN,doubleH[]）
{
inti=0;
double*D=NULL;
D=（double*）malloc（N*sizeof（double））;
for（i=0;i{
D[i]=Y[i]/H[i];
}
returnD;
}
intmain（）
{
doubleX[MAX]={2,4,3,1},Y[MAX]={1,2.5,3.4,4.2},S[MAX][MAX]={0},temp=0.0,M[MAX]={0};
intN=MAX-1,i=0,k=0;
doubleA[MAX-1-2]={0},B[MAX-1-1]={0},C[MAX-1-1]={0};
doubledx0=0.2,dxn=1.0;
double*H=NULL,*D=NULL,*U=NULL;
cout<<"求解经过点（0，0.0），（1，0.5），（2，2.0）和（3，1.5）"<<'\n'<<"而且一阶导数边界条件S'（0）=0.2和S'（3）=-1的三次压紧样条曲线"<H=diff（X,MAX）;
D=divide（diff（Y,MAX）,N,H）;
for（i=1;iA[i]=H[i+1];
for（i=0;iB[i]=2*（H[i]+H[i+1]）;
for（i=1;iC[i]=H[i+1];
U=diff（D,N）;
for（i=0;iU[i]=U[i]*6;
B[0]=B[0]-H[0]/2;
U[0]=U[0]-3*（D[0]-dx0）;
B[N-2]=B[N-2]-H[N-1]/2;
U[N-2]=U[N-2]-3*（dxn-D[N-1]）;
for（k=2;k<=N-1;k++）
{
temp=A[k-2]/B[k-2];
B[k-1]=B[k-1]-temp*C[k-2];
U[k-1]=U[k-1]-temp*U[k-2];
}
M[N-1]=U[N-2]/B[N-2];
for（k=N-2;k>=1;k--）
M[k]=（U[k-1]-C[k-1]*M[k+1]）/B[k-1];
M[0]=3*（D[0]-dx0）/H[0]-M[0]/2;
M[N]=3*（dxn-D[N-1]）/H[N-1]-M[N-1]/2;
for（k=0;k<=N-1;k++）
{
S[k][0]=（M[k+1]-M[k]）/（6*H[k]）;
S[k][1]=M[k]/2;
S[k][2]=D[k]-H[k]*（2*M[k]+M[k+1]）/6;
S[k][3]=Y[k];
}
cout<<"求得的三次压紧样条曲线的矩阵S为:
"<<'\n';
for（i=0;i{for（k=0;k{
cout<}
cout<}
cout<for（i=0;icout<<"在区间（0,1）上的样条为:
"<<"y="<return0;
}
6M次多项式曲线拟合，据计算结果，用Matlab画图并观察拟合效果
6.1算法说明：
设有N个点，横坐
标是确定的。
最小二乘抛物线的系数表示为
求解A，B和C的线性方程组为
6.2待拟合多项式曲线的四个点为：
（1，4）（2，3）（3，2）（4，1）
6.3算法流程图:
图6-1M次多项式曲线拟合
6.4程序调试：
编译组建并运行程序。
6.5程序运行运行界面及运行结果：
图6-2曲线拟合
6.6源程序代码:
#include
#include
#defineMAX20
usingnamespacestd;
voidinv（doubleX[MAX][MAX],intn,doubleE[MAX][MAX]）
{inti=0,j=0,k=0;
doubletemp=0.0;
for（i=0;i{
for（j=0;j<
升级会员 