数值计算课程设计..

数值计算课程设计说明书题目：典型数值算法的C++语言程序设计

院（系）：理学院

专业班级：数学1xx

学号： xxxx12010134

学生姓名： x

指导教师： xxxx

2014年7月 6日

陕西科技大学

数值计算课程设计任务书

理学院信息与计算科学/应用数学专业信息12*/数学12* 班级学生：xxx

题目：典型数值算法的C++语言程序设计

课程设计从2014 年 5 月26 日起到2014 年7 月 6 日

1、课程设计的内容和要求（包括原始数据、技术要求、工作要求等）：

每人需作10个算法的程序、必做6题、自选4题。

对每个算法要求用C++语言进行编程。

必选题：

1、高斯列主元法解线性方程组

2、牛顿法解非线性方程组

3、经典四阶龙格库塔法解一阶微分方程组

4、三次样条插值算法（压紧样条）用C++语言进行编程计算

依据计算结果，用Matlab画图并观察三次样条插值效果。

5、龙贝格求积分算法

6、M次多项式曲线拟合，据计算结果，用Matlab画图并观察拟合效果。

自选题：自选4道其他数值算法题目.每道题目重选次数不得超过5次.

2、对课程设计成果的要求〔包括图表、实物等硬件要求〕：

2.1 提交课程设计报告

按照算法要求，应用C++语言设计和开发算法程序，提交由：①每个算法说明；②与算法相应的程序设计说明（程序中的主要变量语义说明，变量的数据类型，数据在内存中组织和存储结构说明，各函数模块的主要流程图，函数功能说明，函数的形参说明，函数的调用方法说明）；③程序调试与实例运行记录(包括程序调试和修改记录、测试结论、运行结果记录)，每个算法的源程序代码编入附录构成的课程设计报告。

2.2 课程设计报告版式要求

目录的要求：居中打印目录二字，（四号黑体，段后1行），字间空一字符；

章、节、小节及其开始页码（字体均为小四号宋体）。节向右缩进两个字符，小节及以后标题均向右缩进四个字符。目录中应包含正文部分每个算法章节标题、设计体总结、无序号的“参考文献资料”，目录的最后一项是“附录”

正文的要求：算法说明论述清楚，公式符号撰写规范，流程图图符规范，计算正确，文字简练通顺，插图简洁规范，书写整洁。文中图、表按制图要求绘制，程序调试和运行情况记录详实。

打印版面要求：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）……顺序表示，如果编号及其后内容新起一个段落，则编号前空两个中文字符。

曲线图表要求：所有曲线、图表、线路图、流程图、程序框图、示意图等不准徒手画，必须按国家规定标准或工程要求绘制（应采用计算机绘图）。

课程设计说明书（报告）中图表、公式要求如下：

（a）图：图的名称采用中文，中文字体为五号宋体，图名在图片下面。引用图应在图题右上角标出文献来源。图号以章为单位顺序编号。格式为：图1-1，空一字符后接图名，比如第1章第5个图是关于高斯列主元法解方程组算法图，图的下方的图号图名应为：图1-5 高斯列主元法解方程组算法图。

（b）表格：表的名称及表内文字采用中文，中文字体为五号宋体，表名在表格上面。表号以章为单位顺序编号，表内必须按规定的符号标注单位。格式为：表1-1，空一字符后接表格名称。

（c）公式：公式书写应在文中另起一行，居中排列。公式序号按章顺序编号。字体为五号宋体，序号靠页面右侧。格式为：（1-1）……。

设计体会及今后的改进意见：设计总结要写出算法理解，编程经验等技术性、学术性总结；体会要简洁、真实、深刻，切忌空话、大话，客套话和矫揉造作之词。改进意见要合理、中肯。

参考文献的要求：另起一页，居中打印参考文献四字（四号黑体，段前段后1行），字间空一字符；另起一行，按论文中参考文献出现的先后顺序用阿拉伯数字连续编号（参考文献应在正文中注出）；参考文献中每条项目应齐全（字体均为小四号宋体）。（格式：[编号]作者.论文或著作名称.期刊名或出版社.出版时间）。（期刊应注明第几期、起止页数（包括论著））。参考文献中条目要符合科技文献引用文献条目书写的国家标准规范。

2.3 设计报告装订顺序与规范

