正文数值计算课程设计正文档.docx
《正文数值计算课程设计正文档.docx》由会员分享，可在线阅读，更多相关《正文数值计算课程设计正文档.docx（43页珍藏版）》请在冰点文库上搜索。
正文数值计算课程设计正文档
1.经典四阶龙格库塔法解一阶微分方程组
1．1经典四阶龙格库塔法解一阶微分方程组的算法说明
写出初值（t0,x0,y0），微分区间[a,b]和微分方程f（t,x,y）和g（t,x,y）,步长h;利用
和
生成近似值序列
；
1．2经典四阶龙格库塔法解一阶微分方程组的程序设计说明
图1-1流程图
1．3经典四阶龙格库塔法解一阶微分方程组的程序的调试及运行过程
例1．1、用RK4方法求解区间[00.2]上一阶微分方程组
且步长h=0.02时的龙格库塔解。
调用程序运行结果如下：
图1-2界面图
所以本题的解如上图所示。
1．4经典四阶龙格库塔法解一阶微分方程组的程序源代码
#include
usingnamespacestd;
intmain（）
{
doublef（doublett,doublexx,doubleyy）;
doubleg（doublett,doublexx,doubleyy）;
doublef1,f2,f3,f4,h,a,b,t0,x0,m,t,x,t1,x1,g1,g2,g3,g4,y0,y,y1;
cout<<"请输入初值t0,x0,y0:
"<cin>>t0>>x0>>y0;
cout<<"请输入区间：
"<cin>>a>>b;
cout<<"请输入步长：
"<cin>>h;
m=（b-a）/h;
cout<<"t="<t=t0;
x=x0;
y=y0;
m=（b-a）/h;
for（m;m>0;m--）
{t1=t+h;
f1=f（t,x,y）;
g1=g（t,x,y）;
f2=f（t+h/2,x+h*f1/2,y+h*g1/2）;
g2=g（t+h/2,x+h*f1/2,y+h*g1/2）;
f3=f（t+h/2,x+h*f2/2,y+h*g2/2）;
g3=g（t+h/2,x+h*f2/2,y+h*g2/2）;
f4=f（t+h,x+h*f3,y+h*g3）;
g4=g（t+h,x+h*f3,y+h*g3）;
x1=x+h*（f1+2*f2+2*f3+f4）/6;
y1=y+h*（g1+2*g2+2*g3+g4）/6;
cout<<"t="<t=t1;
x=x1;
y=y1;
}
return0;
}
doublef（doublett,doublexx,doubleyy）
{
doublet;
t=xx+2*yy;
return（t）;
}
doubleg（doublett,doublexx,doubleyy）
{
doublet;
t=3*xx+2*yy;
return（t）;}
2.高斯列主元法解线性方程组
2．1高斯列主元法解线性方程组的算法说明
首先写出线性方程组的增广矩阵；然后寻找最大主元，并交换主元行；接着开始消元形成三角矩阵，最后进行上三角矩阵回代求解。
2．2高斯列主元法解线性方程组的程序设计说明
图2-1流程图
2．3、高斯列主元法解线性方程组的程序的调试及运行过程
例2．1：
解下面非线性方程组的解，
调用程序运行结果如下：
图2-2界面图
所以方程组的解为：
2．4、高斯列主元法解线性方程组的源程序代码
#include
#include
#include
usingnamespacestd;
intmain（）
{constintn=4;//请在此处修改维数n,本程序适用于系数A为n*n的非奇异矩阵
doublea[n][n+1];
doublex[n];
doubleg,t,p,Q,r,s;
inti,w,k,z;
intj=0,h=0,q=0;
cout<<"请在程序中输入方程的维数n,本次程序适合的方程组维数为n="<cout<<"请输入增广矩阵[AB]：
"<for（i=0;ifor（j=0;jcin>>a[i][j];
cout<<"原始增广矩阵为：
"<for（i=0;i{for（j=0;jcout<cout<}
for（h=0;h{
for（i=h;i{
for（q=0;h+qif（a[h][h]for（z=0;z+h{
s=a[h][h+z];
a[h][h+z]=a[h+q][h+z];
a[h+q][h+z]=s;
}
t=a[i+1][h]/a[h][h];
for（j=0;ja[i+1][j]=-t*a[h][j]+a[i+1][j];//消元
}
cout<<"消元后的增广矩阵"<for（i=0;i{
for（j=0;jcout<cout<}
if（a[h][h]==0）
cout<<"!
!
此方程组的系数矩阵是奇异的，解不唯一，不适用于高斯回代列主消元法求解。
"<cout<<"上三角矩阵的增广矩阵为：
"<for（i=0;i{for（j=0;jcout<cout<for（j=n-1;j>=0;j--）//回代
{if（j{Q=0;
for（w=n-1;w>j;w--）
Q=Q+a[j][w]*x[w];
x[j]=（-Q+a[j][n]）/a[j][j];}
else
x[j]=a[j][j+1]/a[j][j];}
cout<<"方程组的解分别如下:
"<for（k=0;kcout<<"x"<return0;}
3.牛顿法解非线性方程组
3．1牛顿法解非线性方程组的算法说明
第一步：
写出函数向量和雅可比行列式矩阵及初始值
。
第二步：
求出函数向量
=
和雅可比行列式矩阵
=
。
第三步：
建立线性方程组
解出
。
第四步：
比较
是否满足所要求的误差限，若满足就输出
，否则重复第一步至第三步。
3．2牛顿法解非线性方程组的程序设计说明
图3-1流程图
3．3牛顿法解非线性方程组的程序的调试及运行过程
例3．1：
用牛顿迭代法求解非线性方程组，设有非线性方程组
初始值
误差小于0.0001；
调用程序运行结果如下：
图3-2界面图
经过三次迭代后的x=1.90068和y=0.311219为满足题意的方程组的解。
3．4牛顿法解非线性方程组的源程序代码
#include
#include
usingnamespacestd;
#defineN2
voidjuzhen（float（*A）[N]）
{
intk,i,j;floatQ,Z;intK1[N],M1[N];floatB[N],C[N];
for（k=0;k{
Q=0;
for（i=k;i{
for（j=k;j{
if（fabs（*（*（A+i）+j））>fabs（Q））
{
Q=*（*（A+i）+j）;
K1[k]=i;
M1[k]=j;
}
}
}
if（fabs（Q）>=0.001）
{
if（K1[k]!
=k）
{
for（j=0;j{
Z=*（*（A+K1[k]）+j）;
*（*（A+K1[k]）+j）=*（*（A+k）+j）;
*（*（A+k）+j）=Z;
}
}
elseif（M1[k]!
=k）
{
for（i=0;i{
Z=*（*（A+i）+M1[k]）;
*（*（A+i）+M1[k]）=*（*（A+i）+k）;
*（*（A+i）+k）=Z;
}
}
for（j=0;j{
if（j!
=k）
{
B[j]=*（*（A+k）+j）;
C[j]=-*（*（A+j）+k）/Q;}
else{B[j]=1;C[j]=1/Q;}
*（*（A+k）+j）=0;
*（*（A+j）+k）=0;
}
for（i=0;i{
for（j=0;j*（*（A+i）+j）=*（*（A+i）+j）+C[i]*B[j];}
}
}
for（k=N-1;k>=0;k--）
{
if（K1[k]!
=k）
for（i=0;i{
Z=*（*（A+i）+K1[k]）;
*（*（A+i）+K1[k]）=*（*（A+i）+k）;
*（*（A+i）+k）=Z;
}
elseif（M1[k]!
=k）
{
for（j=0;j{
Z=*（*（A+M1[k]）+j）;
*（*（A+M1[k]）+j）=*（*（A+k）+j）;
*（*（A+k）+j）=Z;
}
}
}
}
intmain（）
{
inti,j,g=0;
floatzz=0;
floatX0[N+1],X1[N+1],Y[N+1],F[N+1],A[N][N];
cout<<"请输入初值：
"<for（i=1;i<=N;i++）
{
cin>>X0[i];}
loop1:
F[1]=X0[1]*X0[1]-2*X0[1]-X0[2]+0.5;
F[2]=X0[1]*X0[1]+4*X0[2]*X0[2]-4;
g=g+1;
cout<A[0][0]=2*X0[1]-2;
A[0][1]=-1;
A[1][0]=2*X0[1];
A[1][1]=8*X0[2];
juzhen（A）;
for（i=1;i<=N;i++）
{
Y[i]=0;
for（j=1;j<=N;j++）
Y[i]=Y[i]-A[i-1][j-1]*F[j];
X1[i]=X0[i]+Y[i];
}
for（i=1;i<=N;i++）
zz=fabs（X0[i]-X1[i]）+zz;
if（zz>0.0001）
{
zz=0;
for（i=1;i<=N;i++）
{
X0[i]=X1[i];
cout<}
gotoloop1;
}
for（i=1;i<=N;i++）
cout<F[1]=X0[1]*X0[1]-2*X0[1]-X0[2]+0.5;
F[2]=X0[1]*X0[1]+4*X0[2]*X0[2]-4;
for（i=1;i<=N;i++）
cout<return0;}
4.龙贝格求积分算法
4．1龙贝格求积分的算法说明
写出积分的上下限b,a以及被积函数
，建立一个二维数组R[J][J],在数组中第一列元素R[J][0]可由递归梯形公式求的，其它元素由R[J][K]=R[J][K-1]+（R[J][K-1]-R[J-1][K-1]）/（4^k-1）,而积分
R[J][J],当|R[J][J]-R[J-1][J-1]|4．2龙贝格求积分的程序设计说明
图4-1流程图
4．3龙贝格求积分的程序的调试及运行过程
例4．1：
用龙贝求积分的方法求
，允许误差限eps=0.0001
调用程序运行结果如下：
图4-2界面图
因此本题所求的积分值为0.333333。
4．4龙贝格求积分的源程序代码
#include
#include
usingnamespacestd;
intmain（）
{
constintn=10;
doublef（doublex）;
doubleR[n][n],a,b,eps,h,t;
inti,j,k=1;
cout<<"请输入积分的下限a,上限b及误差限eps:
"<cin>>a>>b>>eps;
h=b-a;
R[0][0]=（h/2）*（f（a）+f（b））;
for（i=1;i{
h=h/2;
t=0;
for（j=1;jt=t+f（a+j*h）;
k++;
R[i][0]=（h/2）*（f（a）+2*t+f（b））;
}
for（j=1;jfor（i=1;iR[i][j]=R[i][j-1]+（R[i][j-1]-R[i-1][j-1]）/（pow（4,j）-1）;
i=0;
do
i++;
while（fabs（R[i][i]-R[i-1][i-1]）cout<<"所求积分的值为：
"<return0;
}
doublef（doublex）
{
doublem;
m=x*x;
return（m）;
}
5.三次样条插值算法
5．1三次样条插值算法的算法说明
针对三次样条的端点约束，第一：
三次紧压样条确定
，
，构造方程
解出
；第二：
nature三次样条
.；第三：
外推
到端点构造方程
解出
；第四：
是靠近端点的常量，
;第五：
在每个端点处指定
，
，接着通过公式
；从而求出三次多项式。
5．2三次样条插值算法的调试及运行
例5．1、求三次紧压样条曲线，经过点（0，0），（1，1.5），（2，2.5），（3，3.5），且一阶导数的边界条件为
。
调用程序运行结果如下：
图5-1界面图
所以S矩阵就为所求插值函数的对应系数。
5．3三次样条插值算法的源程序代码
#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;iD[i]=Y[i]/H[i];
returnD;
}
intmain（）
{doubleX[max]={0,1,2,3},Y[max]={0,0.5,2.0,1.5};
doubles[max][max]={0},temp=0.0,M[max]={0};
intN=max-1,i=0,k=0;
doubleA[max-3]={0},B[max-2]={0},C[max-2]={0};
doubledx0=0.2,dxn=-1.0;
double*H=NULL,*D=NULL,*U=NULL;
cout<<"输入边界条件：
s'（x0）和s'（x3）:
"<cin>>dx0>>dxn;
cout<<"输入过点的x坐标：
"<for（i=0;icin>>X[i];
cout<<"输入过点的y坐标：
"<for（i=0;icin>>Y[i];
H=diff（X,max）;
D=divide（diff（Y,max）,N,H）;
for（i=1;iA[i]=H[i+1];
for（i=1;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[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为：
"<for（i=0;i{
for（k=0;kcout<cout<}
return0;
}
6.M次多项式曲线拟合算法
6．1M次多项式曲线拟合算法的说明
根据泰勒定理任何函数都可以展开成一个和它非常接近的多项式的组合，所以当我们得到函数N个点时，我们就可以通过多项式近似拟合出函数的图像。
6．2M次多项式曲线拟合算法的程序设计说明
图6-1流程图
6．3M次多项式曲线拟合算法的程序调试及运行过程
例6．1、用下面4个数据点进行M次曲线拟合得到多项式函数。
数据如：
（-3,3），（0,1），（2,1），（4,3）。
调用程序运行结果如下：
图6-2界面图
所以拟合曲线为
，用MATLAB画图表现拟合程度如下:
图6-3函数拟合图形
6．4M次多项式曲线拟合算法的源代码
#include
#include
#defineNUM4//节点的个数
#defineN3//拟合的多项式的次数
voidf1（double*,double*）;//将初始x[i]的值的各幂次方存储在一个二维数组里面
voidf2（double*,double*,double[]）;//计算正规方程组的系数矩阵
voidDirectLU（doublea[N][N+1],double[]）;//列主元LU分解
voidswap（double&,double&）;//交换两个变量的值
intmain（）
{
inti,j;
doublex[NUM];//存储原始的节点
doubley[NUM];
cout<<"输入x的值:
"<for（i=0;icin>>x[i];
cout<<"输入y的值:
"<for（i=0;icin>>y[i];
doublea[N][NUM];
doubleb[N][N+1];//正规方程组的系数矩阵和右侧矩阵
doublec[N];
f1（a[0],x）;
f2（a[0],b[0],y）;//计算正规方程组的系数矩阵
DirectLU（b,c）;//列主元LU分解
cout<<"拟合函数的系数分别为:
"<for（i=0;icout<<"a["<return0;
}
voidf1（double*a,double*x）
{
inti,j,k;
doubletemp;
for（i=0;ifor（j=0;j{
temp=1;
for（k=0;ktemp*=x[j];
*（a+i*NUM+j）=temp;
}
}
voidf2（double*a,double*b,doubley[]）
{
inti,j,k;
doubletemp2;
for（i=0;i{
for（j=0;j{
temp2=0;
for（k=0;ktemp2+=*（a+i*NUM+k）*（*（a+j*NUM+k））;
*（b+i*（N+1）+j）=temp2;
}
temp2=0;
for（k=0;k{
temp2+=y[k]*（*（a+i*NUM+k））;
*（b+i*（N+1）+N）=temp2;
}
}
}
voidswap（double&a,double&b）
{
a=a+b;
b=a-b;
a=a-b;
}
voidDirectLU（doublea[N][N+1],doublex[]）
{
inti,r,k,j;
doubles[N],t[N];
doublemax;
for（r=0;r{
max=0;
j=r;
for（i=r;i{
s[i]=a[i][r];
升级会员 