MATLAB基礎(chǔ)及應(yīng)用曹敦虔編著第5章數(shù)值計(jì)算本章學(xué)習(xí)重點(diǎn)掌握內(nèi)容：插值與擬合方程求根極值問題數(shù)值解法數(shù)值積分與微分方程數(shù)值解第5章數(shù)值計(jì)算5.1 插值與擬合5.2 非線性方程求根5.3 迭代法求線性方程組數(shù)值解5.4 無約束一維極值問題5.5 無約束多維極值5.6 數(shù)值積分5.7 微分方程數(shù)值解第5章數(shù)值計(jì)算數(shù)值計(jì)算是指使用數(shù)字計(jì)算機(jī)求數(shù)學(xué)問題近似解的方法與過程,數(shù)值運(yùn)算的研究領(lǐng)域包括方程數(shù)值解、數(shù)值逼近、最優(yōu)化方法、數(shù)值微分和數(shù)值積分、常微分方程數(shù)值解法、數(shù)值代數(shù)、積分方程數(shù)值解法、偏微分方程數(shù)值解法、計(jì)算幾何、計(jì)算概率統(tǒng)計(jì)等.隨著計(jì)算機(jī)的廣泛應(yīng)用和發(fā)展,許多計(jì)算領(lǐng)域的問題,如計(jì)算物理、計(jì)算力學(xué)、計(jì)算化學(xué)、計(jì)算經(jīng)濟(jì)學(xué)等都可歸結(jié)為數(shù)值計(jì)算問題.本章著重介紹MATLAB在解決數(shù)值逼近、方程求根、線性方程組數(shù)值解、最優(yōu)化方法、數(shù)值積分、微分方程數(shù)值解等問題的思路和方法.5.1插值與擬合插值和擬合是兩種重要的數(shù)據(jù)逼近方法.插值是尋求一個(gè)函數(shù),使函數(shù)在某些已知點(diǎn)的值等于已知的數(shù)值,即函數(shù)圖像經(jīng)過已知的數(shù)據(jù)點(diǎn).擬合是尋求一個(gè)函數(shù),使函數(shù)在某些已知點(diǎn)的值盡量接近已知的數(shù)值.兩種方法適用范圍不同,插值是假定已知的數(shù)據(jù)是精確的,所以要求函數(shù)經(jīng)過數(shù)據(jù)點(diǎn),而擬合是假定已知的數(shù)據(jù)存在誤差,所以只需接近數(shù)據(jù)點(diǎn)即可.插值擬合5.1.1多項(xiàng)式插值

(5.2)插值5.1.1多項(xiàng)式插值

例5.1使用多項(xiàng)式對(duì)下面的數(shù)據(jù)進(jìn)行插值。由于總共有5個(gè)點(diǎn),所以使用4次多項(xiàng)式進(jìn)行插值.x136913y25374%插值基函數(shù)[x^n,x^(n-1),...,1]f=@(x)repmat(x,1,n+1).^repmat(n:-1:0,n+1,1);A=f(x);%生成系數(shù)矩陣c=A\y;%求解線性方程組xx=1:0.1:13;yy=polyval(c,xx);%求多項(xiàng)式在xx處的值plot(x,y,'o',xx,yy);5.1.2分段多項(xiàng)式插值當(dāng)數(shù)據(jù)點(diǎn)較多時(shí),整體的多項(xiàng)式插值結(jié)果是一個(gè)高次多項(xiàng)式,容易出現(xiàn)“龍格現(xiàn)象”,導(dǎo)致插值函數(shù)在插值節(jié)點(diǎn)之間發(fā)生劇烈的波動(dòng),造成較大的誤差.為避免出現(xiàn)“龍格現(xiàn)象”,可使用次數(shù)較低的多項(xiàng)式,對(duì)數(shù)據(jù)點(diǎn)分段進(jìn)行插值,在相鄰兩段之間滿足一定的連續(xù)性、光滑性條件,這種方法稱為分段多項(xiàng)式插值,常用的方法有分段線性插值、樣條函數(shù)插值.龍格現(xiàn)象5.1.2分段多項(xiàng)式插值

5.1.2分段多項(xiàng)式插值下面函數(shù)用于實(shí)現(xiàn)分段線性插值:functionyy=linepiece(x,y,xx)%分段線性插值函數(shù)%x,y插值數(shù)據(jù)%xx插值函數(shù)自變量%yy插值函數(shù)因變量yy=zeros(size(xx));fori=1:length(x)-1yy=yy+(xx>=x(i)&xx<x(i+1)).*((xx-x(i+1))./(x(i)-x(i+1))*y(i)+(xx-x(i))./(x(i+1)-x(i))*y(i+1));endyy=yy+(xx==x(end)).*((xx-x(end))./(x(end-1)-x(end))*y(end-1)+(xx-x(end-1))./(x(end)-x(end-1))*y(end));例5.2對(duì)例5.1中的數(shù)據(jù)進(jìn)行分段線性插值.只需要調(diào)用linepiece函數(shù)即可完成分段線性插值.x=[136913]';y=[25374]';xx=1:0.1:13;yy=linepiece(x,y,xx);%調(diào)用分段線性插值函數(shù)plot(x,y,'o',xx,yy);分段線性插值5.1.2分段多項(xiàng)式插值2.樣條函數(shù)插值樣條函數(shù)是指分段的、具有一定整體光滑性的一類函數(shù).一般情況下是指多項(xiàng)式樣條函數(shù),即每一段函數(shù)都是的一個(gè)多項(xiàng)式,在相鄰兩段之間連接處要求連續(xù)且光滑.樣條函數(shù)5.1.2分段多項(xiàng)式插值

(5.8)5.1.2分段多項(xiàng)式插值

(5.9)5.1.2分段多項(xiàng)式插值

(5.10)

x136913y25374

T這樣就得到了一個(gè)16元線性方程組，寫成矩陣形式為Ac=b

，其中求解得插值結(jié)果函數(shù)圖像如下圖所示.插值函數(shù)f(x)的圖像5.1.2分段多項(xiàng)式插值

5.1.2分段多項(xiàng)式插值

(5.12)5.1.2分段多項(xiàng)式插值

