牛顿迭代法 matlab程序(解线性方程组)

作者：佚名来源：转载发布时间：2009-3-7 16:55:53

减小字体 增大字体

1．功能

本程序采用牛顿法，求实系数高次代数方程

f(x)=a0xn+a1xn-1+…+an-1x+an=０　(an≠０ ) 　　　　(1)

的在初始值x0附近的一个根。

2.使用说明

(1)函数语句

Y=NEWTON_1(A,N,X0,NN,EPS1)

调用M文件newton_1.m。

(2)参数说明

A n+1元素的一维实数组，输入参数，按升幂存放方程系数。

N　　 整变量，输入参数，方程阶数。

X0 实变量，输入参数，初始迭代值。

NN　 整变量，输入参数，允许的最大迭代次数。

EPS1　实变量，输入参数，控制根的精度。

3．方法简介

解非线性方程f(x)=0的牛顿法是把非线性方程线性化的一种近似方法。把f(x)在x0点附近展开成泰勒级数

f(x)=f(x0)+(x－x0)fˊ(x0)+(x－x0)2 +…

取其线性部分，作为非线性方程f(x)=0的近似方程，则有

f(x0)+fˊ(x0)(x－x0)=0

设fˊ(x0)≠0则其解为

x1=x0－f(x0)/fˊ(x0)

再把f(x)在x1附近展开成泰勒级数，也取其线性部分作f(x)=0的近似方程。若f(x1)≠0，则得

x2=x1－f(x1)/fˊ(x1)

这样，得到牛顿法的一个迭代序列

xn+1=xn－f(xn)/fˊ(xn)

4．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(ˊnn is fullˊ)；

return；

end

end

y=x(i)；

i

程序中调用的n_f.m和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

function y=n_df(a,n,x)

y=0.0；

for i=1:n

y=y+a(i)*(n+１-i)*xˆ(n-i)；

end

5．程序附注

(1)程序中调用n_f.m和n_df.m文件。n_f.m是待求根的实数代数方程的函数,n_df.m是方程一阶导数的函数。由使用者自己编写。

(2)牛顿迭代法的收敛速度：如果f(x)在零点附近存在连续的二阶微商，ξ是f(x)的一个重零点，且初始值x0充分接近于ξ，那么牛顿迭代是收敛的，其收敛速度是二阶的，即平方收敛速度。

6．例题

用牛顿法求下面方程的根

f(x)=x3+2x2+10x-20

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

7．运行结果

>>a=[1,2,10,-20] ；

>>n=3；

>>x0=1；

>>nn=1000；

>>eps1=1e-8；

>>y=newton_1(a,n,x0,nn,eps1)

y=

1.1373e+000

i=

6

function fp = newton_interpolation(x,y,p)

% Script for Newtons Interpolation.

% Muhammad Rafiullah Arain

% Mathematics & Basic Sciences Department

% NED University of Engineering & Technology - Karachi

% Pakistan.

% ---------

% x and y are two Row Matrices and p is point of interpolation

%

% Example

% >> x=[1,2,4,7,8]

% >> y=[-9,-41,-189,9,523]

% >> newton_interpolation(x, y, 5)

% OR

% >> a = newton_interpolation(x, y, 5)

n = length(x);

a(1) = y(1);

for k = 1 : n - 1

d(k, 1) = (y(k+1) - y(k))/(x(k+1) - x(k));

end

for j = 2 : n - 1

for k = 1 : n - j

d(k, j) = (d(k+1, j - 1) - d(k, j - 1))/(x(k+j) - x(k));

end

end

d

for j = 2 : n

a(j) = d(1, j-1);

end

Df(1) = 1;

c(1) = a(1);

for j = 2 : n

Df(j)=(p - x(j-1)) .* Df(j-1);

c(j) = a(j) .* Df(j);

end

fp=sum(c);

展开阅读全文