南京邮电大学 数值计算实践报告

数值计算实践

I 、方程求根

一、实验目的

熟悉和掌握Newton 法,割线法,抛物线法的方法思路，并能够在matlab 上编程实现

二、问题描述

(1).给定一个三次方程,分别用Newton 法,割线法,抛物线法求解. 方程的构造方法:

(a)根:方程的根为学号的后三位乘以倒数第二位加1再除以1000. 假设你的学号为B06060141,则根为141*(4+1)/1000=0.564

(b)方程:以你的学号的后三位数分别作为方程的三次项,二次项,一次项的系数,根据所给的根以及三个系数确定常数项. 例如:

你的学号是B06060141,则你的方程是x 3+4x 2+x+a 0=0的形式. 方程的根为0.564,因此有

0.5643+4*0.5642+0.564+a0=0,于是a0=-2.015790144 你的方程为x 3+4x 2+x-2.015790144=0.

(2)假设方程是sinx+4x 2+x+a0=0的形式（三个系数分别是学号中的数字），重新解决类似的问题

(3)构造一个五次方程完成上面的工作.

四次方程的构造:将三次多项式再乘以(x-p*)2得到对应的五次多项式(p*为已经确定的方程的根,显然,得到的五次方程有重根).

(4)将（2）中的方程同样乘以(x-p*)得到一个新的方程来求解

注:(1)Newton 法取0.5为初值,割线法以 0,1为初值,抛物线法以0,0.5,1为初值, (2)计算精度尽量地取高.

终止准则：根据ε<--||1n n p p 来终止

(3)可供研究的问题：

（一）ε的取值不同对收敛速度有多大的影响

（二）将注（1）中的初值该为其它的初值，对收敛性以及收敛速度有无影响 （三）能否求出方程的所有的根 （4）实验报告的撰写

实验报告包含的内容：（一）实验目的（二）问题描述（三）算法介绍（包括基本原理）（四）程序（五）计算结果（六）结果分析（七）心得体会

三、算法介绍

在本问题中，我们用到了newton 法，割线法，抛物线法。 1.Newton 法迭代格式为： )

()

('1k k k k x f x f x x -

=+ 当初值与真解足够靠近，newton 迭代法收敛，对于单根，newton 收敛速度很快，对于重根，收敛较慢。

2.割线法：为了回避导数值的计算，使用上的差商代替，得到割线法迭代公式：

)