5.1.2分段多項(xiàng)式插值利用deBoor-Cox遞推公式計(jì)算出來的基函數(shù)稱為B樣條基函數(shù).需要注意的是,雖然在定義函數(shù)序列Bm時(shí)需要對(duì)區(qū)間進(jìn)行擴(kuò)展,但最終基函數(shù)的定義域仍取[a,b].例如,對(duì)于給定區(qū)間[0,1]上的剖分?n={0,0.25,0.5,0.75,1}的B樣條基函數(shù)圖形如圖所示.B樣條基函數(shù)圖形5.1.2分段多項(xiàng)式插值

(5.15)例5.3對(duì)下面的數(shù)據(jù)使用3次B樣條函數(shù)進(jìn)行插值.解：首先根據(jù)deBoor-Cox遞推公式編寫如下程序：x136913y25374functiony=BSpline(X,i,m,x)%X分點(diǎn)，i基函數(shù)序號(hào)，m次數(shù)，x求值點(diǎn)%y與x對(duì)應(yīng)的插值函數(shù)值X0=X;i0=i;h=X(2)-X(1);n=length(X)-1;q=0;

if(i<1)X=[((i-1):1:-1)*h+X(1),X];%在左邊擴(kuò)展點(diǎn)

q=1-i;

elseif(i+m>n)X=[X,X(end)+(1:(i+m)-n)*h];%在右邊擴(kuò)展點(diǎn)endif(m==0)y=x>=X(i+q)&x<X(i+q+1);%B_{i0}(x)elsey=(x-X(i+q))./(X(i+m+q)-X(i+q)).*BSpline(X0,i0,m-1,x)+(X(i+q+m+1)-x)./(X(i+m+1+q)-X(i+1+q)).*BSpline(X0,i0+1,m-1,x);%B_{ip}(x)end利用上面的函數(shù)可以計(jì)算出B樣條基函數(shù)在各x分點(diǎn)的函數(shù)值：x=[1,3,6,9,13];n=length(x)-1;%分段函數(shù)段數(shù)m=3;%使用3次樣條插值B=[];

%構(gòu)造系數(shù)矩陣fori=-floor((m-1)/2)-1:n+ceil((m-1)/2)-1

B=[B,BSpline(x,i,m,x)'];endBB=0.16670.69050.1429000000.25710.64290.1000000000.18750.66250.1500000000.22860.62860.1429000000.07410.59260.3333將上面計(jì)算出來的矩陣B作為系數(shù)矩陣，y作為右邊項(xiàng)，構(gòu)建線性方程組Ba=y，求解即可得到樣條插值函數(shù)f(x)的系數(shù)。y=[2,5,3,7,4]';a=linsolve(B,y)'%求解方程a=

01.41367.16780.28639.77595.52800下面就可以計(jì)算插值函數(shù)f(x)在各區(qū)間[1,13]內(nèi)的各細(xì)分點(diǎn)的函數(shù)值，并畫圖。xx=x(1):0.1:x(end);yy=zeros(size(xx));fori=-floor((m-1)/2)-1:n+ceil((m-1)/2)-1yy=yy+a(i+1-(-floor((m-1)/2)-1))*BSpline(x,i,m,xx);endplot(x,y,'o',xx,yy);注意，上面的程序并沒有使用邊界條件，所以系數(shù)矩陣B是5行7列的矩陣，方程Ba=y并沒有唯一解，程序中使用linsolve函數(shù)求出來的是最小二乘解。另外，還可以使用MATLAB提供的spline或者interp1函數(shù)來實(shí)現(xiàn)樣條函數(shù)插值。5.1.3曲線擬合

5.1.3曲線擬合

5.1.3曲線擬合

5.1.3曲線擬合

x0.0000.1000.2000.3000.4000.5000.6000.7000.8000.9001.000y0.4471.9783.2806.1607.0807.3407.6609.5609.4809.30011.20x=(0:0.1:1)';y=[0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2]';f=@(x)[x.^2,x,ones(size(x))];%擬合基函數(shù)x^2,x,1A=f(x);%系數(shù)矩陣a=(A'*A)^(-1)*(A'*y)%求解線性方程組Aa=yyy=polyval(a,x);%求擬合函數(shù)在x處的函數(shù)值plot(x,y,'+',x,yy);%畫圖MATLAB程序?yàn)椋簲M合結(jié)果5.1.3曲線擬合還可以使用MATLAB提供的polyfit函數(shù)進(jìn)行多項(xiàng)式擬合。更一般的非線性擬合可以使用lsqcurvefit、lsqnonlin函數(shù)以及應(yīng)用程序CurveFitting來進(jìn)行。5.2非線性方程求根在解決科學(xué)、工程問題中,常常會(huì)遇到求解線性方程或非線性方程,對(duì)于線性方程和部分非線性方程,人們已經(jīng)能夠求出解析解,但仍然有很多非線性方程無法求出解析解,這時(shí)就需要使用數(shù)值解法來求解.數(shù)值解法有多種,每種都有一定的使用條件、求解精度和計(jì)算量,要根據(jù)實(shí)際問題選擇合適的方法.非線性方程的根可以是實(shí)數(shù),也可以是復(fù)數(shù),本節(jié)介紹的數(shù)值解法只限于求方程的實(shí)數(shù)根.5.2.1有限窮舉法窮舉法的基本思想是對(duì)所有可能的情況都一一進(jìn)行考察,從而找出問題的解.有限窮舉法是先對(duì)求解區(qū)間[a,b]進(jìn)行剖分,得到(x0,x1,x2,?,xn),對(duì)于預(yù)先設(shè)定的精度ε>0和判別參數(shù)0<δ<b?a,依次計(jì)算yi

=|f(xi)|,i=0,1,?,n.ε-εx'1x'2f(x)如果以yi=0為條件判斷xi是否是方程的根,則由于計(jì)算精度問題無法求出根.而如果僅以|yi|<ε就判定xi是方程的根,則在根的附近會(huì)出現(xiàn)很多滿足條件的“偽根”.如圖所示,若xi∈(x'1,x'2),則有|f(xi)|<ε,如果將其中的所有xi都判定為f(x)的根顯然是不合適的。5.2.1有限窮舉法判定條件：(1)

若y0=y1=?=yn

<ε,則所有xi都是方程的根.否則根據(jù)(2)、(3)、(4)進(jìn)行判斷.(2)

