数值课设

沈阳航空航天大学

数值分析课程设计

（设计程序）

题目 编写复合辛普森求积方法的M文件

班级 / 学号 143414010106 学 生 姓 名 林琳 指 导 教 师 刘勇进

沈阳航空航天大学

课 程 设 计 任 务 书

课 程 名 称 数值分析课程设计 院（系） 理学院 专业 信息与计算科学 班级 1434140101 学号 143414010106 姓名 林琳 课程设计题目 编写复合辛普森求积方法的M文件 课程设计时间: 2015 年 12 月 21 日至 2016 年 1 月 1日

课程设计的要求及内容：

[要求]

1. 学习态度要认真，要积极参与课程设计，锻炼独立思考能力； 2. 严格遵守上机时间安排；

3. 按照MATLAB编程训练的任务要求来编写程序； 4. 根据任务书来完成课程设计论文；

5. 报告书写格式要求按照沈阳航空航天大学“课程设计报告撰写规范”； 6. 报告上交时间：课程设计结束时上交报告； 7. 严禁抄袭行为，一旦发现，课程设计成绩为不及格。

一、 通识题

1. 作出下列函数的图形：

(1) 编写M文件对下面函数进行分段曲面绘制：

(2)

及

；

.

要求：a/ 分别考虑

b/

2. 用数值方法计算

其中n分别取100, 200, 1000。

3. 求由参数方程4. 作出函数求由方程5. 二、

，

确定的函数的二阶导数所确定的隐函数的导数

的曲线在区间

上的两条水平切线。。

。 。

编写复合辛普森求积方法的M文件

要求： 1） 了解复合辛普森求积方法的基本原理；

2） 绘出程序流程图；

3）编写复合辛普森求积方法的M文件； 4）利用编写M文件，计算?sin?x?dx；

0? 5）给出程序注释。

指导教师 年 月 日 负责教师 年 月 日 学生签字 年 月 日

沈阳航空航天大学

课 程 设 计 成 绩 评 定 单

课 程 名 称 数值分析课程设计 院（系） 理学院 专业 信息与计算科学 课程设计题目 编写复合辛普森求积方法的M文件 学号 143414010106 姓名 林琳 指导教师评语：

课程设计成绩

指导教师签字 年 月 日

数值分析课程设计 沈阳航空航天大学 课程设计用纸 目 录

目 录

一、 正文 ............................................................ 1

1通识题

1.1作出下列函数的图形………………………………………………………………………...1

1.2用数值方法计算

1.3求由参数方程

1.4作出函数求由方程1.5

，其中n分别取100, 200, 1000。…………………..4

确定的函数的二阶导数所确定的隐函数的导数

的曲线在区间

。…...4 。…....4

上的两条水平切线。………………..4

2编写复合辛普森求积方法的M文件……………………………………………………………5

要求： 1） 了解复合辛普森求积方法的基本原理；…………………………5 2） 绘出程序流程图；…………………………………………………7 3）编写复合辛普森求积方法的M文件

4）利用编写M文件，计算?sin?x?dx；……………………………..7

0? 5）给出程序注释..…………………………………………………….13

二、总结 ............................................................... 9 参考文献 .............................................................. 10 附 录 ????????????????????????????????11

第 I 页

数值分析课程设计 沈阳航空航天大学 课程设计用纸 正 文

一、 正文

一1.（1）编写M文件对下面函数进行分段曲面绘制：

这道题让我们画的是函数的三维立曲面，有matlab,我们知道

三维画图指令是surf(),suef()是用来绘制三维曲面图指令在 这道题里z是平方的形式，画图要开平方，在这里我们用指令 real(),real()是取实部。源程序见附录1.1.1

(2)

要求：a/ 分别考虑

b/

.

；

及

分析：在这里要注意反正切函数在这里输入为atan(x)，另外还用到了二维画图指令

第 1 页

数值分析课程设计 沈阳航空航天大学 课程设计用纸 正 文

plot.

源程序见附录1.1.2.1 函数图像为：

21.510.50-0.5-1-1.5-2-50-40-30-20-1001020304050

这里的关键是将函数表示出来，程序中我们用for循环 实现函数的输入，plot()实现画图. 源程序见附录1.1.2.2 函数图像如图： n=5,

n=10时

第 2 页

数值分析课程设计 沈阳航空航天大学 课程设计用纸 正 文

n=20时

n=50时

第 3 页

数值分析课程设计 沈阳航空航天大学 课程设计用纸 正 文

2.用数值方法计算

，

其中n分别取100, 200, 1000。

分析：此函数由n项累和组成，在程序中我们用for函数实现s函数，定义全局变量（c）在程序中我们取c=9[100,200,1000],求得n分别为100,200,1000时分别对应的函数值向量s

源程序见附录1.1.2. n=100 1.6350

n=200 1.6399

n=500 1.6439 3.求由参数方程

确定的函数的二阶导数

。