()()

)((111--+---

=k k k k k k k x f x f x x x f x x

割线法的收敛阶虽然低于newton 法，但迭代以此只需计算一次函数值，不需计算其导数，所以效率高，实际问题中经常应用。

3.抛物线法：可以通过三点做一条抛物线，产生迭代序列的方法称为抛物线法。其迭代公式为：

)

)(](,,[)](,[][)(1112-----+-+=-k k k k k k k k x x x x x x x f x x x x f x f x p k 其中 ],,[],,[211---k k k k k x x x f x x f 是一阶均差和二阶均差。 收敛速度比割线法更接近于newton 法。

对于本问题的解决就以上述理论为依据。终止准则为：ε<--||1n n x x 本题中所有精度取1e-8。

四、程序计算结果 问题一

根据所给的要求，可知待求的方程为：03356.0223=-++x x x 牛顿法

建立newton_1.m 的源程序，源程序代码为：

function y=newton_1(a,n,x0,nn,eps1) x(1)=x0; b=1; i=1;

while(abs(b)>eps1*x(i))

i=i+1;

x(i)=x(i-1)-n_f(a,n,x(i-1))/n_df(a,n,x(i-1));

b=x(i)-x(i-1);

if(i>nn)error

return;

end

end

y=x(i);

建立n_f.m的源程序来求待求根的实数代数方程的函数，源程序代码为：function y=n_f(a,n,x)% 待求根的实数代数方程的函数

y=0.0;

for i=1:(n+1)

y=y+a(i)*x^(n+1-i);

end

建立n_df.m的源程序来方程的一阶导数的函数，源程序代码为：function y=n_df(a,n,x)% 方程的一阶导数的函数

y=0.0;

for i=1:n

y=y+a(i)*(n+1-i)*x^(n-i);

end

在matlab软件中执行下列语句并得到的最终结果截图：

割线法

建立gexian.m的源程序，源程序代码为

function x=gexian(f,x0,x1,e)

if nargin<4,e=1e-4;end

y=x0;x=x1;i=0;

while abs(x-y)>e

i=i+1;

z=x-(feval(f,x)*(x-y))/(feval(f,x)-feval(f,y));

y=x;

x=z;

end

i

在matlab软件中执行下列语句并得到的最终结果截图：

抛物线

建立paowuxian.m的源程序，源程序代码为：

function x=paowuxian(f,x0,x1,x2,e)

if nargin<4,e=1e-4;end

x=x2;y=x1;z=x0;i=0;

while abs(x-y)>e

i=i+1;

h1=y-z;

h2=x-y;

c1=(feval(f,y)-feval(f,z))/h1;

c2=(feval(f,x)-feval(f,y))/h2;

d=(c1-c2)/(h1+h2);

w=c2+h2*d;

xi=x-(2*feval(f,x))/(w+(w/abs(w))*sqrt(w^2-4*feval(f,x)*d));

z=y;

y=x;

x=xi;

end

i

在matlab软件中执行下列语句并得到的最终结果截图：

研究一：只改变初值

由上述结果可知，方程的解在0.2附近，所以将牛顿法为0.2；割线法的初值设为0，0.4；抛物线法的初值设为0，0.2，0.4；

牛顿法

根据问题1中牛顿法的程序，在matlab软件中执行下列语句并得到的最终结果截图：

割线法

根据问题1中割线程序，在matlab软件中执行下列语句并得到的最终结果截图：

抛物线法

根据问题1中抛物线法程序，在matlab软件中执行下列语句并得到的最终结果截图：

研究二只改变精度

将精度由1e-8改为1e-50和1e-100观察迭代次数有何变化

牛顿法：

根据问题1中的牛顿法的程序，在matlab软件中执行下列语句并得到的最终结果截图：

精度为1e-50时

精度为1e-100时

割线法

根据问题1中的割线法的程序，在matlab软件中执行下列语句并得到的最终结果截图：

精度为1e-50时

精度为1e-100时

抛物线法

根据问题1中的抛物线法的程序，在matlab软件中执行下列语句并得到的最终结果截图：

精度为1e-50时

精度为1e-100时

、

研究结论

在只改变初值时，当初值定得越靠近初值，迭代次数就越少。

在只改变精度时，当精度越来越大时，迭代次数并几乎不变。

综上所述，初值对迭代次数的影响比较大，精度对迭代次数影响不大。

问题二

问题描述

根据所给的要求，可知待求的方程为：03356.02sin 2=-++x x x

问题分析

仍然利用（1）中方法求解这一问题，并利用图解法找到初值，通过观察图像，将newton 法初值设为：0.1，割线法初值设为：0,0.2。抛物线法初值设为：0,0.1,0.2。 图像见下图：

00.10.20.30.40.50.60.70.80.91

-0.5

00.511.522.533.5

4

Newton 法

问题一的牛顿法的求解只适用于线性方程，所以在问题二中用其他方法来求解方程。

建立newton1.m 源程序，源程序代码为： function x=newton1(fn,dfn,x0,e) if nargin<4,e=1e-4;end

x=x0;x0=x+2*e;i=1；

while abs(x0-x)>e

i=i+1；

x0=x;

x=x0-feval(fn,x0)/feval(dfn,x0);

end

在matlab软件中执行下列语句并得到最终结果截图

割线法

利用问题一中的割线法程序，在matlab软件中执行下列语句

抛物线法

利用问题一中的抛物线法程序，在matlab 软件中执行下列语句

问题三

问题描述

按照题目要求对五次方程进行构造为：0)224.0)(3356.02(223=--++x x x x

问题分析

仍然利用一中方法求解这一问题，并利用图解法找到初值，通过观察图像，将牛顿法的两组初值为0以及0.5，割线法初值设为：0,0.5以及0，0.2。抛物线法初值设为：两组初值为0,0.5,1以及0,0.25,0.5。

00.050.10.150.20.250.3

-18

-16-14-12-10-8-6-4-202

-3

牛顿法

利用问题二中的程序，在matlab 软件中执行下列语句： 初值为0时

初值为0.5时

割线法

利用问题一中割线法的程序，在matlab软件中执行下列语句：初值为0，0.5时

初值为0，0.2时

抛物线法

利用问题一中抛物线法的程序，在matlab软件中执行下列语句：初值为0，0.5，1时

初值为0，0.25，0.5时

问题四

问题描述

根据题目要求对方程进行构造为：0)224.0)(3356.02(sin 2=--++x x x x

问题分析

仍然利用问题一中方法求解这一问题，并利用图解法找到初值，通过观察图像，newton 法初值选取了两组初值为0以及0.5，割线法初值设为：0，0.5和0，0.3。抛物线法初值设为：两组初值为0,0.2,0.4以及0，0.5，1。

0.050.10.150.20.250.3

-0.01

00.010.020.030.040.050.060.070.08

牛顿法

利用问题二中的程序，在matlab 软件中执行下列语句： 初值为0时

初值为0.5时

割线法

利用问题一中割线法的程序，在matlab软件中执行下列语句：初值为0，0.5时

初值为0，0.3时

抛物线法

利用问题一中抛物线法的程序，在matlab软件中执行下列语句：

初值为0，0.1，0.2时

初值为0，0.5，1时

五、计算结果及分析

Newton法割线法抛物线法问题一0.224 0.224 0.224

问题二0.1466 0.1466 0.1466 问题三0.224 0.224 0.224

问题四0.224和

0.146577258557101 0.224和

0.1466

0.224和

0.1466

结果分析

将Newton法，割线法，抛物线法进行比较可以看到在本文题中，三种方法计算得到的最终结果基本相同，但是迭代步数有较大差别，综合看来Newton法迭代步数最少，割线法次之，抛物线法最次。在各个问题的研究中，我通常都会采用不同的初值，发现不同初值会对应不同的迭代次数，另外针对问题一，我选用了不同的精度，发现迭代次数并没有很大的变化，因而一个初值的选定可以对迭代次数产生很大的影响，而精度的改变对迭代次数的影响很小。

对每个算法单独来看，显然选择初值不同对于迭代步数影响较大，对于找到根也会有影响。因此应该先通过画图确定根的大致位置，给出在其附近的初值。

六、心得体会

在实现这三个算法的过程中，本身编程较易实现，最重要的是对算法本身的理解，只有真正理解算法的含义才能更快更好的实现程序。

II、离散型最小二乘和连续型最小二乘问题一、实验目的

掌握并能够利用离散型最小二乘和连续性最小二乘求解问题。

二、问题描述

1：以函数错误!未找到引用源。生成的6个数据点（0.25,23.1）,（1.0,1.68）, （1.5,1.0）,（2.0,0.84）,（2.4,0.826）,（5.0,1.2576）为例，运行程序得到不同阶对应的曲线拟合产生的多项式函数。

2. 例题：计算

f(x)=exp(x)在[-1，1]上的二、三次最佳平方逼近多项式。

并画图进行比较。

三、问题分析

问题1是离散最小二乘问题。

最小离散最小二乘就是根据一批有误的数据（错误!未找到引用源。）

i=1,2，…，n 确定参数，并通过均方误差来比较曲线拟合的优劣，在本题中通过画图来比较不同阶方程拟合效果的优劣。 选择两种方法实现离散最小二乘。

方法一，建立normal equation(法方程组)，求解k 次多项式系数。法方程组构造方法：

∑∑∑∑∑∑∑∑∑∑∑∑===+===+=======+++=+++=+++m

i n

i

i m

i n

i

n m

i n i

m

i n

i m

i i

i m

i n i

n m

i i m

i i m

i i

i m i n

i n m i i m i i x y x a x a x a x y x a x a x a x y x a x a x a 0

20

1

10

001

1

2

10

1

00

001

100

0.... . . . . . ......

方法二：由于在matlab 中存在ployfit 函数，可以即为方便的用k 次多项式拟合。

问题2是一个连续型最小二乘法的应用实例。

对于给定的函数],[)(b a C x f ∈，如果存在

*01(){(),(),,()}

n S x Span x x x ???∈L

使得

[]22

*

()()()min ()()()b

b

a

a

a x b

x f x S x dx x f x s x dx

ρρ≤≤??-=-???

?

则称S *(x )是f (x )在集合01{(),(),,()}n Span x x x ???L 中的最佳平方逼近函数。

显然，求最佳平方逼近函数)()(0**

x a x S j n

j j ??=∑=的问题可归结为求它的系数

*

*1*0,,,n

a a a Λ，使多元函数 dx x a x f x a a a I j n j j b

a

n 2

010)()()(),,,(??

?

???-=∑?

=?ρΛ

取得极小值，也即点(*

*1*0,,,n

a a a Λ)是I (a 0, …，a n )的极点。由于I (a 0, a 1, …，a n )是关于a 0, a 1, …，a n 的二次函数，利用多元函数取得极值的必要条件，