若y0<ε且y0<y1,則x0是方程的根.(3)

當(dāng)0<i<n時(shí),若yi<ε,且yi<yi?1,yi<yi+1,則xi是方程的根;若yi=yi+1=?=yi+k<ε,且yi<yi?1或yi<yi+k+1,若xi+k?xi≥δ,則xi,xi+1,?,xi+k都是方程的根,若xi+k?xi<δ,則選擇xi,xi+1,?,xi+k中的任一個(gè)作為方程的根;(4)

若yn<ε且yn<yn?1,則xn是方程的根.5.2.1有限窮舉法因?yàn)閥3

<ε

,y3<

y2,y3<y4,所以x3是根因?yàn)閥14

<ε

,y14

<y13,y14

<y15,所以x14是根f(x)x3x145.2.1有限窮舉法例5.5使用窮舉法求方程xsin(x)=0在區(qū)間[?10,10]內(nèi)的實(shí)根.f=@(x)x.*sin(x);RootExhaustive(f,-10,10,1E-4,2E-4,1E-2)%調(diào)用有限窮舉法進(jìn)行求解ans=-9.4248-6.2832-3.141603.14166.28329.42485.2.1有限窮舉法當(dāng)區(qū)間剖分的步長較大時(shí),有限窮舉法可能會(huì)“漏掉”一些根.可以對(duì)有限窮舉法進(jìn)行改進(jìn),若發(fā)現(xiàn)相鄰兩個(gè)函數(shù)值f(xi)和f(xi+1)異號(hào),則對(duì)區(qū)間[xi,xi+1]進(jìn)行細(xì)分.細(xì)分方法有多種,二分法就是一種比較好的方法.例5.6使用改進(jìn)的窮舉法求方程xsin(x)=0在區(qū)間[?10,10]內(nèi)的實(shí)根.f=@(x)x.*sin(x);x1=RootExhaustive(f,-10,10,1E-2,2E-4,1E-5)%調(diào)用有限窮舉法進(jìn)行求解x2=RootExhaustive2(f,-10,10,1E-2,2E-4,1E-5)%調(diào)用改進(jìn)的有限窮舉法進(jìn)行求解x1=0x2=-9.4248-6.2832-3.141603.14166.28329.4248

5.2.2二分法設(shè)函數(shù)f(x)在[a,b]上連續(xù),且f(a)·f(b)<0,則f(x)在[a,b]內(nèi)一定有實(shí)根.如果f(x)在[a,b]內(nèi)只有唯一的單實(shí)根,則可以使用下面的方法求解:取x0=(a+b)/2,則x0將區(qū)間[a,b]分為前后兩個(gè)小區(qū)間,如果f(a)·f(x0)≤0,則根在區(qū)間[a,x0]上,令a1=a,b1=x0;否則根在[x0,b]上,令a1=x0,b1=b.在新的區(qū)間[a1,b1]上用同樣的方法搜索.由于新的區(qū)間每次減半,并且保證方程的根總是在新的區(qū)間里面,所以當(dāng)區(qū)間寬度小于預(yù)先給定的精度ε時(shí),區(qū)間內(nèi)的任意一點(diǎn)都是滿足精度要求的近似解.5.2.1有限窮舉法f(x)ab(a+b)/2f(x)a1b1f(x)a1b1(a1+b1)/2f(x)b2a2f(a)·f((a+b)/2)<0f(b1)·f((a1+b1)/2)<05.2.2二分法例5.7使用二分法求方程x3?x?1=0在區(qū)間[1,1.5]內(nèi)的實(shí)根.f=@(x)x.^3-x-1;RootBisection(f,1,1.5,0.01)%調(diào)用二分法進(jìn)行求解%RootBisectionRecursive(f,1,1.5,0.01)%二分法的遞歸實(shí)現(xiàn)ans=1.3242二分法求方程的根5.2.3不動(dòng)點(diǎn)迭代法5.2.3不動(dòng)點(diǎn)迭代法不動(dòng)點(diǎn)迭代法求方程f(x)=0的根是將方程等價(jià)地改寫成x=g(x)的形式,然后構(gòu)造序列xk+1=g(xk)如果序列{xk}收斂于x*,則x*是方程x=g(x)的解,從而也是f(x)=0的解.對(duì)于預(yù)先給定的精度ε,必然存在正整數(shù)N,對(duì)任意的n>N時(shí),有|

xk

?x*

|<

ε也就是說,如果只是求近似解,對(duì)任意的xn(n>N)都是滿足精度的解.5.2.3不動(dòng)點(diǎn)迭代法

functionx=RootIteration(g,x0,e)%使用不動(dòng)點(diǎn)迭代法求x=g(x)的根%g迭代函數(shù)%x0迭代初值%e精度%x根if(nargin<3),e=1e-5;end

%設(shè)定默認(rèn)精度x=g(x0);holdonplot([x0xx],[xxg(x)],'--r');%作圖,僅用于演示while(abs(x-x0)>=e)x0=x;x=g(x0);%迭代

plot([x0xx],[xxg(x)],'--r');%作圖,僅用于演示endholdoffend不動(dòng)點(diǎn)迭代法程序代碼如下：5.2.3不動(dòng)點(diǎn)迭代法例5.8使用不動(dòng)點(diǎn)迭代法求方程x?e?x=0的根.g=@(x)exp(-x);fplot(g,[-1,4]);%畫函數(shù)圖形holdonplot([-14],[-14],'g');%畫直線x=yx=RootIteration(g,3)%調(diào)用牛頓法求根xlabel('x');ylabel('g(x)');holdoffx=0.5671不動(dòng)點(diǎn)求根過程5.2.4牛頓法

5.2.4牛頓法牛頓法是一種用切線近似曲線的方法,它的幾何意義是過點(diǎn)(xk,f(xk))作函數(shù)的切線,與x軸相交,交點(diǎn)的橫坐標(biāo)為,將此作為新的近似解即可得到牛頓迭代公式.牛頓法的幾何意義functionx=RootNewton(f,x0,e)%使用牛頓法求f(x)=0的根%f函數(shù)句柄%x0迭代初值%e精度%x根if(nargin<3),e=1e-5;end

