牛顿科茨matlab,牛顿-科茨（Newton-Cotes）公式算法.doc

牛顿-科茨(Newton-Cotes)公式算法.doc

牛顿－科茨(Newton-Cotes)公式算法：

Step 1：判断type类型，1转Step 2,2转Step 3,3转Step 4；否则输出值为0；

Step 2：计算科茨公式；

Step 3：计算牛顿-科茨六点公式；

Step 4：计算牛顿-科茨七点公式；

流程图：

开始type=2type=1计算科茨六点公式计算科茨公式type=3计算科茨七点公式输入type

开始

type=2

type=1

计算科茨六点公式

计算科茨公式

type=3

计算科茨七点公式

输入type

N N N

输出结果为0输出结果输出结果输出结果

输出结果为0

输出结果

输出结果

输出结果

定义function函数如下：

function I = NewtonCotes(f,a,b,type)

%type = 1 科茨公式

%type = 2 牛顿－科茨六点公式

%type = 3 牛顿－科茨七点公式

I=0;

switch type

case 1,

I=((b-a)/90)*(7*subs(sym(f),findsym(sym(f)),a)+...

32*subs(sym(f),findsym(sym(f)),(3*a+b)/4)+...

12*subs(sym(f),findsym(sym(f)),(a+b)/2)+...

32*subs(sym(f),findsym(sym(f)),(a+3*b)/4)+7*subs(sym(f),findsym(sym(f)),b));

case 2,

I=((b-a)/288)*(19*subs(sym(f),findsym(sym(f)),a)+...

75*subs(sym(f),findsym(sym(f)),(4*a+b)/5)+...

50*subs(sym(f),findsym(sym(f)),(3*a+2*b)/5)+...

50*subs(sym(f),findsym(sym(f)),(2*a+3*b)/5)+...

75*subs(sym(f),findsym(sym(f)),(a+4*b)/5)+19*subs(sym(f),findsym(sym(f)),b));

case 3,

I=((b-a)/840)*(41*subs(sym(f),findsym(sym(f)),a)+...

216*subs(sym(f),findsym(sym(f)),(5*a+b)/6)+...

27*subs(sym(f),findsym(sym(f)),(2*a+b)/3)+...

272*subs(sym(f),findsym(sym(f)),(a+b)/2)+...

27*subs(sym(f),findsym(sym(f)),(a+2*b)/3)+...

216*subs(sym(f),findsym(sym(f)),(a+5*b)/6)+41*subs(sym(f),findsym(sym(f)),b));

end

运算的输入方式及结果为：

>> y=NewtonCotes('sin(x)',0,10,1)

y =

3.7613

>> y=NewtonCotes('sin(x)',0,10,2)

y =

2.7865

>> y=NewtonCotes('sin(x)',0,10,3)

y =

1.5296

２、编程解决以下科学计算问题。

1.编程解决以下科学计算问题。

求编写求解方程函数文件，然后调用该函数文件求30000个点上的x，y，最后在所有的(x、y)坐标处标记一个点(不要连线)绘出图形。这种图形又称埃农(Henon)引力线图，他将迭代迟来的随机点吸引到一起，最后得出貌似连贯的引力线图。

流程图：

开始

开始

创建初始矩阵，并赋初值(0，0)，固定点n＝30000

创建初始矩阵，并赋初值(0，0)，固定点n＝30000

将第K行的值通过迭代关系赋值给第(K＋1)行变量为k，k<3000

将第K行的值通过迭代关系赋值给第(K＋1)行

变量为k，k<3000

N