典型数值算法的说明 C++语言程序的设计开发.docx
《典型数值算法的说明 C++语言程序的设计开发.docx》由会员分享，可在线阅读，更多相关《典型数值算法的说明 C++语言程序的设计开发.docx（44页珍藏版）》请在冰豆网上搜索。
典型数值算法的说明C++语言程序的设计开发
数值方法课程设计说明书
题目：
典型数值算法的C++语言程序设计
***************************************
学号：
************
院（系）：
理学院
专业：
数学与应用数学091班
***********************
2011年6月15日
陕西科技大学
数值计算课程设计任务书
理学院应用数学专业数学091班级学生：
晏瑞
题目：
典型数值算法的C++语言程序设计
课程设计从2011年5月20日起到2011年6月25日
1、课程设计的内容和要求（包括原始数据、技术要求、工作要求等）：
每人需作10个算法的程序、必做6题、自选4题。
对每个算法要求用C++语言进行编程。
必选题：
1、经典四阶龙格库塔法解一阶微分方程组
2、高斯列主元法解线性方程组
3、牛顿法解非线性方程组
4、龙贝格求积分算法
5、三次样条插值算法（压紧样条）用C++语言进行编程计算
依据计算结果，用Matlab画图并观察三次样条插值效果。
6、M次多项式曲线拟合，据计算结果，用Matlab画图并观察拟合效果。
自选题：
自选4道其他数值算法题目.每道题目重选次数不得超过5次.
2、对课程设计成果的要求〔包括图表、实物等硬件要求〕：
1）提交课程设计报告
按照算法要求，用C++语言设计和开发应用程序，提交由算法说明；程序设计说明；系统技术文档（包括系统各模块主要流程图，软件测试方案与测试记录、软件调试和修改记录、测试结论、运行情况记录），系统使用说明书，源程序代码为附录构成的课程设计报告。
2）课程设计报告版式要求
打印版面要求：
A4纸，页边距：
上2cm，下2cm，左2.5cm、右2cm；字体：
正文宋体、小四号；行距：
固定值20；页眉1.5cm，页脚1.75cm；页码位于页脚居中打印；奇数页页眉“数值计算课程设计”，偶数页页眉“算法名称”，页眉宋体小5号；段落及层次要求：
每节标题以四号黑体左起打印（段前段后各0.5行），节下为小节，以小四号黑体左起打印（段前段后各0.5行）。
换行后以小四号宋体打印正文。
节、小节分别以1、1.1、1.1.1依次标出，空一字符后接各部分的标题。
当论文结构复杂，小节以下的标题，左起顶格书写，编号依次用
（1）、
（2）……或1）、2）……顺序表示。
字体为小四号宋体。
对条文内容采用分行并叙时，其编号用（a）、（b）……或a）、b）……顺序表示，如果编号及其后内容新起一个段落，则编号前空两个中文字符。
3）设计报告装订顺序与规范
封面
数值计算课程设计任务书
数值计算设计课程设计报告正文
设计体会及今后的改进意见
参考文献（资料）
左边缘装订
3、课程设计工作进度计划：
时间
设计任务及要求
第16周
编写和调试程序并按要求撰写设计报告
指导教师：
日期：
教研室主任：
日期：
1.经典四阶龙格库塔法解一阶微分方程
1.1算法说明：
龙格-库塔（Runge-Kutta）方法是一种在工程上应用广泛的高精度单步算法。
由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。
该算法是构建在数学支持的基础之上的。
对于一精度的欧拉公式有：
yi+1=yi+h*K1
K1=f（xi,yi）
当用点xi处的斜率近似值K1与右端点xi+1处的斜率K2的算术平均值作为平均斜率K*的近似值，那么就会得到二阶精度的改进欧拉公式：
yi+1=yi+h*（K1+K2）/2
K1=f（xi,yi）
K2=f（xi+h,yi+h*K1）
依次类推，如果在区间[xi,xi+1]内多预估几个点上的斜率值K1、K2、……Km，并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。
经数学推导、求解，可
以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法：
yi+1=yi+h*（K1+2*K2+2*K3+K4）/6
K1=f（xi,yi）
K2=f（xi+h/2,yi+h*K1/2）
K3=f（xi+h/2,yi+h*K2/2）
K4=f（xi+h,yi+h*K3）
通常所说的龙格-库塔法是指四阶而言的，我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式
1.2，算法程序：
#include
#include
#defineM10
usingnamespacestd;
intmain（）
{
longdoublefeval（longdouble,longdouble,longdouble）;
longdoublefeval1（longdouble,longdouble,longdouble）;
longdoublef[4],g[4];
longdoubleh,a,b,xa,ya;
longdoublex[M+1],y[M+1],T[M+1];
cout<<"请输入区间左右端点a,b:
";
cin>>a>>b;
cout<//xa,ya是初值
cout<<"请输入函数在左端点的初值xa,ya:
";
cin>>xa>>ya;
cout<//h是步长
h=（b-a）/M;
x[0]=xa;
y[0]=ya;
for（inti=0;i<=M;i++）
{
T[i]=a+h*i;//给T赋值
}
for（intk=0;k<=M;k++）
{
f[0]=feval（T[k],x[k],y[k]）;
g[0]=feval1（T[k],x[k],y[k]）;
f[1]=feval（T[k]+h/2,x[k]+h/2*f[0],y[k]+h/2*g[0]）;
g[1]=feval1（T[k]+h/2,x[k]+h/2*f[0],y[k]+h/2*g[0]）;
f[2]=feval（T[k]+h/2,x[k]+h/2*f[1],y[k]+h/2*g[1]）;
g[2]=feval1（T[k]+h/2,x[k]+h/2*f[1],y[k]+h/2*g[1]）;
f[3]=feval（T[k]+h,x[k]+h*f[2],y[k]+h*g[2]）;
g[3]=feval1（T[k]+h,x[k]+h*f[2],y[k]+h*g[2]）;
x[k+1]=x[k]+h/6*（f[0]+2*f[1]+2*f[2]+f[3]）;
y[k+1]=y[k]+h/6*（g[0]+2*g[1]+2*g[2]+f[3]）;
}
cout<for（k=0;k<=M;k++）
{
cout<}
return0;
}
longdoublefeval（longdoublet,longdoublex,longdoubley）
{
longdoublef;f=x+2*y;returnf;
}
longdoublefeval1（longdoublet,longdoublex,longdoubley）
{
longdoublef;f=3*x+2*y;returnf;
}
1.3，运行结果
2.高斯列主元法解线性方程组
2.1算法说明：
首先将线性方程组做成增广矩阵，对增广矩阵进行行变换。
对于元素
，在第i列中，第i行及以下的元素选取绝对值最大的元素，将该最大元素所在的行与第i行交换，然后采用高斯消元法用新得到的
所在的第i行消去第i行以下的元素。
依次进行直到
。
从而得到上三角矩阵。
2.2算法程序：
本程序包含enter.h，remove.h，judge.h及main.cpp这四个函数
1，main.cpp函数
#include
#include
#include"enter.h"//包含名为enter的头函数
#include"remove.h"//包含名为remove的头函数
#include"judge.h"//包含名为judge的头函数
intmain（）
{
introw,col;
cout<<"请输入方程的个数:
";
cin>>row;
cout<<"请输入未知数的个数加一:
";
cin>>col;
doublea[100][100],**p1,*p[100];
for（inti=0;ip[i]=a[i];
p1=p;
enter（p1,row,col）;//调用enter函数
remove（p1,row,col）;//调用remove函数
intnumber=judge（p1,row,col）;//调用judge函数
if（number==0）
{
cout<<"------------------结果----------------"<cout<<"此方程组无解"<}
if（number==-1）
{
cout<<"------------------结果----------------"<cout<<"此方程组有无穷多解"<}
if（number==1）
{
cout<<"------------------结果----------------"<cout<<"此方程组有唯一解"<doublex[100];
x[col-2]=a[row-1][col-1]/a[row-1][row-1];//回代法求解上三角矩阵
for（i=col-3;i>=0;i--）
{
doublet=a[i][col-1];
for（intj=col-2;j>i;j--）
t=t-a[i][j]*x[j];
x[i]=t/a[i][i];
}
for（i=0;icout<<"x"<cout<}
return0;
}
2，enter.h函数
//输入系数矩阵
voidenter（double**q,intm,intn）
{
cout<<"请按行输入未知数前面的系数和等式右边的常数：
"<for（inti=0;ifor（intj=0;jcin>>q[i][j];
}
3，remove.h函数
//寻找列主元，并移动该行，最后化为上三角矩阵
voidremove（double**q,intm,intn）
{
intmin=m;
if（n-1min=n-1;
for（inti=0;i{
intk=i;
doublemax=q[i][i];
for（intj=i+1;jif（fabs（q[j][i]）>fabs（max））
{
max=q[j][i];
k=j;
}//找到第i列从a[i][i]开始的绝对值最大的元素
if（k!
=i）
for（intj=0;j{
doublemat=q[i][j];
q[i][j]=q[k][j];
q[k][j]=mat;
}//通过换行以保证主对角线上的元素是其所在位置及以下元素中绝对值最大的一个
intt=0;
for（j=i;j{
if（q[j][i]==0）
t++;
}
//在消元前判断a[i][i]及其所在列以下元素是否都为零，不都为零再进行消元
if（t!
=m-i）
{
for（intj=i+1;jfor（intk=n-1;k>=i;k--）
q[j][k]=q[j][k]-q[i][k]*q[j][i]/q[i][i];
}
}
for（i=0;ifor（intj=i+1;jq[j][i]=0;//保证经消元后得到上三角阵
}
4，judge.h函数
//用增广矩阵的秩判断该矩阵所对应的方程组的解的情况
intjudge（double**q,intm,intn）
{
intr1=m;
for（inti=0;i{
intg=0;
for（intj=0;jif（q[i][j]==0）
g++;
if（g==n）
r1=r1-1;//求出增广矩阵的秩
}
intr2=m;
for（i=0;i{
inth=0;
for（intj=0;jif（q[i][j]==0）
h++;
if（h==n-1）
r2=r2-1;//求出系数矩阵的秩
}
if（r1==n-1&&r2==n-1）
return1;
elseif（r1==r2&&r1return-1;//据不同的秩的情况返回不同的值
elsereturn0;
}
2.3运行结果：
3.牛顿法解非线性方程组
3.1算法说明
设
已知。
第1步：
计算函数
第2步：
计算雅可比矩阵
第3步：
求线性方程组
的解
。
第4步：
计算下一点
重复上述过程。
3.2算法程序：
#include
#include
#defineN3//非线性方程组中方程个数、未知量个数
#defineEpsilon0.0001//差向量1范数的上限
#defineMax3//最大迭代次数
usingnamespacestd;
constintN2=2*N;
intmain（）
{
voidff（floatxx[N],floatyy[N]）;//计算向量函数的因变量向量yy[N]
voidffjacobian（floatxx[N],floatyy[N][N]）;//计算雅克比矩阵yy[N][N]
voidinv_jacobian（floatyy[N][N],floatinv[N][N]）;//计算雅克比矩阵的逆矩阵inv
voidnewdundiedai（floatx0[N],floatinv[N][N],floaty0[N],floatx1[N]）;//由近似解向量x0计算近似解向量x1
floatx0[N]={0,0,0},y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errornorm;
inti,j,iter=0;
//如果取消对x0的初始化，撤销下面两行的注释符，就可以由键盘向x0读入初始近似解向量
//for（i=0;i//cin>>x0[i];
cout<<"初始近似解向量：
"<for（i=0;icout<cout<do
{
iter=iter+1;
cout<<"第"<//计算向量函数的因变量向量y0
ff（x0,y0）;
//计算雅克比矩阵jacobian
ffjacobian（x0,jacobian）;
//计算雅克比矩阵的逆矩阵invjacobian
inv_jacobian（jacobian,invjacobian）;
//由近似解向量x0计算近似解向量x1
newdundiedai（x0,invjacobian,y0,x1）;
//计算差向量的1范数errornorm
errornorm=0;
for（i=0;i{
errornorm=errornorm+fabs（x1[i]-x0[i]）;
if（errornormbreak;
}
for（i=0;i{
x0[i]=x1[i];
}
}while（iterreturn0;
}
voidff（floatxx[N],floatyy[N]）
{
floatx,y,z;
inti;
x=xx[0];
y=xx[1];
z=xx[2];
yy[0]=x*x-x+y*y+z*z-5;
yy[1]=x*x+y*y-y+z*z-4;
yy[2]=x*x+y*y+z*z+z-6;
cout<<"向量函数的因变量向量是：
"<for（i=0;i{
cout<cout<}
cout<}
voidffjacobian（floatxx[N],floatyy[N][N]）
{
floatx,y,z;
inti,j;
x=xx[0];
y=xx[1];
z=xx[2];
//jacobianhaven*nelement
yy[0][0]=2*x-1;
yy[0][1]=2*y;
yy[0][2]=2*z;
yy[1][0]=2*x;
yy[1][1]=2*y-1;
yy[1][2]=2*z;
yy[2][0]=2*x;
yy[2][1]=2*y;
yy[2][2]=2*z+1;
cout<<"雅克比矩阵是：
"<for（i=0;i{
for（j=0;jcout<cout<}
cout<}
voidinv_jacobian（floatyy[N][N],floatinv[N][N]）
{floataug[N][N2],L;
inti,j,k;
cout<<"开始计算雅克比矩阵的逆矩阵：
"<for（i=0;i{for（j=0;jaug[i][j]=yy[i][j];
for（j=N;jif（j==i+N）aug[i][j]=1;
elseaug[i][j]=0;
}
for（i=0;i{for（j=0;jcout<cout<}
cout<for（i=0;i{
for（k=i+1;k{L=-aug[k][i]/aug[i][i];
for（j=i;jaug[k][j]=aug[k][j]+L*aug[i][j];
}
}
for（i=0;i{for（j=0;jcout<cout<}
cout<for（i=N-1;i>0;i--）
{
for（k=i-1;k>=0;k--）
{L=-aug[k][i]/aug[i][i];
for（j=N2-1;j>=0;j--）
aug[k][j]=aug[k][j]+L*aug[i][j];
}
}
for（i=0;i{for（j=0;jcout<cout<}
cout<for（i=N-1;i>=0;i--）
for（j=N2-1;j>=0;j--）
aug[i][j]=aug[i][j]/aug[i][i];
for（i=0;i{for（j=0;jcout<cout<for（j=N;jinv[i][j-N]=aug[i][j];
}
cout<cout<<"雅克比矩阵的逆矩阵：
"<for（i=0;i{for（j=0;jcout<cout<}
cout<}
voidnewdundiedai（floatx0[N],floatinv[N][N],floaty0[N],floatx1[N]）
{
inti,j;
floatsum=0;
for（i=0;i{sum=0;
for（j=0;jsum=sum+inv[i][j]*y0[j];
x1[i]=x0[i]-sum;
}
cout<<"近似解向量：
"<for（i=0;icout<cout<}
3.3运行结果：
结果为：
x=-1.36628;y=-0.366281;z=2.36628
4.龙贝格求积分算法
4.1算法说明
生成
的逼近表
，并以
为最终解来逼近积分
升级会员 