%設(shè)定默認(rèn)精度symst;%定義符號(hào)變量df=@(x)subs(diff(f(t)),t,x);%求導(dǎo)g=@(x)x-f(x)/double(df(x));%牛頓公式holdonx=g(x0);plot([x0x0x],[0f(x0)0],'--r');%作圖,僅用于演示while(abs(x-x0)>=e)%如果相鄰兩次結(jié)果差別小于精度要求,則結(jié)束

x0=x;%保存上一次結(jié)果

x=g(x0);%迭代

plot([x0x0x],[0f(x0)0],'--r');%作圖,僅用于演示endholdoff牛頓法程序代碼如下:5.2.4牛頓法例5.9使用牛頓法求方程f(x)=x2-4sin(x)的根.f=@(x)x.^2-4*sin(x);fplot(f,[1,5]);%畫函數(shù)圖形holdonplot([15],[00],'k');%畫橫坐標(biāo)x=RootNewton(f,4)%調(diào)用牛頓法求根holdonplot(x,f(x),'+r');%顯示根的位置x=1.93385.2.5雙點(diǎn)弦截法

5.2.5雙點(diǎn)弦截法

雙點(diǎn)弦截法的幾何意義f(x)xk-1xkxk+15.2.5雙點(diǎn)弦截法例5.10使用雙點(diǎn)弦截法求方程ex?3x2

=0的根.f=@(x)exp(x)-3*x.^2;%待求解方程x=RootSecant(f,0,1)%使用雙點(diǎn)弦截法求解x

=0.9100雙點(diǎn)弦截法求根過程5.3迭代法求線性方程組數(shù)值解當(dāng)線性方程組的規(guī)模較大時(shí),使用直接法需要的計(jì)算量太大,有時(shí)無法在有效的時(shí)間內(nèi)求解出來,此時(shí)可采用迭代法求近似解.迭代法的基本思想是:構(gòu)造一個(gè)可收斂到精確解的迭代公式,從事先設(shè)定的初值開始迭代,逐步逼近精確解,當(dāng)滿足精度要求時(shí)停止迭代,此時(shí)即可得到近似解.常見的迭代格式有Jacobi迭代、Gauss-Seidel迭代和超松馳迭代法.5.3.1Jacobi迭代法

5.3.1Jacobi迭代法

5.3.1Jacobi迭代法

Jacobi迭代法求線性方程組的解的程序如下:functionx=LineSolveJacobi(A,b,eps)%使用雅可比迭代法求線性方程組的解%A系數(shù)矩陣%b方程右邊向量%eps精度%x近似解if(nargin<3)eps=1e-5;%設(shè)置默認(rèn)精度end[~,n]=size(A);x=rand(n,1);

%從隨機(jī)向量開始迭代while(1)x0=x;

fori=1:ns=A(i,:)*x0-A(i,i)*x0(i);x(i)=(b(i)-s)/A(i,i);%使用Jacobi迭代

end

if(max(abs(x-x0))<eps)%滿足終止條件則退出

break;

end

if(any(abs(x)>1e10))%檢測是否發(fā)散

x=zeros(size(b));disp('迭代過程不收斂');

return;

endend

A=[10,-2,3,-2;1,-15,-2,-1;2,4,8,1;-3,-4,9,20];%系數(shù)矩陣b=[3;-11;16;10];

%右邊項(xiàng)x1=LineSolveJacobi(A,b)

%使用Jacobi迭代法求解x1=-0.20470.49201.8376-0.25925.3.2Gauss-Seidel迭代法

Gauss-Seidel迭代法求線性方程組的解的程序如下:functionx=LineSolveGaussSeidel(A,b,eps)%使用Gauss-Seidel迭代法求線性方程組的解%A系數(shù)矩陣%b方程右邊向量%eps精度%x近似解if(nargin<3)eps=1e-5;%設(shè)置默認(rèn)精度end[~,n]=size(A);x=rand(n,1);%從隨機(jī)向量開始迭代while(1)x0=x;

fori=1:ns=A(i,1:i-1)*x(1:i-1)+A(i,i+1:n)*x0(i+1:n);x(i)=(b(i)-s)/A(i,i);%使用Gauss-Seidel迭代

end

if(max(abs(x-x0))<eps),break;end

if(any(abs(x)>1e10))%檢測是否發(fā)散

x=zeros(size(b));disp('迭代過程不收斂');

return;

endend

A=[10,-2,3,-2;1,-15,-2,-1;2,4,8,1;-3,-4,9,20];%系數(shù)矩陣b=[3;-11;16;10];

%右邊項(xiàng)x2=LineSolveGaussSeidel(A,b)

%使用Gauss-Seidel迭代法求解x2=-0.20470.49201.8376-0.25925.3.3超松馳迭代法

超松馳迭代法求線性方程組的解的程序如下:functionx=LineSolveSOR(A,b,w,eps)%使用松馳迭代法求線性方程組的解%A系數(shù)矩陣%b方程右邊向量%w松馳因子%eps精度,當(dāng)相鄰兩次迭代結(jié)果的差小于eps時(shí)結(jié)束%x近似解if(nargin<3)w=1.3;%設(shè)置默認(rèn)松馳因子endif(nargin<4)eps=1e-5;%設(shè)置默認(rèn)精度end[~,n]=size(A);x=rand(n,1);

%從隨機(jī)向量開始迭代while(1)x0=x;

fori=1:ns=A(i,1:i-1)*x(1:i-1)+A(i,i:n)*x0(i:n);x(i)=x0(i)+w*(b(i)-s)/A(i,i);

end

if(max(abs(x-x0))<eps)%滿足終止條件則退出

break;

end

if(any(abs(x)>1e10))%檢測是否發(fā)散x=zeros(size(b));

disp('迭代過程不收斂');

return;

endend

A=[10,-2,3,-2;1,-15,-2,-1;2,4,8,1;-3,-4,9,20];%系數(shù)矩陣b=[3;-11;16;10];

%右邊項(xiàng)x3=LineSolveSOR(A,b)