封面

数值计算课程设计任务书

目录

数值计算课程设计报告正文

设计体会及今后的改进意见

参考文献(无需加目录序号)

附录（无需加目录序号)

左边缘装订

3、课程设计工作进度计划：

时间设计任务及要求

第19周编写和调试程序并按要求撰写设计报告

目录

1高斯列主元法解线性方程组 (1)

1.1算法说明 (1)

1.2算例 (1)

1.3程序代码 (1)

2牛顿法解非线性方程组 (2)

2.1算法说明 (2)

2.1算例： (3)

2.2程序代码 (5)

3经典四阶龙格库塔法解一阶微分方程组 (9)

3.1算法说明 (9)

3.2算例 (10)

3.3程序代码 (10)

4三次样条插值算法（压紧样条） (12)

4.1算法说明 (12)

4.2算例 (12)

4.3 C++程序代码 (13)

4.4 matlab程序代码 (16)

5龙贝格求积分算法 (17)

5.1算法说明 (17)

5.2算例 (17)

5.3 程序代码 (17)

6 M次多项式曲线拟合 (18)

6.1算法说明 (18)

6.2算例 (19)

6.3程序代码 (19)

7拉格朗日插值解多项式 (22)

7.1算法说明 (22)

7.2算例 (22)

7.3程序代码 (23)

8二分法求解非线性方程 (24)

8.1算法说明 (24)

8.1算例 (24)

8.2程序代码 (24)

9 不动点迭代 (26)

9.1算法说明 (26)

9.2算例 (26)

9.3程序代码 (26)

10复化梯形求积分公式 (27)

10.1算法说明 (27)

10.2算例 (28)

10.3程序代码 (28)

11设计体会 (29)

参考文献 (30)

1高斯列主元法解线性方程组

1.1算法说明

将线性方程组做成增广矩阵，对增广矩阵进行行变换。对第1列元素，在第i 行及以下的元素中选取绝对值最大的元素，将该元素最大的行与第i 行交换，然后采用高斯消元法将新得到的消去第i 行以下的元素。一次进行直到。从而得到上三角矩阵。再对得到的上三角矩阵进行回代操作，即可以得到方程组的解。 1.2算例

课本99页例3.16，求解方程组的解

运行结果如图1-1所示。

图1-1

1.3程序代码

#include #include

using namespace std;

const N=20; float a[N][N]; int m;

6

2332022428340213424321432143214

321=+++-=+++=+++=+++x x x x x x x x x x x x x x x

x

int main()