分析：在这道题中我们要用到求导指令diff（y,x）,和定义符号变量的指令sysm(). 源程序见附录1.1.3. 执行过程

f2 =

- (6*t + 2)/(1/(t + 1) - 1) - (3*t^2 + 2*t)/((1/(t + 1) - 1)^2*(t + 1)^2) 4.作出函数求由方程

所确定的隐函数的导数

。

分析：本题是求二阶导数，x,y是关于t的参数方程，求二阶导数我们在这里用两次diff（）指令.

源程序见附录1.1.4 运行过程和结果：

x*exp(y/x)*(1/x + (y*exp(-y/x))/x^2) 5.

的曲线在区间

上的两条水平切线。。

分析：这是一个作图题在图上我们要画函数f(x)的图像，还要求出函数的切线，做

出切线的图像，我们用导数等于零求出切点的横坐标，带入原方程求出切线.程序中我们主要用到了求导指令diff(x),作图指定plot(). 源程序见附录1.1.5. 运行过程与结果：

第 4 页

数值分析课程设计 沈阳航空航天大学 课程设计用纸 正 文

二.编写复合辛普森求积方法的M文件

要求： 1） 了解复合辛普森求积方法的基本原理；

2） 绘出程序流程图；

3）编写复合辛普森求积方法的M文件； 4）利用编写M文件，计算?sin?x?dx；

0? 5）给出程序注释。

1·复合辛普森求积公式原理 辛普森求积公式

设将积分区间[a,b]划分为n等份，步长h=(b-a)/n,选取等距节点xk?a?kh构造出的插值型求积公式In?(b?a)?Ck(n)f(xk) 称为Newton-Cotes公式，式中Ck称为柯

(n)k?0n特斯系数。按Ak?x?a?th，则有

?balk(x)d(x),k?0,1,2,...,n，引进变量

第 5 页

数值分析课程设计 沈阳航空航天大学 课程设计用纸 正 文

C(n)knnhnnt?j(?1)n?k?(t?j)dt?dt?nk!(n?k)!?0?....................(2) b?a?0j?0k?jj?0j?kj?k由于是多项式的积分，柯特斯系数的计算不会遇到实质性的困难。 当n=1时，C(1)0?C(1)1?1，这时的求积公式就是我们所熟悉的梯形公式。 2当n=2时，按（2）式，这时的柯特斯系数为

(2)C0?C1(2)(2)C2121(t?1)(t?2)dt?4?06124???t(t?2)dt?206 121??t(t?1)dt?4062b?aa?b[f(a)?4f()?f(b)]就为辛普森求积公式。?f(x)k

162 将区间[a,b]划分为n等份，分点xk=a+kh,h=(b-a)/n,k=0,1,...,n,在每个子区间

相应的求积公式S?[xk,xk+1](k=0,1,?,n-1)上采用辛普森公式，若记

则有

记

=为复合辛普森公式，

余项为,。

于是当f（x）

，与复合辛普森公式余项=,

第 6 页

数值分析课程设计 沈阳航空航天大学 课程设计用纸 正 文

2.程序流程图

开始 输入n,a,b,g

定义未知变量x，将g定义为函数f(x)，计算h的值 调用函数利用for循环计算和号的值 利用复合辛普森公式，带入数值，计算S的值

输出S值

结束

3.计算?sin?x?dx

0?> Untitled5(2,0,pi,@(x)sin(x)) f =

@(x)g(x) ans =

2.0046

>> Untitled5(3,0,pi,@(x)sin(x)) f =

@(x)g(x)

第 7 页

数值分析课程设计 沈阳航空航天大学 课程设计用纸 正 文

ans =

2.0009

>> Untitled5(5,0,pi,@(x)sin(x)) f =

@(x)g(x) ans =

2.0001

>> Untitled5(7,0,pi,@(x)sin(x)) f =

@(x)g(x) ans =

2.0000

>> Untitled5(10,0,pi,@(x)sin(x)) f =

@(x)g(x) ans =

2.0000 >>

第 8 页

数值分析课程设计 沈阳航空航天大学 课程设计用纸 正 文

二、总结

通过这两周的MATLAB课程设计，我学到了很多知识，也了解到MATLAB是

一个运用于很多领域的软件。初次接触MATLAB，感到自己还有很多不足，对一些程序的设计还没有足够精通，但是通过两周的学习，在MATLAB设计方面，觉得自己还是有所进步。在课设过程中，也遇到过很多问题，有时候程序运行不出来或者结果错误，仅仅只是因为忘记加分号，或者是多加了括号，或者是字母的大小写不一致等，这让我意识到MATLAB的严谨性，也在此方面更加注意，避免此类错误发生。有时是因为函数编写思路不正确，通过上网查找，去图书馆查找，以及请教同学等方法来解决这些问题。MATLAB以矩阵计算为基础，在科学研究以及工程设计邻域有着十分广泛的应用，我以后会继

续学习，掌握更多的知识，更好的运用MATLAB。

完稿日期:2016 年1 月 4日

第 9 页