%使用超松馳迭代法求解x3=-0.20470.49201.8376-0.2592上面使用了3種方法進(jìn)行求解,求解結(jié)果都相同.一般情況下,Jacobi迭代法要比Gauss-Seidel迭代法的迭代次數(shù)多,SOR迭代法的迭代次數(shù)與松馳因子有關(guān),需要進(jìn)行多次嘗試才能得到比較好的因子.5.4無約束一維極值問題無約束一維最小值問題可簡單表述為minf(x),x∈[a,b]滿足上面問題的解x*稱為函數(shù)f(x)在區(qū)間[a,b]上的一個(gè)最小值點(diǎn).如果f(x)在x*的一個(gè)鄰域D有定義,對(duì)任意x∈D,都有f(x*)≤f(x),則稱x*是函數(shù)f(x)的一個(gè)極小值點(diǎn).最小值點(diǎn)是全局的,極小值是局部的.同樣的方法也可以定義最大值點(diǎn)和極大值點(diǎn).5.4無約束一維極值問題

5.4.1有限窮舉法有限窮舉法是對(duì)求解區(qū)間[a,b]進(jìn)行剖分,得到(x0,x1,x2,?,xn),對(duì)每一個(gè)分點(diǎn)xi,依次計(jì)算yi=f(xi),求最小的yi及與之對(duì)應(yīng)的xi即可得到最小值點(diǎn).窮舉法的優(yōu)點(diǎn)是對(duì)函數(shù)的要求不高,只要函數(shù)連續(xù)即可.缺點(diǎn)是計(jì)算量較大,精度較低.有限窮舉法的程序如下:function[minx,miny]=MinExhaustive(f,a,b,e)%使用有限窮舉法求函數(shù)最小值%f目標(biāo)函數(shù)(句柄)%[a,b]

搜索區(qū)間%e求解精度%minx最小值點(diǎn)%miny最小值if(nargin<4),e=1E-2;end

%默認(rèn)精度為1E-2x=a:e:b;%對(duì)區(qū)間進(jìn)行剖分y=f(x);%計(jì)算y值[miny,mini]=min(y);%求最小值及對(duì)應(yīng)的下標(biāo)minx=x(mini);%求最小的y值對(duì)應(yīng)的x值fplot(f,[a,b]);

%繪制函數(shù)圖形holdonplot(minx,miny,'+r');%顯示極小值點(diǎn)holdoffend注:程序中首先求出所有剖分分點(diǎn)的函數(shù)值,然后利用內(nèi)置函數(shù)min()求出最小值miny及對(duì)應(yīng)的下標(biāo)mini,再根據(jù)下標(biāo)獲取對(duì)應(yīng)的極小值點(diǎn)minx.作圖代碼僅用于演示,在實(shí)際使用時(shí)可刪除.例5.12使用窮舉法求函數(shù)f(x)=?(x+2sinx

+3cosx)在區(qū)間[?4,4]內(nèi)的最小值點(diǎn).f=@(x)-x-2*sin(x)-3*cos(x);MinExhaustive(f,-4,4)ans=0.8700求解結(jié)果5.4.2進(jìn)退法如果f(x)在區(qū)間[a,b]內(nèi)連續(xù)且只存在一個(gè)極值點(diǎn)時(shí),可以使用進(jìn)退法求解.給定初值x0∈[a,b],步長h(0<h<x0?a且h<b?x0).(1)

如果f(x0)<f(x0+h),則后退到x0?kh,使得f(x0?kh)>f(x0),k為滿足條件的最小的正整數(shù),新的搜索區(qū)間為[x0?kh,x0+h].如果找不到這樣的k,則新的搜索區(qū)間為[a,x0+h].(2)

如果f(x0)>f(x0+h),則前進(jìn)到x0+kh,使得f(x0)>

f(x0+kh),k為滿足條件的最小的正整數(shù).新的搜索區(qū)間為[x0,x0+kh].如果找不到這樣的k,則新的搜索區(qū)間為[x0,b].(3)如果新的搜索區(qū)間寬度小于預(yù)先設(shè)定的精度,結(jié)束搜索過程;否則縮小步長h,在新的區(qū)間內(nèi)繼續(xù)搜索;進(jìn)退法的程序如下:function[x,y]=MinAdvanceRetreat(f,a,b,e)%使用進(jìn)退法求函數(shù)最小值%f目標(biāo)函數(shù)(句柄)%[a,b]搜索區(qū)間%e求解精度%x最小值點(diǎn)%y最小值if(nargin<4),e=1E-5;end

%默認(rèn)精度為1E-5fplot(f,[ab]);

%作函數(shù)圖形holdon;x=a+(b-a)*rand(1);

%在搜索區(qū)間內(nèi)隨機(jī)生成一個(gè)初始點(diǎn)h=(b-a)/5;

%步長為區(qū)間寬度的1/5while(b-a>e)y=f(x);

if(y<f(x+h))

%后退，搜索滿足f(x-k*h)>=f(x)的最小的正整數(shù)kx1=x-h:-h:a;y1=f(x1);k=find(y1>=y,1);

if(~isempty(k))a=x1(k);

endb=x+h;%修改搜索區(qū)間

else

%前進(jìn)，搜索滿足f(x+k*h)>=f(x)的最小的正整數(shù)kx1=x+h:h:b;y1=f(x1);k=find(y1>=y,1);

if(~isempty(k))b=x1(k);

enda=x;%修改搜索區(qū)間

endx=(a+b)/2;

%取區(qū)間中點(diǎn)的新的解

h=(b-a)/5;

%收縮步長

text(a,f(a),'[');text(b,f(b),']');%顯示搜索區(qū)間

%pause(0.5);

%暫停endy=f(x);plot(x,y,'+r');%顯示極小值點(diǎn)holdoffend例5.13使用進(jìn)退法求函數(shù)f(x)=x2+2x+1在區(qū)間[?3,1]上的極小值點(diǎn).f=@(x)x.^2+2*x+1;MinAdvanceRetreat(f,-3,1)ans=-1.0000求解結(jié)果5.4.3黃金分割法黃金分割法也叫0.618法,也是只適合f(x)在區(qū)間[a,b]內(nèi)連續(xù)且只存在一個(gè)極值點(diǎn)的情形.基本思想是,利用黃金分割點(diǎn)產(chǎn)生[a,b]內(nèi)的兩個(gè)分點(diǎn)x1,x2,即x1=a+0.382(b?a)x2=a+0.618(b?a)(1)如果f(x1)<f(x2),則搜索區(qū)間為[a,x2];(2)如果f(x1)>f(x2),則搜索區(qū)間為[x1,b];然后在新的搜索區(qū)間用同樣的方法繼續(xù)搜索,直到區(qū)間寬度小于預(yù)先給定的精度為止.黃金分割法的程序如下:function[x,y]=MinGoldenSection(f,a,b,e)%使用黃金分割法求函數(shù)最小值%f目標(biāo)函數(shù)(句柄)%[a,b]搜索區(qū)間%e求解精度%x最小值點(diǎn)%y最小值if(nargin<4),e=1E-5;end

