牛頓迭代法-matlab程序(解線性方程組)

1、牛頓迭代法 matlab程序（解線性方程組）作者：佚名  來(lái)源：轉(zhuǎn)載  發(fā)布時(shí)間：2009-3-7 16:55:53減小字體 增大字體 1功能本程序采用牛頓法，求實(shí)系數(shù)高次代數(shù)方程f(x)=a0xn+a1xn-1+an-1x+an=（an ） （1）的在初始值x0附近的一個(gè)根。2.使用說(shuō)明（1）函數(shù)語(yǔ)句Y=NEWTON_1(A,N,X0,NN,EPS1)調(diào)用M文件newton_1.m。（2）參數(shù)說(shuō)明A n+1元素的一維實(shí)數(shù)組，輸入?yún)?shù)，按升冪存放方程系數(shù)。N 整變量，輸入?yún)?shù)，方程階數(shù)。X0 實(shí)變量，輸入?yún)?shù)，初始迭代值。NN 整變量，輸入?yún)?shù)，允許的最

2、大迭代次數(shù)。EPS1實(shí)變量，輸入?yún)?shù)，控制根的精度。3方法簡(jiǎn)介解非線性方程f(x)=0的牛頓法是把非線性方程線性化的一種近似方法。把f(x)在x0點(diǎn)附近展開(kāi)成泰勒級(jí)數(shù)f(x)=f(x0)+(xx0)f(x0)+(xx0)2 +取其線性部分，作為非線性方程f(x)=0的近似方程，則有f(x0)+f(x0)(xx0)=0設(shè)f(x0)0則其解為x1=x0f(x0)/f(x0)再把f(x)在x1附近展開(kāi)成泰勒級(jí)數(shù)，也取其線性部分作f(x)=0的近似方程。若f(x1)0，則得x2=x1f(x1)/f(x1)這樣，得到牛頓法的一個(gè)迭代序列xn+1=xnf(xn)/f(xn)4newton_1.m程序fun

3、ction 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；endendy=x(i)；i程序中調(diào)用的n_f.m和n_df.m文件如下：function y=n_df(a,n,x)%方程一階導(dǎo)數(shù)的函數(shù)y=0.0;for i=1:ny=y+a(i)*(n+1-i)*x(n-i);endfunction y=n_df(a

4、,n,x)y=0.0；for i=1:ny=y+a(i)*(n+-i)*x(n-i)；end5程序附注（1）程序中調(diào)用n_f.m和n_df.m文件。n_f.m是待求根的實(shí)數(shù)代數(shù)方程的函數(shù),n_df.m是方程一階導(dǎo)數(shù)的函數(shù)。由使用者自己編寫。（2）牛頓迭代法的收斂速度：如果f(x)在零點(diǎn)附近存在連續(xù)的二階微商，是f(x)的一個(gè)重零點(diǎn)，且初始值x0充分接近于，那么牛頓迭代是收斂的，其收斂速度是二階的，即平方收斂速度。6例題用牛頓法求下面方程的根f(x)=x3+2x2+10x-20 y=y+a(i)*(n+1-i)*x(n-i);7運(yùn)行結(jié)果>>a=1,2,10,-20 ；>>

5、n=3；>>x0=1；>>nn=1000；>>eps1=1e-8；>>y=newton_1(a,n,x0,nn,eps1)y=1.368808107821373e+000i=6function fp = newton_interpolation(x,y,p)% Script for Newton's Interpolation.% Muhammad Rafiullah Arain% Mathematics & Basic Sciences Department% NED University of Engineering &

6、; 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);endfor 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); endenddfor j = 2 : n a(j) = d(1, j-1);