数值分析课程设计 沈阳航空航天大学 课程设计用纸 参考文献

参考文献

[1] 李庆阳,王能超,易大义.数值分析.5版.北京:清华大学出版社,2008.12

[2] 张霞萍,MATLAB 8.X程序设计及典型应用.西安:西安电子科技大学出版社, 2014.7

第 10 页

数值分析课程设计 沈阳航空航天大学 课程设计用纸 附 录

附 录

源程序: （1.1.1）

x=-2:0.02:2;%从-2，增量为0.02 y=x; %从-2，增量为0.02

[X,Y]=meshgrid(x,y);%对X.Y取x,y的网格坐标矩阵 Z=zeros(201);%取Z为201*1个的零矩阵 for i=1:201 %利用for循环从1到201 for j=1:201%利用for循环从1到201 if X(i,j)+Y(i,j)>1

Z(i,j)=0.5*exp(-0.75*Y(i,j).^2-3.75*X(i,j).^2-1.5*X(i,j)); %如果在X+Y>1范围内将X(i,j), Y(i,j)利用相对应的表达式带入Z零矩阵相对应的点中

else if X(i,j)+Y(i,j)<=-1

Z(i,j)=0.5*exp(-0.75*Y(i,j).^2-3.75*X(i,j).^2+1.5*X(i,j)); %否则，如果在X+Y<=1范围内将X(i,j), Y(i,j) 利用相对应的表达式带入Z零矩阵相对应的点中

else Z(i,j)=0.8*exp(-Y(i,j).^2-6*X(i,j).^2);

%否则将X(i,j), Y(i,j) 利用相对应的表达式带入Z零矩阵其余的点中

end end

end %利用条件语句写出题中给的分段函数

mesh(X,Y,Z);%作出三维曲面图 end

（1.1.2.1）

x=-50:1:50; %从-50到50，增量为1 y=atan(x); %将带入表达式中 plot(x,y); %作出图像

（1.1.2.2）

function y=myfun(n) %建立变量为n的函数y x=-50:1:50 %从-50到50，增量为1 for k=1:n %用for循环从1到n p=0;

p=p+(-1)^(k-1)/(2*k-1).*x.^(2*k-1); %带入表达式

第 11 页

数值分析课程设计 沈阳航空航天大学 课程设计用纸 附 录

end

plot(x,p); %作出图像

（1.2）

function s=l(n) %建立变量为n的函数s s=0;

for k=1:n %对k从1一直取n到 s=s+1./k^2; %带入表达式 end s （1.3）

syms t %定义符号变量t

y=t^3+t^2; %带入t对y的表达式 x=t-log(1+t); %带入t对x的表达式 f1=diff(y)/diff(x); %求y对x的一阶导数. f2=diff(f1)/diff(x); %求y对x的二阶导数. f2

(1.4)

syms x y %定义符号变量x y f=log(x)+exp(-y/x); %写出f(x)表达式 diff(f,'x'); %求f对x的偏导数. diff(f,'y'); %求f对y的偏导数. z=(-1)* diff(f,'x')/diff(f,'y') %求y对x的导数

（1.5）：

syms x ; %定义符号变量x. f=4*(x.^3)-6*(x.^2)+x-2; df=diff(f,x);

x0=solve(df,x); %解函数df的零点. y1=4*(x0(1).^3)-6*(x0(1).^2)+x0(1)-2;

y2=4*(x0(2).^3)-6*(x0(2).^2)+x0(2)-2;%将零点分别带入表达式其结果为y1,y2 x=0:0.002:1;

b=4*(x.^3)-6*(x.^2)+x-2; plot(x,b)

hold on %作出原函数图像

plot(x,y1,'r') %将x0(1),y1带入切线表达式中 hold on %作出一条切线图像

plot(x,y2,'r') %将x0(2),y2带入切线表达式中 hold on %作出另一条切线图像 二

第 12 页

数值分析课程设计 沈阳航空航天大学 课程设计用纸 附 录

function ans=Untitled5(n,a,b,g)

%建立变量为n,a,b,g的函数ans; n为被区间被分个数；a,b为被积区间；g为被积函数 f=@(x)g(x) %定义句柄g为函数f(x) h=((b-a)/n); %计算h

m=0; %将0赋值给m

for i=1:n-1 %用for循环从1到n-1 k1=a+i*h; %将a+i*h赋值给k1

m=m+feval(f,k1); %调用函数f(x)在k1点的值，利用for循环计算m值 end

t=0; %将0赋值给t

for j=0:n-1 %用for循环从0到n-1

k2=a+(1/2+j)*h; %将a+(1/2+j)*h赋值给k2

t=t+feval(f,k2); %调用函数f(x)在k2点的值，利用for循环计算t值 end

p=feval(f,a)+feval(f,b); %调用函数计算f(a),f(b)的值 ans=h/6*(p+4*t+2*m); %带入公式计算积分值

第 13 页

本文来源：https://www.bwwdw.com/article/0s76.html