%默認(rèn)精度為1E-5fplot(f,[ab]);

%作函數(shù)圖形holdon;while(b-a>e)x1=a+0.382*(b-a);

%分點(diǎn)1x2=a+0.618*(b-a);%分點(diǎn)2

if(f(x1)<f(x2))b=x2;

elsea=x1;

endtext(a,f(a),'[');text(b,f(b),']');%顯示搜索區(qū)間

%pause(0.5);

%暫停endx=(a+b)/2;y=f(x);plot(x,y,'+r');

%顯示極小值點(diǎn)holdoffend例5.14使用黃金分割法求函數(shù)f(x)=x2+2x+1在區(qū)間[?3,1]上的極小值點(diǎn).f=@(x)x.^2+2*x+1;MinGoldenSection(f,-3,1)ans=-1.0000求解結(jié)果5.4.4斐波那契法

5.4.4斐波那契法

斐波那契法的程序如下:function[minx,miny]=MinFibonacci(f,a,b,e)%使用斐波那契法求函數(shù)最小值%f目標(biāo)函數(shù)(句柄)%[a,b]搜索區(qū)間%e求解精度%minx最小值點(diǎn)%miny最小值if(nargin<4),e=1E-5;end

%默認(rèn)精度為1E-5fplot(f,[a,b]);holdon;

%繪制函數(shù)圖形text(a,f(a),'[');text(b,f(b),']');%顯示搜索區(qū)間%%計(jì)算Fibonacc數(shù)列前n項(xiàng)F(1)=1;F(2)=1;n=2;N=(b-a)/e;while(F(n)<N)n=n+1;F(n)=F(n-1)+F(n-2);endk=n;a1=a+F(k-2)/F(k)*(b-a);b1=a+F(k-1)/F(k)*(b-a);ya1=f(a1);yb1=f(b1);%計(jì)算分點(diǎn)的函數(shù)值while(k>2)

%迭代n-2次

k=k-1;

if(ya1<yb1)%當(dāng)左分點(diǎn)函數(shù)值小于右分點(diǎn)函數(shù)值時(shí),將搜索區(qū)間設(shè)在[a,b1]b=b1;%設(shè)置新區(qū)間的右端點(diǎn)

b1=a1;%新區(qū)間右分點(diǎn)為原區(qū)間的左分點(diǎn)

yb1=ya1;a1=a+F(k-2)/F(k)*(b-a);%計(jì)算新區(qū)間的左分點(diǎn)

ya1=f(a1);

else

%當(dāng)左分點(diǎn)函數(shù)值大于等于右分點(diǎn)函數(shù)值時(shí),將搜索區(qū)間設(shè)在[a1,b]a=a1;%設(shè)置新區(qū)間的左端點(diǎn)

a1=b1;%新區(qū)間左分點(diǎn)為原區(qū)間的右分點(diǎn)

ya1=yb1;b1=a+F(k-1)/F(k)*(b-a);%計(jì)算新區(qū)間的右分點(diǎn)

yb1=f(b1);

end

text(a,f(a),'[');text(b,f(b),']');

%顯示搜索區(qū)間endminx=(a+b)/2;%取最終區(qū)間的中點(diǎn)為最小值點(diǎn)miny=f(minx);plot(minx,miny,'+r');holdoff%顯示極小值點(diǎn)end斐波那契法的程序如下(遞歸版):function[minx,miny]=MinFibonacci2(f,a,b,e)%使用斐波那契法求函數(shù)最小值%f目標(biāo)函數(shù)(句柄)%[a,b]搜索區(qū)間%e求解精度%minx最小值點(diǎn)%miny最小值if(nargin<4),e=1E-5;end

%默認(rèn)精度為1E-5if(b-a<=e)%當(dāng)b-a<=e時(shí)結(jié)束遞歸過程

minx=(a+b)/2;%取最終區(qū)間的中點(diǎn)為最小值點(diǎn)

miny=f(minx);

return;else

%%計(jì)算Fibonacc數(shù)列前n項(xiàng)F(1)=1;F(2)=1;n=2;

while(F(n)<(b-a)/e)n=n+1;F(n)=F(n-1)+F(n-2);

end

%%a1=a+F(n-2)/F(n)*(b-a);%計(jì)算兩個(gè)分點(diǎn)

b1=a+F(n-1)/F(n)*(b-a);

if(f(a1)<f(b1))%當(dāng)左分點(diǎn)函數(shù)值小于右分點(diǎn)函數(shù)值時(shí),將搜索區(qū)間設(shè)在[a,b1]b=b1;

else

%當(dāng)左分點(diǎn)函數(shù)值大于等于右分點(diǎn)函數(shù)值時(shí),將搜索區(qū)間設(shè)在[a1,b]a=a1;

end[minx,miny]=MinFibonacci2(f,a,b,e);%遞歸endend遞歸版的程序比較簡潔，但執(zhí)行效率較低.例5.14使用斐波那契法求函數(shù)f(x)=x2+2x+1在區(qū)間[?3,1]上的極小值點(diǎn).f=@(x)x.^2+2*x+1;MinFibonacci(f,-3,1)ans=-1.0000求解結(jié)果5.4.5牛頓法

牛頓法的程序如下:function[x,y]=MinNewton(f,a,b,e)%使用牛頓法求函數(shù)最小值%f目標(biāo)函數(shù)(句柄)%[a,b]搜索區(qū)間%e求解精度%x最小值點(diǎn)%y最小值if(nargin<4),e=1E-5;end

%默認(rèn)精度為1E-5fplot(f,[ab]);%繪制函數(shù)圖形holdon;symst;%定義符號(hào)變量symsdf=diff(f(t),t);%求一階導(dǎo)函數(shù)symsdf2=diff(symsdf,t);%求二階導(dǎo)函數(shù)df1=@(x)subs(symsdf,t,x);df2=@(x)subs(symsdf2,t,x);x=(a+b)/2;%使用區(qū)間中點(diǎn)做為迭代初值dy=double(df1(x));while(abs(dy)>e)x=x-dy/double(df2(x));%迭代公式