{ int i,j;

int c,k,n,p,r;

float x[N],l[N][N],s,d;

cout<<"下面请输入未知数的个数m=";

cin>>m;

cout<

cout<<"请按顺序输入增广矩阵a:"<

for(i=0;i

{

for(j=0;j

cin>>a[i][j];

}

for(i=0;i

{ for(j=i;j

c=(fabs(a[j][i])>fabs(a[i][i]))?j:i; /*找列最大元素*/

for(n=0;n

{s=a[i][n]; a[i][n]=a[c][n]; a[c][n]=s;} /*将列最大数防在对角线上*/ for(p=0;p

cout<

cout<

for(k=i+1;k

{ l[k][i]=a[k][i]/a[i][i];

for(r=i;r

a[k][r]=a[k][r]-l[k][i]*a[i][r]; }

}

x[m-1]=a[m-1][m]/a[m-1][m-1];

for(i=m-2;i>=0;i--)

{ d=0;

for(j=i+1;j

d=d+a[i][j]*x[j];

x[i]=(a[i][m]-d)/a[i][i]; /*求解*/

}

cout<<"该方程组的解为:"<

for(i=0;i

cout<<"x["<

return 0;

}

2牛顿法解非线性方程组

2.1算法说明

P已知，第一步计算函数

设

k

?

?

?

???=),(),()(21k k k k k q p f q p f P F 第二步计算雅可比矩阵

??

???

?

???

???????????=),(),,(),(),,()(2211k k k k k k k k k q p f x q p f x q p f y

q p f x P J 第三步求线性方程组)()(k k P F P P J -=?的解P ?。第四步计算下一点P P P k k ?+=+1重复以上过程。 2.1算例：

设有非线性方程组（课本137页例3.32）

设初始值)25.0,00.2(),(00=q p ，用牛顿法计算),(),,(),,(332211q p q p q p 。

解：迭代3次后的近似解向量为[1.9068,0.311219];具体运行结果如图2-1，2-2，2-3所示。

图2-1

4

405.02022-+=+--=y x y x x

图2-2

图2-1

2.2程序代码

#include

#include

#define N 2 //非线性方程组中方程的个数、未知量的个数

#define Epsilon 0.0001 //差向量的上限

#define Max 100 //最大迭代次数

using namespace std;

const int N2=2*N;

int main()

{

void ff(float xx[N],float yy[N]); //计算向量函数的因变量向量yy[N]

void ffyakebi(float xx[N],float yy[N][N]);//计算雅可比矩阵yy[N][N]

void inv_yakebi(float yy[N][N],float inv[N][N]);//计算雅可比矩阵的逆矩阵inv void newton(float x0[N],float inv[N][N],float y0[N],float x1[N]);//由近似解向量x0求近似解向量x1

float x0[N]={2.0,0.25},y0[N],yakebi[N][N],invyakebi[N][N],x1[N],errornorm;

int i,iter=0;

cout<<"初始近似解向量:"<

for(i=0;i

cin>>x0[i];

for(i=0;i

cout<

cout<

do

{

iter=iter+1;

cout<<"第"<

ff(x0,y0);

ffyakebi(x0,yakebi);

inv_yakebi(yakebi,invyakebi);

newton(x0,invyakebi,y0,x1);

errornorm=0;

for(i=0;i

errornorm=errornorm+fabs(x1[i]);

if(errornorm

for(i=0;i

x0[i]=x1[i];

}while(iter

return 0;

}

void ff(float xx[N],float yy[N])

{

float x,y;

int i;

x=xx[0];

y=yy[1];

yy[0]=x*x-2*x-y+0.5;

yy[1]=x*x+4*y*y-4;

cout<<"向量函数的因变量向量是:"<

for(i=0;i

cout<

cout<

cout<

}

void ffyakebi(float xx[N],float yy[N][N]) {

float x,y;

int i,j;

x=xx[0];

y=xx[1];

yy[0][0]=2*x-2;

yy[0][1]=-1;

yy[1][0]=2*x;

yy[1][1]=8*y;

cout<<"雅可比矩阵是:"<

for(i=0;i

cout<

cout<

}

cout<

}

void inv_yakebi(float yy[N][N],float inv[N][N]) {

float aug[N][N2],L;

int i,j,k;

cout<<"开始计算雅可比矩阵的逆矩阵:"<

for(i=0;i

{for(j=0;j

aug[i][j]=yy[i][j];

for(j=N;j

if(j==i+N)

aug[i][j]=1;

else

aug[i][j]=0;

}

for(i=0;i

{for(j=0;j

cout<

cout<

}

cout<

for(i=0;i

{

for(k=i+1;k

{L=-aug[k][i]/aug[i][i];

for(j=1;j

aug[k][j]=aug[k][j]+L*aug[i][j];} }

{for(j=0;j

cout<

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;j

cout<

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;j

cout<

cout<

for(j=N;j

inv[i][j-N]=aug[i][j];

}

cout<

cout<<"雅可比矩阵的逆矩阵是:"<

for(i=0;i

{ for(j=0;j

cout<

}

cout<

}

void newton(float x0[N],float inv[N][N],float y0[N],float x1[N]) { int i,j; float sum=0; for(i=0;i

for(j=0;j

}

cout<<"近似解向量:"<

cout<

}

3经典四阶龙格库塔法解一阶微分方程组

3.1算法说明

4阶龙格-库塔方法(RK4)可模拟N=4的泰勒方法的精度。这种算法可以描述为，自初始点00(,)t y 开始，利用下面的计算方法生成近似序列

3.2算例：

用4阶龙格-库塔方法计算区间[0.0,0.2]上式（3）的数值解，采用10个子区间，步长h=0.02。（课本400页例9.15）求解运行结果如图3-1所示。

???????+=+=y x dt

dy y x dt

dx

232初值??

?==4)0(6)0(y x

图3-1

3.3程序代码 #include #include using namespace std;

void LG(double (*f)(double t,double x,double y), double (*g)(double t,double x,double y),

double chuzhi[3],double resu[3],double h)

{ double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1; t0=chuzhi[0];x0=chuzhi[1];y0=chuzhi[2]; f1=f(t0,x0,y0); g1=g(t0,x0,y0);

f2=f(t0+h/2,x0+h*f1/2,y0+h*g1/2); g2=g(t0+h/2,x0+h*f1/2,y0+h*g1/2); f3=f(t0+h/2,x0+h*f2/2,y0+h*g2/2);

g3=g(t0+h/2,x0+h*f2/2,y0+h*g2/2);

f4=f(t0+h,x0+h*f3,y0+h*g3);

g4=f(t0+h,x0+h*f3,y0+h*g3);

x1=x0+h*(f1+2*f2+2*f3+f4)/6;

y1=y0+h*(g1+2*g2+2*g3+g4)/6;

resu[0]=t0+h;resu[1]=x1;resu[2]=y1;

}

int main()

{

double f(double t,double x,double y);

double g(double t,double x,double y);

double chuzhi[3],resu[3];

double a,b,S;

double t,step;

int i;

cout<<"输入微分方程的初值:";

cin>>chuzhi[0]>>chuzhi[1]>>chuzhi[2];

cout<<"输入所求微分方程组的微分区间[a,b]:";

cin>>a>>b;

cout<<"输入所求微分方程组所分解子区间的个数step:";

cin>>step;

S=(b-a)/step;

cout<

for(i=0;i

{

LG(f,g,chuzhi,resu,S);

cout<

chuzhi[0]=resu[0];chuzhi[1]=resu[1];chuzhi[2]=resu[2];

}

return 0;

}

double f(double t,double x,double y)

{

double dx; dx=x+2*y; return(dx);

}

double g(double t,double x,double y) { double dy; dy=3*x+2*y; return(dy);

}

4三次样条插值算法（压紧样条）

4.1算法说明

表5-1 三次样条插值算法说明表

策略描述

包含0m 和N m 的方程

(i)三次紧压样条，确定0()S x '，()n S x '（如果导

数已知，这是“最佳选择”）

100003

(())2

m m d S x h '=

-- 1

11

3(())2

N N N N N m m S x d h ---'=

--

(ii)natural 三次样条（一条“松弛曲线”）

00m =，0N m =

(iii)外挂()S x ''到端点

021011

()h m m m m h -=-

11212

()

N N N N N N h m m m m h ------=-

(iv) ()S x ''是靠近端点的常量

01m m =，1N N m m -=

(v)在每个端点处指定()S x ''

00()m S x ''=，()N N m S x ''=

4.2算例

求三次紧压样条曲线，以经过点（0,0），（1,0.5）,(2,2.0)和（3,1.5），且一阶导数的边界条件为S ’(0)=0.2和S ’(3)=-1。（课本222页例5.7） 求解运行结果界面如下：

图4-1 三次样条插值算法（压紧样条）

Matlab 作图如图4-2所示。

图4-2

4.3 C++程序代码 #include #include using namespace std; const int max = 50;

float x[max], y[max], h[max]; float c[max], a[max], fxym[max];

00.51 1.52 2.53

0.5

1

1.5

2

2.5S1

S2

S3

float f(int x1, int x2, int x3)

{

float a=(y[x3]-y[x2])/(x[x3]-x[x2]);

float b=(y[x2]-y[x1])/(x[x2]-x[x1]);

return (a-b)/(x[x3]-x[x1]);

} //求差分

void cal_m(int n)//用追赶法求解出弯矩向量M

{

float B[max];

B[0] = c[0] / 2;

for(int i = 1; i < n; i++)

B[i] = c[i] / (2 - a[i]*B[i-1]);

fxym[0] = fxym[0] / 2;

for(i = 1; i <= n; i++)

fxym[i] = (fxym[i] - a[i]*fxym[i-1]) / (2 - a[i]*B[i-1]); for(i = n-1; i >= 0; i--)

fxym[i] = fxym[i] - B[i]*fxym[i+1];

}

void printout(int n);

int main()

{

int n,i; char ch;

do{ cout<<"输入x的最大下标:";

cin>>n;

for(i = 0; i <= n; i++)

{

cout<<"Please put in X"<

cin>>x[i];

cout<<"Please put in Y"<

cin>>y[i];

}