Matlab解非線性超定方程組

1、文檔來源為：從網(wǎng)絡(luò)收集整理,word版本可編輯.歡迎下載支持Matlab解非線性超定方程組3x+2/(5+y)=6,4x+4/(5+y)=7 , 9x+4/(8+y)=12 11x+2/(4+y)=15x,y是未知數(shù) clc;clear;%其實樓主的問題可以等效為求最小值的問題，我使用的指標(biāo)是典型的平方和最 小x0=1,1；f=(x)(3*x(1)+2/(5+x(2)-6)A2+(4*x(1)+4/(5+x(2)-7)A2+(9*x(1)+4/(8+x(2)-12)A2+(11*x(1)+2/(4+x(2)-15)A2;x,fval=fminsearch(f,x0)=求解線性方程組solve,

2、 linsolve例：A=5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1;%矩陣的行之間用分號隔開，元素之間用逗號或空格B=3;1;1;0X=zeros(4,1);% 建立一個4元列向量X=linsolve(A,B)diff (fun , var, n):對表達(dá)式fun中的變量var求n階導(dǎo)數(shù)。例如：F=sym (u(x,y)*v(x,y), ) ; %sym ()用來定義一個符號表達(dá)式diff(F); %matlab區(qū)分大小寫pretty(ans) %pretty ():用習(xí)慣書寫方式顯示變量；ans是答案表達(dá)式非線性方程求解fsolve(fun,x0,options)其中

3、fun為待解方程或方程組的文件名;文檔來源為：從網(wǎng)絡(luò)收集整理,word版本可編輯.歡迎下載支持x0位求解方程的初始向量或矩陣；option為設(shè)置命令參數(shù)建立文件fun.m ：function y=fun(x)y=x(1)-0.5*sin(x(1)-0.3*cos(x(2), .x(2) - 0.5*cos(x(1)+0.3*sin(x(2);clear;x0=0.1,0.1;fsolve(fun,x0,optimset(fsolve)注：.為續(xù)行符m文件必須以function為文件頭，調(diào)用符為;文件名必須與定義的函數(shù)名相同； fsolve ()主要求解復(fù)雜非線性方程和方程組，求解過程是一個逼近

4、過程。Matlab求解線性方程組AX=B 或 XA=B在MATLAB中，求解線性方程組時，主要采用前面章節(jié)介紹的除法運算符 ”和如：X=AB表示求矩陣方程AX = B的解；X=B/A表示矩陣方程XA=B的解。對方程組X = AB ,要求A和B用相同的行數(shù)，X和B有相同的列數(shù)，它的行數(shù) 等于矩陣A的列數(shù)，方程X = B/A同理。如果矩陣A不是方陣，其維數(shù)是mx n,則有：m = n恰定方程，求解精確解；mn超定方程，尋求最小二乘解；mm。則方程組沒有精 確解，此時稱方程組為超定方程組。線性超定方程組經(jīng)常遇到的問題是數(shù)據(jù)的曲 線擬合。對于超定方程，在 MATLAB中，利用左除命令(x=Ab)來尋求

5、它的最 小二乘解；還可以用廣義逆來求，即 x=pinv(A),所得的解不一定滿足 Ax=b , x 只是最小二乘意義上的解。左除的方法是建立在奇異值分解基礎(chǔ)之上，由此獲得的解最可靠；廣義逆法是建立在對原超定方程直接進(jìn)行householder變換的基礎(chǔ)上，其算法可靠性稍遜與奇異值求解，但速度較快；【例7】求解超定方程組A=2 -1 3;3 1 -5;4 -1 1;1 3 -13 A=文檔來源為：從網(wǎng)絡(luò)收集整理,word版本可編輯.歡迎下載支持-1 31 -5-1 1 1 3 -13b = 3 0 3 -6； rank(A) ans= 3 x1=Ab x1 = 1.00002.00001.0000

6、 x2=pinv(A)*b x2= 1.0000 2.0000 1.0000A*x1-b ans= 1.0e-014 -0.0888 -0.0888 -0.1332 0可見x1并不是方程Ax=b的精確解，用x2=pinv(A)*b所得的解與x1相同。 三.欠定方程組欠定方程組未知量個數(shù)多于方程個數(shù)，但理論上有無窮個解。MATLAB將尋求一個基本解，其中最多只能有 m個非零元素。特解由列主元qr分解求得。文檔來源為：從網(wǎng)絡(luò)收集整理,word版本可編輯.歡迎下載支持【例8】解欠定方程組A=1 -2 1 1;1 -2 1 -1;1 -2 1 5 A=1 -2 1 11 -2 1 -11 -2 1 -

7、11 -2 1 5b=1 -1 5x1=AbWarning:Rank deficient,rank=2 tol=4.6151e-015 x1 = 0-0.000001.0000 x2=pinv(A)*b x2= 0-0.00000.00001.0000四.方程組的非負(fù)最小二乘解在某些條件下，所求的線性方程組的解出現(xiàn)負(fù)數(shù)是沒有意義的。 雖然方程組可以 得到精確解，但卻不能取負(fù)值解。在這種情況下，其非負(fù)最小二乘解比方程的精 確解更有意義。在MATLAB中，求非負(fù)最小二乘解常用函數(shù) nnls,其調(diào)用格式 為：X=nnls(A,b)返回方程Ax=b的最小二乘解，方程的求解過程被限制在 x的 條件下；文

8、檔來源為：從網(wǎng)絡(luò)收集整理,word版本可編輯.歡迎下載支持X=nnls(A,b,TOL)指定誤差TOL來求解，TOL的默認(rèn)值為TOL=max(size(A)*norm(A,1)*eps ,矩陣的一1范數(shù)越大，求解的誤差越大；X,W=nnls(A,b)當(dāng) x(i)=0 時，w(i)0 時，w(i)0,同時返回一 個雙向量Wo【例9】求方程組的非負(fù)最小二乘解A=3.4336 -0.5238 0.6710-0.5238 3.2833 -0.73020.6710 -0.7302 4.0261;b=-1.000 1.5000 2.5000;X,W=nnls(A,b)X=00.65630.6998W=-3

9、.6820-0.0000-0.0000 x1=Abx1 =-0.35690.57440.7846A*X-bans=1.12580.1437-0.1616A*x1-b ans=文檔來源為：從網(wǎng)絡(luò)收集整理,word版本可編輯.歡迎下載支持1.0e-0.15-0.22200.44410關(guān)于采用matlab進(jìn)行指定非線性方程擬合的問題(1)派1。優(yōu)化工具箱的利用函數(shù) 描述LSQLIN有約束線性最小二乘優(yōu)化LSQNONNEG 非負(fù)約束線性最小二乘優(yōu)化問題當(dāng)有約束問題存在的時候，應(yīng)該采用上面的方法代替Polyfit與反斜線()。具體例子請參閱優(yōu)化工具箱文檔中的相應(yīng)利用這兩個函數(shù)的例子。d,非線性曲線擬合利

10、用MATLAB 的內(nèi)建函數(shù)函數(shù)名描述FMINBND只解決單變量固定區(qū)域的最小值問題FMINSEARCH 多變量無約束非線性最小化問題( Nelder-Mead方法)。下面給出一個小例子展示一下如何利用FMINSEARCH1.首先生成數(shù)據(jù) t=0:.1:10； t=t(:)； Data=40*exp(-.5*t)+rand(size(t); % 將數(shù)據(jù)加上隨機(jī)噪聲.寫一個m文件，以曲線參數(shù)作為輸入，以擬合誤差作為輸出function sse=myfit(params,Input,Actural_Output)A=params(1);lamda=params(2);Fitted_Curve=A.*

11、exp(-lamda*Input);Error_Vector=Fitted_Curve-Actural_Output;文檔來源為：從網(wǎng)絡(luò)收集整理,word版本可編輯.歡迎下載支持%當(dāng)曲線擬合的時候，一個典型的質(zhì)量評價標(biāo)準(zhǔn)就是誤差平方和sse=sum(Error_Vector.A2);%當(dāng)然，也可以將 sse寫作：sse=Error_Vector(:)*Error_Vector(:);.調(diào)用 FMINSEARCH Strarting=rand(1,2); options=optimset(Display,iter); Estimates=fiminsearch(myfit,Strarting,o

12、ptions,t,Data); plot(t,Data,*); hold on plot(t,Estimates(1)*exp(-Estimates(2)*t),r);Estimates將是一個包含了對原數(shù)據(jù)集進(jìn)行估計的參數(shù)值的向量。附圖見后面：FMINSEARCH通常能夠用來解決不連續(xù)情況，特別是如果他們不出現(xiàn)在解的附近的時候。它得到的通常也是局部解。FMINSEARCH只能夠最小化實數(shù)值(也就是說，解的域必須只能包括實數(shù)，函數(shù)的輸出只能夠為實數(shù)值)。當(dāng)感興趣的是復(fù)數(shù)變量的域的時候，他們必須被分割為實部與虛部。X2.MATLAB 的FIGURE窗口：最基本的擬合界面與數(shù)據(jù)統(tǒng)計工具M(jìn)ATLAB

13、通過基本的擬合界面也支持基本曲線擬合。利用這個界面，你可以快速地在簡單 易用的環(huán)境中實現(xiàn)許多基本的曲線擬合。這個界面可以實現(xiàn)以下功能：a.通過比樣條插值(spline interpolant)、hermite插值、或者是高達(dá)10階的多項式插值實文檔來源為：從網(wǎng)絡(luò)收集整理,word版本可編輯.歡迎下載支持現(xiàn)數(shù)據(jù)的擬合；b.對給定數(shù)據(jù)同時實現(xiàn)多樣插值的繪制；c.繪制殘差圖；d.檢查擬合結(jié)果的殘差的數(shù)值；e.通過內(nèi)插值或者外推插值評價一個擬合結(jié)果；f .對擬合結(jié)果和殘差的模進(jìn)行圖形繪制；g.將擬合結(jié)果保存入 MATLAB工作空間。開發(fā)你的擬合應(yīng)用的時候，你可以通過基本擬合( Basic Fittin

14、g )界面，也可以通過命令行函數(shù)，也可以同時作用。你可以通過基本擬合界面只能夠?qū)崿F(xiàn)2-D數(shù)據(jù)的擬合。然而，如果你用subplot繪制多個數(shù)據(jù)集，只要有至少一個數(shù)據(jù)集是2D的，那么就可以用基本擬合界面?？梢酝ㄟ^如下步驟激活基本擬合界面：繪制數(shù)據(jù)；從figure窗口的 Tools菜單條下面選擇 Basic Fitting 菜單項；有關(guān)Basic Fitting界面的更多信息，請查閱 MATLAB幫助文檔的相應(yīng)部分。注意：對于 HP, IBM 以及 SGI 平臺，MATLAB6.0 (R12.0)以及 MATLAB6.1(R12.1)的基 本擬合界面不受支持。數(shù)據(jù)統(tǒng)計界面可以用來對圖形中的每個數(shù)據(jù)集

15、進(jìn)行統(tǒng)計量的計算?？梢酝ㄟ^如下步驟將數(shù)據(jù)統(tǒng)計界面激活： 1.制數(shù)據(jù)；從figure窗口的 Tools菜單條下面選擇 Data Statistics菜單項；關(guān)于采用matlab進(jìn)行指定非線性方程擬合的問題(2)一。優(yōu)化工具箱函數(shù)LSQNONLIN解決非線性最小二乘法問題，包括非線性數(shù)據(jù)擬合問題LSQCURVEFIT 解決非線性數(shù)據(jù)擬合問題下面給出利用這兩個函數(shù)的例子：LSQNONLIN :利用這個函數(shù)最小化連續(xù)函數(shù)只能夠找到句柄解。下面的例子說明利用LSQNONLIN函數(shù)用下面的函數(shù)進(jìn)行擬合：9文檔來源為：從網(wǎng)絡(luò)收集整理,word版本可編輯.歡迎下載支持f = A + B exp(C*x)+D*

16、exp(E*x)對數(shù)據(jù)集x與y進(jìn)行擬合，其中y是在給定x的情況下的期望輸出(可以是方程給出數(shù)組，也可以是單獨數(shù)據(jù)組成的數(shù)組)。為了解決這個問題，先建立下面的名為fit_simp.m的函數(shù)，它利用數(shù)據(jù)x與y,將他們作為優(yōu)化輸入?yún)?shù)傳遞給LSQNONLIN。利用給定的數(shù)據(jù)x計算f的值，再與原始數(shù)據(jù)y進(jìn)行比較。經(jīng)驗值與實際計算出的值之間的差異作為輸出值返回。LSQNOLIN函數(shù)就是最小化這些差的平方和。function diff = fit_simp(x,X,Y)%此函數(shù)被LSQNONLIN調(diào)用% x是包含等式系數(shù)的向量% X與Y是作為操作數(shù)傳遞給IsnonlinA = x(1);B = x(2);

17、C = x(3);D = x(4);E = x(5);diff = A + B.*exp(C.*X)+D.*exp(E.*X)-Y;下面的腳本是利用上面定義的fit_simp.m函數(shù)的一個例子：%定義你打算擬合的數(shù)據(jù)集合 X=0:.01:.5; Y=2.0.*exp(5.0.*X)+3.0.*exp(2.5.*X)+1.5.*rand(size(X);%初始化方程系數(shù) X0=1 1 1 1 1;% 設(shè)置用中等模式(memdium-scale )算法 options=optimset(LargescaIe,off);%通過調(diào)用LSQNONLIN重現(xiàn)計算新的系數(shù) x=IsqnonIin(fit_s

18、imp,X0,options,X,Y);%調(diào)用LSQNONLIN結(jié)果輸出表明擬合是成功的Optimization terminated successfully:Gradient in the search direction less than tolFunGradient less than 10*(tolFun+tolX)10文檔來源為：從網(wǎng)絡(luò)收集整理,word版本可編輯.歡迎下載支持%繪制原始數(shù)據(jù)與新的計算的數(shù)據(jù) Y_new=x(1)+x(2).*exp(x(3).*X)+x(4).*exp(x(5).*X); plot(X,Y,+r,X,Y_new,b);注意：LSQNONLIN 只

19、可以處理實數(shù)變量。在處理包括復(fù)數(shù)變量的實例的擬合的時候，數(shù)據(jù)集應(yīng)該被切分成實數(shù)與虛數(shù)部分。下面給出一個例子演示如何對復(fù)數(shù)參數(shù)進(jìn)行最小二乘 擬合。為了擬合復(fù)數(shù)變量，你需要將復(fù)數(shù)分解為實數(shù)部分與虛數(shù)部分，然后把他們傳遞到函數(shù)中去，這個函數(shù)被LEASTSQ作為單個輸入調(diào)用。首先，將復(fù)數(shù)分解為實部與虛部兩個向量。其次,將這兩個向量理解成諸如第一部分是實部、第二部分是虛部。在 MATLAB函數(shù)中，重新裝 配復(fù)數(shù)數(shù)據(jù)，并用你想擬合的復(fù)數(shù)方程計算。將輸出向量分解實部與虛部，將這兩部分連接 為一個單一的輸出向量傳遞回LEASTSQ。下面，給出一個例子演示如何根據(jù)兩個復(fù)數(shù)指數(shù)擬合實數(shù)X與丫。建立方程：funct

20、ion zero = fit2(x,X,Y)%根據(jù)輸入x重建復(fù)數(shù)輸入 cmpx = x(1:4)+i.*x(5:8);%利用復(fù)數(shù)計算函數(shù)zerocomp = cmpx(1),*exp(cmpx(2).*X) + cmpx(3),* exp(cmpx(4).*X)-Y;%將結(jié)果轉(zhuǎn)換成一個列向量%其中第一部分是實部，第二部分是虛部numx = length(X); % 實部長度zero=real(zerocomp); % 實部zero(numx+1:2*numx)=imag(zerocomp); % 虛部為了評價計算這個函數(shù)，需要X與Y數(shù)據(jù)集。LSQNONLIN將根據(jù)它擬合出下面方程中的參數(shù)a,

21、b, c與d：Y = a*exp(b*X)+c*exp(d*X);其中，a, b, c與d是復(fù)數(shù)變量。 X=0:.1:5; Y=sin(X); Y=Y+.1*rand(size(Y)-.05;11文檔來源為：從網(wǎng)絡(luò)收集整理,word版本可編輯.歡迎下載支持 cmpx0=1 i 2 2*i; x0(1:4)=real(cmpx0); x0(5:8)=imag(cmpx0); x=leastsq(fit2,x0,X,Y); cmpx=x(1:4)+i.*x(5:8); Y1=real(cmpx(1).*exp(cmpx(2).*X)+cmpx(3).*exp(cmpx(4).*X); plot(X

22、,Y1,r); hold on plot(X,Y,+);0 O.S 11522 5335 A 4.56二。LSQCURVEFIT :利用此函數(shù)可以在最小二乘意義上解決非線性曲線擬合(數(shù)據(jù)擬合) 問題。也就是說，給定輸入數(shù)據(jù)xdata,以及觀測的輸出數(shù)據(jù) ydata,找到系數(shù)x,使得函數(shù)F(x, xdata)能夠最好的擬合向量值。LSQCURVEFIT 利用與LSQNONLIN 相同的算法。它的目的在于專門為數(shù)據(jù)擬合問題提供一個接口。在擬合的時候，2維、3維或者N維參數(shù)擬合是沒有什么差別的。下面給出一個3維參數(shù)擬合的例子。待擬合函數(shù)是：z = a1*y.*x.A2+a2*sin(x)+a3*y.

23、A3;12文檔來源為：從網(wǎng)絡(luò)收集整理,word版本可編輯.歡迎下載支持建立的myfun.m的函數(shù)如下：function F = myfun(a, data);x = data(1,:);y = data(2,:);F= a(1)*y.*x.A2+a(2)*sin(x)+a(3)*y.A3;下面的腳本展示了這么利用上面的函數(shù)： xdata= 3.6 7.7 9.3 4.1 8.6 2.8 1.3 7.9 10.0 5.4; ydata= 16,5 150.6 263.1 24.7 208.5 9.9 2.7 163.9 325.0 54.3; zdata= 95,09 23.11 60.63 4

24、8.59 89.12 76.97 45.68 1,84 82.17 44.47; data=xdata; ydata; a0= 10,10,10; % 初識揣測 a, resnorm = lsqcurvefit(myfun,a0,data,zdata)Maximum number of function evaluations exceeded;increase options.MaxFunEvalsa = 0.0088 -34.2886 -0.0000resnorm = 2.2636e+004 format long aa = 0.007493 -34.919983 -0.000009 option=optimset(MaxFunEvals,800); a, resnorm = lsqcurvefit(myfun,a0, data, zdata, , , option)Optimization terminated successfully:Relative function value changing by less than OPTIONS.TolFuna = 0.009653 -20.111138 -0.00000