dy=double(df1(x));plot(x,f(x),'o');

%pause;

%暫停endy=f(x);plot(x,y,'+r');holdoffend程序中使用了diff函數(shù),diff(expr,t)可以求expr的導(dǎo)數(shù),其中expr是關(guān)于t的符號(hào)表達(dá)式,t是符號(hào)變量.subs表示代入并計(jì)算結(jié)果,subs(expr,t,x)表示將表達(dá)式expr中的變量t替換成x,并計(jì)算表達(dá)式的值.例5.16使用牛頓法求函數(shù)f(x)=x4+2x+1在區(qū)間[?2,1]上的極小值點(diǎn).f=@(x)x.^4+2*x+1;MinNewton(f,-2,1)ans=-0.7937求解結(jié)果5.4.6割線法

割線法的程序如下:function[x1,y]=MinSecant(f,a,b,e)%使用割線法求函數(shù)最小值%f目標(biāo)函數(shù)(句柄)%[a,b]搜索區(qū)間%e求解精度%x1最小值點(diǎn)%y最小值if(nargin<4),e=1E-5;end

%默認(rèn)精度為1E-5fplot(f,[ab]);%繪制函數(shù)圖形holdon;symst;%定義符號(hào)變量symsdf=diff(f(t),t);%求一階導(dǎo)函數(shù)df1=@(x)subs(symsdf,t,x);x0=(a+b)/3;%初值1x1=2*(a+b)/3;%初值2dy0=double(df1(x0));dy1=double(df1(x1));while(abs(dy1)>e)dy1=double(df1(x1));x2=x1-(x1-x0)*dy1/(dy1-dy0);%迭代公式

plot(x2,f(x2),'o');x0=x1;%更新

x1=x2;dy0=dy1;endy=f(x1);plot(x1,y,'+r');%顯示解的位置holdoffend例5.17使用割線法求函數(shù)f(x)=x4+2x+1在區(qū)間[?2,1]上的極小值點(diǎn).f=@(x)x.^4+2*x+1;MinSecant(f,-2,1)ans=-0.7937求解結(jié)果5.5無約束多維極值無約束多維極值問題與無約束一維極值問題類似,都是求解函數(shù)f(x)的極值,不同的是這里的x是向量.其一般形式為

minf(x),x∈Rn求解這類問題的方法非常多,總體上分為直接法和使用導(dǎo)數(shù)計(jì)算的間接法.直接法有模式搜索法、Rosenbrock法、單純形法、Powell法等,間接法有最速下降法、共軛梯度法、牛頓法、修正牛頓法、擬牛頓法、信賴域法、顯示最速下降法等.5.5.1模式搜索法(Hooke-Jeeves法)模式搜索法主要由兩種移動(dòng)過程組成:探測移動(dòng)和模式移動(dòng).探測移動(dòng)是沿坐標(biāo)軸方向的移動(dòng),模式移動(dòng)則是沿相鄰兩個(gè)探測點(diǎn)連線方向上移動(dòng).兩種移動(dòng)交替進(jìn)行,目的是尋找使得函數(shù)下降的最佳方向.5.5.1模式搜索法(Hooke-Jeeves法)模式搜索法的算法如下:(1)給定初始點(diǎn)x(0),初始步長δ>0,加速系數(shù)α>0,收縮系數(shù)θ∈(0,1)及精ε度.令k=0;(2)置x(k

+1)=x(k).從x(k)出發(fā),依次沿坐標(biāo)軸方向探測移動(dòng):正向探測:若f(x(k)+δej)<f(x(k)),則令x(k+1)=x(k)+

δej,否則做負(fù)向探測;負(fù)向探測:若f(x(k)?

δej)<f(x(k)),則令x(k+1)=x(k)?

δej;其中ej為平行于第j個(gè)坐標(biāo)軸方向的單位向量,j=1,2,?,n;(3)若f(x(k+1))<f(x(k)),對(duì)x(k+1)沿p(k)

=x(k+1)?x(k)方向作模式移動(dòng),x(k+1)

=x(k+1)

+α

p(k);(4)當(dāng)f(x(k+1))≥f(x(k))時(shí),如果步長δ>ε則減小步長δ=θδ,k=k+1,轉(zhuǎn)到第(2)步繼續(xù)搜索;當(dāng)步長δ<ε時(shí)結(jié)束搜索過程;模式搜索法程序代碼如下:function[X,k]=MinHookeJeeves(f,X0,delta,sita,alpha,eps)%使用模式搜索法求多元函數(shù)極小值%f目標(biāo)函數(shù)%X0初始點(diǎn)%delta初始步長%sita收縮系數(shù)%alpha加速系數(shù)%e求解精度%X極小值點(diǎn)%k迭代次數(shù)if(nargin<6),eps=1E-5;end

%默認(rèn)精度為1E-5k=0;n=length(X0);%根據(jù)x0的維數(shù)獲取函數(shù)f的維數(shù)while(1)k=k+1;y0=f(X0);X1=X0;y1=f(X1);

fori=1:nt=zeros(size(X0));t(i)=delta;%探測步長

if(f(X1+t)<y1)%如果向前的結(jié)果更好

X1=X1+t;%更新位置

y1=f(X1);%更新函數(shù)值

elseif(f(X1-t)<y1)%如果向后的結(jié)果更好

X1=X1-t;%更新位置

y1=f(X1);%更新函數(shù)值

end

end

if(y1<y0)%如果新的位置更好

if(f(X0+alpha*(X1-X0))<y1)X0=X0+alpha*(X1-X0);%更新位置

elseX0=X1;

end

else

if(delta<eps)%如果滿足精度要求

X=X0;%將X0作為最優(yōu)解

break;%結(jié)束

enddelta=sita*delta;%縮小搜索步長

endend

