典型数值算法的说明 C++语言程序的设计开发Word文件下载.docx
《典型数值算法的说明 C++语言程序的设计开发Word文件下载.docx》由会员分享，可在线阅读，更多相关《典型数值算法的说明 C++语言程序的设计开发Word文件下载.docx（44页珍藏版）》请在冰豆网上搜索。
固定值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
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
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<
iostream>
iomanip>
#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;
endl;
//xa,ya是初值
请输入函数在左端点的初值xa,ya:
xa>
ya;
//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<
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]）;
setw（12）<
T"
x"
y"
for（k=0;
T[k]<
x[k]<
y[k]<
return0;
longdoublefeval（longdoublet,longdoublex,longdoubley）
longdoublef;
f=x+2*y;
returnf;
}
longdoublefeval1（longdoublet,longdoublex,longdoubley）
f=3*x+2*y;
1.3，运行结果
2.高斯列主元法解线性方程组
2.1算法说明：
首先将线性方程组做成增广矩阵，对增广矩阵进行行变换。
对于元素
，在第i列中，第i行及以下的元素选取绝对值最大的元素，将该最大元素所在的行与第i行交换，然后采用高斯消元法用新得到的
所在的第i行消去第i行以下的元素。
依次进行直到
。
从而得到上三角矩阵。
2.2算法程序：
本程序包含enter.h，remove.h，judge.h及main.cpp这四个函数
1，main.cpp函数
iostream.h>
#include<
math.h>
#include"
enter.h"
//包含名为enter的头函数
remove.h"
//包含名为remove的头函数
judge.h"
//包含名为judge的头函数
intmain（）
{
introw,col;
cout<
请输入方程的个数:
cin>
row;
请输入未知数的个数加一:
col;
doublea[100][100],**p1,*p[100];
for（inti=0;
p[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<
------------------结果----------------"
此方程组无解"
//无解时予以提示
}
if（number==-1）
此方程组有无穷多解"
//有无穷多解时予以提示
if（number==1）
此方程组有唯一解"
//有唯一解时，解出方程组的解
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;
col-1;
i+1<
="
x[i]<
\t"
//依次输出方程组的解
cout<
return0;
}
2，enter.h函数
//输入系数矩阵
voidenter（double**q,intm,intn）
{
请按行输入未知数前面的系数和等式右边的常数：
for（inti=0;
m;
for（intj=0;
j<
n;
j++）
cin>
q[i][j];
}
3，remove.h函数
//寻找列主元，并移动该行，最后化为上三角矩阵
voidremove（double**q,intm,intn）
intmin=m;
if（n-1<
m）
min=n-1;
min-1;
intk=i;
doublemax=q[i][i];
for（intj=i+1;
if（fabs（q[j][i]）>
fabs（max））
{
max=q[j][i];
k=j;
}//找到第i列从a[i][i]开始的绝对值最大的元素
if（k!
=i）
for（intj=0;
{
doublemat=q[i][j];
q[i][j]=q[k][j];
q[k][j]=mat;
}//通过换行以保证主对角线上的元素是其所在位置及以下元素中绝对值最大的一个
intt=0;
for（j=i;
if（q[j][i]==0）
t++;
}
//在消元前判断a[i][i]及其所在列以下元素是否都为零，不都为零再进行消元
if（t!
=m-i）
for（intj=i+1;
for（intk=n-1;
k>
=i;
k--）
q[j][k]=q[j][k]-q[i][k]*q[j][i]/q[i][i];
n-2;
q[j][i]=0;
//保证经消元后得到上三角阵
4，judge.h函数
//用增广矩阵的秩判断该矩阵所对应的方程组的解的情况
intjudge（double**q,intm,intn）
intr1=m;
intg=0;
if（q[i][j]==0）
g++;
if（g==n）
r1=r1-1;
//求出增广矩阵的秩
intr2=m;
inth=0;
n-1;
h++;
if（h==n-1）
r2=r2-1;
//求出系数矩阵的秩
if（r1==n-1&
&
r2==n-1）
return1;
elseif（r1==r2&
r1<
n-1）
return-1;
//据不同的秩的情况返回不同的值
elsereturn0;
2.3运行结果：
3.牛顿法解非线性方程组
3.1算法说明
设
已知。
第1步：
计算函数
第2步：
计算雅可比矩阵
第3步：
求线性方程组
的解
第4步：
计算下一点
重复上述过程。
3.2算法程序：
cmath>
#defineN3//非线性方程组中方程个数、未知量个数
#defineEpsilon0.0001//差向量1范数的上限
#defineMax3//最大迭代次数
constintN2=2*N;
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;
N;
//cin>
x0[i];
初始近似解向量：
for（i=0;
cout<
x0[i]<
"
do
iter=iter+1;
第"
iter<
次迭代开始"
//计算向量函数的因变量向量y0
ff（x0,y0）;
//计算雅克比矩阵jacobian
ffjacobian（x0,jacobian）;
//计算雅克比矩阵的逆矩阵invjacobian
inv_jacobian（jacobian,invjacobian）;
newdundiedai（x0,invjacobian,y0,x1）;
//计算差向量的1范数errornorm
errornorm=0;
errornorm=errornorm+fabs（x1[i]-x0[i]）;
if（errornorm<
Epsilon）
break;
for（i=0;
x0[i]=x1[i];
}
}while（iter<
Max）;
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;
向量函数的因变量向量是：
for（i=0;
yy[i]<
voidffjacobian（floatxx[N],floatyy[N][N]）
inti,j;
//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;
雅克比矩阵是：
for（j=0;
yy[i][j]<
voidinv_jacobian（floatyy[N][N],floatinv[N][N]）
{floataug[N][N2],L;
inti,j,k;
开始计算雅克比矩阵的逆矩阵：
{for（j=0;
aug[i][j]=yy[i][j];
for（j=N;
N2;
if（j==i+N）aug[i][j]=1;
elseaug[i][j]=0;
aug[i][j]<
for（k=i+1;
{L=-aug[k][i]/aug[i][i];
for（j=i;
aug[k][j]=aug[k][j]+L*aug[i][j];
for（i=N-1;
0;
{
for（k=i-1;
for（j=N2-1;
for（i=N-1;
for（j=N2-1;
aug[i][j]=aug[i][j]/aug[i][i];
inv[i][j-N]=aug[i][j];
雅克比矩阵的逆矩阵：
inv[i][j]<
voidnewdundiedai（floatx0[N],floatinv[N][N],floaty0[N],floatx1[N]）
floatsum=0;
for（i=0;
{sum=0;
for（j=0;
sum=sum+inv[i][j]*y0[j];
x1[i]=x0[i]-sum;
近似解向量：
x1[i]<
3.3运行结果：
结果为：
x=-1.36628;
y=-0.366281;
z=2.36628
4.龙贝格求积分算法
4.1算法说明
生成
的逼近表
，并以
为最终解来逼近积分
升级会员 