ezmesh(@(x,y)x.^4+y.^4+4*x.*y-4*x-3*y,[-2,2,-2,2]);%繪制曲面holdon;f=@(x)x(:,1).^4+x(:,2).^4+4*x(:,1)*x(:,2)-4*x(:,1)-3*x(:,2);%f(x)是多元函數(shù)，參數(shù)x的每一列為一個(gè)自變量X=MinHookeJeeves(f,[00],1,0.05,1.5)%調(diào)用模式搜索法求最小值plot3(X(1),X(2),f(X),'*r');holdoffview(127,54)%調(diào)整圖形視圖角度求解結(jié)果5.5.2Rosenbrock法Rosenbrock法是一種轉(zhuǎn)軸法,其基本思想是在當(dāng)前點(diǎn)構(gòu)造n個(gè)正交方向,然后在每個(gè)方向做探測移動(dòng),找到函數(shù)值下降最快的方向,移動(dòng)某個(gè)步長,再在新的點(diǎn)構(gòu)造n個(gè)正交方向,如此循環(huán).算法如下:(1)給定初始點(diǎn)x(0),n個(gè)單位正交方向向量dj,j=1,2,?,n,初始步長δ0>0,放大系數(shù)α>1,收縮系數(shù)θ∈(0,1)及精度ε.令k=0;(2)置x(k+1)=x(k),δ=δ0.從x(k)出發(fā),依次沿dj方向作反復(fù)探測移動(dòng).若f(x(k)+δdj)<f(x(k)),則令x(k+1)=x(k)+δdj,δ=αδ;若f(x(k)+δdj)>f(x(k)),則令δ=?θδ;直到n個(gè)方向都在成功后出現(xiàn)了失敗.(3)如果||x(k+1)?x(k)||<ε,停止搜索過程,輸出x(k+1);否則用x(k+1)?x(k)替換n個(gè)單位正交方向向量dj中的某一個(gè),替換后應(yīng)使dj線性無關(guān).使用Schmidt正交化過程建立新的標(biāo)準(zhǔn)正交化向量,令k=k+1,轉(zhuǎn)到第(2)步繼續(xù)搜索.Schmidt正交化過程

f=@(x)x(:,1).^4+x(:,2).^4+4*x(:,1)*x(:,2)-4*x(:,1)-3*x(:,2);X0=[0,0];D=[10;01];delta=[1,1];theta=[0.5,0.5];alpha=[1.5,1.5];eps=1E-5;[X,k]=MinRosenBrock(f,X0,D,delta,theta,alpha,eps)X=

1.2101-0.7720k=

335.5.3最速下降法

最速下降法程序代碼如下:function

[x,k]=MinGradientDescent(f,x0,eps)%使用最速下降法求多元函數(shù)極小值%f目標(biāo)函數(shù)%x0迭代初值%eps預(yù)設(shè)精度%x最優(yōu)解%k迭代次數(shù)n=length(x0);t=sym('t',[1n]);

%定義符號(hào)向量gradf=gradient(f(t),t)';%梯度v=[1,0];k=1;x=x0;while(norm(v)>eps)v=double(-subs(gradf,t,x));%計(jì)算梯度向量在當(dāng)前位置x的值

psi=@(lambda)f(x+lambda*v);%構(gòu)造一維極值問題目標(biāo)函數(shù)

lambda=fminbnd(psi,-1,1);%求解最優(yōu)步長

x=x+lambda*v;%更新當(dāng)前位置

k=k+1;end

f=@(x)x(:,1).^4+x(:,2).^4+4*x(:,1)*x(:,2)-4*x(:,1)-3*x(:,2);X0=[00];eps=1E-5;[x,k]=MinGradientDescent(f,X0,eps)x=

1.2101-0.7720k=

245.5.4Newton法

5.5.4Newton法

Newton法程序代碼如下:function[x,k]=MinMultipleNewton(f,x0,eps)%使用牛頓法求多元函數(shù)極小值%f目標(biāo)函數(shù)%x0迭代初值%eps預(yù)設(shè)精度%x最優(yōu)解%k迭代次數(shù)n=length(x0);t=sym('t',[1n]);%定義符號(hào)向量gradf=gradient(f(t),t);%梯度jacf=jacobian(gradf,t);%雅可比矩陣v=[10];k=1;x=x0;while(norm(v)>eps)v=double(-subs(gradf,t,x));%計(jì)算梯度向量在當(dāng)前位置x的值

J=double(-subs(jacf,t,x));%計(jì)算雅可比矩陣在當(dāng)前位置x的值

p=-inv(J)*v;%搜索方向

x=x+p';%更新當(dāng)前位置

k=k+1;end

f=@(x)x(:,1).^4+x(:,2).^4+4*x(:,1)*x(:,2)-4*x(:,1)-3*x(:,2);X0=[00];eps=1E-5;[x,k]=MinMultipleNewton(f,X0,eps)x=1.2101-0.7720k=1895.6數(shù)值積分

5.6數(shù)值積分

5.6.1矩形法

5.6.1矩形法如圖所示,左矩形法是將小矩形面積之和作為函數(shù)曲線下方面積的近似值,小矩形的高度為小區(qū)間左端點(diǎn)的函數(shù)值,計(jì)算的誤差是曲線下方與小矩形上方曲邊三角形的面積(f(x)為增函數(shù)時(shí)).顯然,剖分?jǐn)?shù)越多,小區(qū)間的寬度越小,精度越高.左矩形法求定積分示意圖左矩形法求定積分的程序如下:functionS=IntLeftBox(f,a,b,n)%使用左矩形法求f(x)在[a,b]上的定積分%f被積函數(shù)%[a,b]積分區(qū)間%n區(qū)間等分?jǐn)?shù)%S積分if(nargin<4),n=50;endfplot(f,[a,b],'k');holdon;h=(b-a)/n;S=0;fori=0:n-1x=a+h*i;y=f(x);S=S+y;fill([x,x+h,x+h,x,x],[0,0,y,y,0],[00.51],'facealpha',0.5);%畫小矩形endS=h*S;holdofffunctionS=IntLeftBox2(f,a,b,n)%使用左矩形法求f(x)在[a,b]上的定積分%f被積函數(shù)%[a,b]積分區(qū)間%n區(qū)間等分?jǐn)?shù)%S積分if(nargin<4),n=50;endh=(b-a)/n;x=linspace(a,b,n+1);y=f(x);S=h*sum(y(1:end-1));或

f=@(x)exp(x.*sin(x));S1=IntLeftBox(f,0,1)%左矩形法