matlab符號計算及其應用_第1頁
matlab符號計算及其應用_第2頁
matlab符號計算及其應用_第3頁
matlab符號計算及其應用_第4頁
matlab符號計算及其應用_第5頁
已閱讀5頁,還剩93頁未讀 繼續(xù)免費閱讀

下載本文檔

版權(quán)說明:本文檔由用戶提供并上傳,收益歸屬內(nèi)容提供方,若內(nèi)容存在侵權(quán),請進行舉報或認領(lǐng)

文檔簡介

MATLAB

軟件及其應用

ApplicationofMatlabLanguage1

MATLAB符號運算

(Symbolic)2Matlab

符號運算介紹

Matlab符號運算是通過符號數(shù)學工具箱(SymbolicMathToolbox)來實現(xiàn)的。Matlab符號數(shù)學工具箱是建立在功能強大的Maple軟件的基礎(chǔ)上的,當Matlab進行符號運算時,它就請求Maple軟件去計算并將結(jié)果返回給Matlab。

Matlab的符號數(shù)學工具箱可以完成幾乎所有得符號運算功能。主要包括:符號表達式的運算,符號表達式的復合、化簡,符號矩陣的運算,符號微積分、符號作圖,符號代數(shù)方程求解,符號微分方程求解等。此外,該工具箱還支持可變精度運算,即支持以指定的精度返回結(jié)果。3Matlab

符號運算特點計算以推理方式進行,因此不受計算誤差累積所帶來的困擾。

符號計算指令的調(diào)用比較簡單,與數(shù)學教科書上的公式相近。

符號計算可以給出完全正確的封閉解,或任意精度的數(shù)值解(封閉解不存在時)。

符號計算所需的運行時間相對較長。4Matlab

符號運算舉例求一元二次方程ax2+bx+c=0

的根>>

solve('a*x^2+b*x+c=0')求的根f(x)=(cos

x)2

的一次導數(shù)>>

x=sym('x');>>

diff(cos(x)^2)計算f(x)=x2

在區(qū)間[a,b]

上的定積分>>

symsabx;>>

int(x^2,a,b)5在進行符號運算時,必須先定義基本的符號對象,可以是符號常量、符號變量、符號表達式等。符號對象是一種數(shù)據(jù)結(jié)構(gòu)。符號對象與符號表達式含有符號對象的表達式稱為符號表達式,Matlab在內(nèi)部把符號表達式表示成字符串,以與數(shù)字變量或運算相區(qū)別。符號矩陣/數(shù)組:元素為符號表達式的矩陣/數(shù)組。6

sym

函數(shù)用來建立單個符號變量,一般調(diào)用格式為:符號對象的建立:sym

和syms符號對象的建立例:>>

a=sym('a')符號變量

=

sym(A)參數(shù)A

可以是一個數(shù)或數(shù)值矩陣,也可以是字符串a(chǎn)

是符號變量b

是符號常量>>

b=sym(1/3)C

是符號矩陣>>

C=sym('[1ab;cd]')7符號對象的建立:sym

和syms符號對象的建立

syms

命令用來建立多個符號變量,一般調(diào)用格式為:syms符號變量1符號變量2...符號變量n例:>>

symsabc>>

a=sym('a');>>

b=sym('b');>>

c=sym('c');8符號表達式的建立:例:建立符號表達式通常有以下2種方法:(1)用sym

函數(shù)直接建立符號表達式。

(2)使用已經(jīng)定義的符號變量組成符號表達式。>>

y=sym('sin(x)+cos(x)')>>x=sym('x');>>

y=sin(x)+cos(x)符號表達式的建立>>symsx;>>

y=sin(x)+cos(x)9Matlab符號運算采用的運算符和基本函數(shù),在形狀、名稱和使用上,都與數(shù)值計算中的運算符和基本函數(shù)完全相同符號對象的基本運算基本運算符普通運算:+

-

*

\

/

^數(shù)組運算:.*

.\

./

.^矩陣轉(zhuǎn)置:'

.'例:>>X=sym('[x11,x12;x21,x22;x31,x32]');>>

Y=sym('[y11,y12,y13;y21,y22,y23]');>>

Z1=X*Y;Z2=X'.*Y;10符號對象的基本運算sin、cos、tan、cot、sec、csc、…asin、acos、atan、acot、asec、acsc、…exp、log、log2、log10、sqrtabs、conj、real、imagrank、det、inv、eig、lu、qr、svddiag、triu、tril、expm三角函數(shù)與反三角函數(shù)、指數(shù)函數(shù)、對數(shù)函數(shù)等基本函數(shù)11查找符號表達式中的符號變量若表達式中有兩個符號變量與x

的距離相等,則ASCII碼大者優(yōu)先。查找符號變量findsym(expr)

按字母順序列出符號表達式expr

中的所有符號變量findsym(expr,N)

按順序列出expr

中離x

最近的N

個符號變量常量pi,

i,

j

不作為符號變量12例:>>

f=sym('2*w-3*y+z^2+5*a')>>

findsym(f)>>

findsym(f,3)>>

findsym(f,1)findsym

舉例13符號表達式的替換subs(f,x,a)用a

替換字符函數(shù)

f

中的字符變量x

a

是可以是數(shù)/數(shù)值變量/表達式或字符變量/表達式若x

是一個由多個字符變量組成的數(shù)組或矩陣,則a

應該具有與x

相同的形狀的數(shù)組或矩陣。

用給定的數(shù)據(jù)替換符號表達式中的指定的符號變量14subs舉例>>

f=sym('2*u');>>

subs(f,'u',2)>>

f2=subs(f,'u','u+2')>>

a=3;>>

subs(f2,'u',a+2)>>

subs(f2,'u','a+2')>>

symsxy>>

f3=subs(f,'u',x+y)>>

subs(f3,[x,y],[1,2])ans=4f2=2*(u+2)ans=14ans=2*((a+2)+2)f3=2*x+2*yans=6例:指出下面各條語句的輸出結(jié)果f=2*u15符號矩陣>>

A=sym('[1+x,sin(x);5,exp(x)]')

使用sym

函數(shù)直接生成

將數(shù)值矩陣轉(zhuǎn)化成符號矩陣

符號矩陣中元素的引用和修改>>

B=[2/3,sqrt(2);5.2,log(3)];>>

C=sym(B)>>

A=sym('[1+x,sin(x);5,exp(x)]');>>

A(1,2)%引用>>

A(2,2)=sym('cos(x)')%重新賦值16六類常見符號運算因式分解、展開、合并、簡化及通分等計算極限計算導數(shù)計算積分符號求和代數(shù)方程和微分方程求解17因式分解因式分解factor(f)>>

symsx;f=x^6+1;>>

factor(f)

factor

也可用于正整數(shù)的分解>>

s=factor(100)>>

factor(sym('12345678901234567890'))大整數(shù)的分解要轉(zhuǎn)化成符號常量18函數(shù)展開函數(shù)展開expand(f)>>

symsx;f=(x+1)^6;>>

expand(f)

多項式展開

三角函數(shù)展開>>

symsxy;f=sin(x+y);>>

expand(f)19合并同類項合并同類項collect(f,v):按指定變量v

進行合并collect(f):按默認變量進行合并>>

symsxy;

>>

f=x^2*y+y*x-x^2+2*x;>>

collect(f)>>

collect(f,y)20函數(shù)簡化函數(shù)簡化y=simple(f):

對f

嘗試多種不同的算法進行簡化,返回其中最簡短的形式[How,y]=simple(f):

y

為f

的最簡短形式,How中記錄的為簡化過程中使用的方法。fRHOW2*cos(x)^2-sin(x)^23*cos(x)^2-1simplify(x+1)*x*(x-1)x^3-xcombine(trig)x^3+3*x^2+3*x+1(x+1)^3factorcos(3*acos(x))4*x^3-3*xexpand21函數(shù)簡化函數(shù)簡化y=simplify(f):

對f

進行簡化>>

symsx;f=sin(x)^2+cos(x)^2;>>

simplify(f)>>

symscalphabeta;

>>

f=exp(c*log(sqrt(alpha+beta)));>>

simplify(f)22函數(shù)簡化舉例>>

symsx;

>>

f=(1/x^3+6/x^2+12/x+8)^(1/3);>>

y1=simplify(f)>>

g1=simple(f)>>

g2=simple(g1)多次使用simple

可以達到最簡表達。例:簡化23分式通分函數(shù)簡化[N,D]=numden(f):

N

為通分后的分子,D

為通分后的分母>>

symsxy;

>>

f=x/y+y/x;>>

[N,D]=numden(f)>>

[n,d]=numden(sym(112/1024))24horner

多項式

horner

多項式:嵌套形式的多項式>>

symsx;

>>

f=x^4+2*x^3+4*x^2+x+1;>>

g=horner(f)例:25計算極限limit(f,x,a):計算limit(f,a):當默認變量趨向于a時的極限limit(f):計算a=0時的極限limit(f,x,a,'right'):計算右極限limit(f,x,a,'left'):計算左極限例:計算,>>

symsxhn;

>>

L=limit((log(x+h)-log(x))/h,h,0)>>

M=limit((1-x/n)^n,n,inf)26計算導數(shù)g=diff(f,v):求符號表達式

f

關(guān)于

v

的導數(shù)g=diff(f):求符號表達式

f

關(guān)于默認變量的導數(shù)g=diff(f,v,n):求

f

關(guān)于

v

n階導數(shù)

diff>>

symsx;>>

f=sin(x)+3*x^2;>>

g=diff(f,x)27計算積分int(f,v,a,b):計算定積分int(f,a,b):計算關(guān)于默認變量的定積分int(f,v):計算不定積分int(f):計算關(guān)于默認變量的不定積分>>

symsx;f=(x^2+1)/(x^2-2*x+2)^2;>>

I=int(f,x)>>

K=int(exp(-x^2),x,0,inf)例:計算和28符號求和>>

symsn;f=1/n^2;>>

S=symsum(f,n,1,inf)>>

S100=symsum(f,n,1,100)symsum(f,v,a,b):求和symsum(f,a,b):關(guān)于默認變量求和

例:計算級數(shù)及其前100項的部分和例:計算函數(shù)級數(shù)>>

symsnx;f=x/n^2;>>

S=symsum(f,n,1,inf)29代數(shù)方程求解solve(f,v):求方程關(guān)于指定自變量的解,f

可以是用字符串表示的方程、符號表達式或符號方程;

solve也可解方程組(包含非線性);得不到解析解時,給出數(shù)值解。30微分方程求解

dsolvey=dsolve('eq1','eq2',...,'cond1','cond2',...,'v')其中

y

為輸出的解,

eq1、eq2、...為微分方程,

cond1、cond2、...為初值條件,v

為自變量例1:求微分方程的通解,并驗證。>>

y=dsolve('Dy+2*x*y=x*exp(-x^2)','x')>>

symsx;diff(y)+2*x*y-x*exp(-x^2)31微分方程求解

幾點說明如果省略初值條件,則表示求通解;如果省略自變量,則默認自變量為t

dsolve('Dy=2*x','x');%

dy/dx=2xdsolve('Dy=2*x');%dy/dt=2x若找不到解析解,則返回其積分形式。微分方程中用D

表示對自變量的導數(shù),如:Dyy';

D2yy'';

D3yy'''32微分方程求解例2:求微分方程滿足初值條件

的特解,并畫出解函數(shù)的圖形。>>

y=dsolve('x*Dy+y-exp(x)=0',...

'y(1)=2*exp(1)','x')>>

ezplot(y);33微分方程求解例3:求微分方程組在初值條件下的特解,并畫出解函數(shù)的圖形。[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0',...'x(0)=1','y(0)=0','t')ezplot(x,y,[0,1.3]);注:解微分方程組時,如果所給的輸出個數(shù)與方程個數(shù)相同,則方程組的解按詞典順序輸出。34其它運算反函數(shù)finverse(f,v):求f

關(guān)于指定變量v

的反函數(shù)finverse(f):求f

關(guān)于默認變量的反函數(shù)>>

symsxt;f=x^2+2*t;>>

g1=finverse(f,x)>>

g2=finverse(f,t)例:計算函數(shù)的反函數(shù)35

MATLAB應用--解線性規(guī)劃

36用MATLAB優(yōu)化工具箱解線性規(guī)劃minz=cX

1、模型:命令:x=linprog(c,A,b)

2、模型:minz=cX

命令:x=linprog(c,A,b,Aeq,beq)注意:若沒有不等式:存在,則令A=[],b=[].373、模型:minz=cX

VLB≤X≤VUB命令:[1]x=linprog(c,A,b,Aeq,beq,VLB,VUB)

[2]x=linprog(c,A,b,Aeq,beq,VLB,VUB,X0)

注意:[1]若沒有等式約束:,則令Aeq=[],beq=[].[2]其中X0表示初始點4、命令:[x,fval]=linprog(…)返回最優(yōu)解x及x處的目標函數(shù)值fval.38解編寫M文件xxgh1.m如下:c=[-0.4-0.28-0.32-0.72-0.64-0.6];A=[0.010.010.010.030.030.03;0.02000.0500;00.02000.050;000.03000.08];b=[850;700;100;900];

Aeq=[];beq=[];

vlb=[0;0;0;0;0;0];vub=[];[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)

39解:編寫M文件xxgh2.m如下:

c=[634];A=[010];b=[50];

Aeq=[111];

beq=[120];

vlb=[30,0,20];

vub=[];[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)40

MATLAB應用—

解非線性規(guī)劃

41一、二次規(guī)劃(QuadraticProgram)

概念42二、Matlab中求解二次規(guī)劃4344轉(zhuǎn)化為matlab求解格式:4546定義

如果目標函數(shù)或約束條件中至少有一個是非線性函數(shù),則最優(yōu)化問題就叫做非線性規(guī)劃問題.四、非線性規(guī)劃的基本概念一般形式:

(1)其中,是定義在Rn

上的實值函數(shù)()n

TnRxxxX?=,,,21L()()???íì===3.,...,2,1

0

m;1,2,...,

0..

ljXhiXgtsji47Matlab求解非線性規(guī)劃問題

其中X為n維變元向量,G(X)與Ceq(X)均為非線性函數(shù)組成的向量。48

1.首先建立M文件fun.m,用來定義目標函數(shù)F(X):functionf=fun(X);f=F(X);MATLAB求解上述問題,基本步驟分三步493.建立主程序.求解非線性規(guī)劃的函數(shù)是fmincon,命令的基本格式如下:

(1)x=fmincon(‘fun’,X0,A,b)

(2)x=fmincon(‘fun’,X0,A,b,Aeq,beq)

(3)x=fmincon(‘fun’,X0,A,b,Aeq,beq,VLB,VUB)

(4)x=fmincon(‘fun’,X0,A,b,Aeq,beq,VLB,VUB,’nonlcon’)(5)x=fmincon(‘fun’,X0,A,b,Aeq,beq,VLB,VUB,’nonlcon’,options)

(6)[x,fval]=fmincon(…)

(7)[x,fval,exitflag]=fmincon(…)(8)[x,fval,exitflag,output]=fmincon(…)輸出極值點M文件迭代的初值參數(shù)說明變量上下限fmincon函數(shù)可能會給出局部最優(yōu)解,這與初值X0的選取有關(guān).501.寫成標準形式:

s.t.

2x1+3x26

s.t.

x1+4x25

x1,x20例512.先建立M-文件fun3.m:

functionf=fun3(x);f=-x(1)-2*x(2)+(1/2)*x(1)^2+(1/2)*x(2)^23.再建立主程序youh2.m:

x0=[1;1];A=[23;14];b=[6;5];

Aeq=[];beq=[];VLB=[0;0];VUB=[];[x,fval]=fmincon('fun3',x0,A,b,Aeq,beq,VLB,VUB)4.運算結(jié)果為:

x=0.76471.0588fval=-2.0294521.先建立M文件fun4.m定義目標函數(shù):

functionf=fun4(x);f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);

x1+x2=0

s.t.1.5+x1x2-x1-x20

-x1x2–10

0例

2.再建立M文件mycon.m定義非線性約束:

function[g,ceq]=mycon(x)g=[x(1)+x(2);1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10];533.主程序youh3.m為:x0=[-1;1];A=[];b=[];Aeq=[11];beq=[0];vlb=[];vub=[];[x,fval]=fmincon('fun4',x0,A,b,Aeq,beq,vlb,

vub,'mycon')4.運算結(jié)果為:

x=-1.22501.2250fval=1.895154例1.先建立M文件fun.m定義目標函數(shù):functionf=fun(x);f=-2*x(1)-x(2);2.再建立M文件mycon2.m定義非線性約束:function[g,ceq]=mycon2(x)g=[x(1)^2+x(2)^2-25;x(1)^2-x(2)^2-7];

553.主程序fxx.m為:

x0=[3;2.5];VLB=[00];VUB=[510];[x,fval,exitflag,output]=fmincon('fun',x0,[],[],[],[],VLB,VUB,'mycon2')56

MATLAB應用—

微分方程解法

57要求目的主要內(nèi)容2、學會用Matlab求微分方程的數(shù)值解.1、學會用Matlab求簡單微分方程的解析解.1、求簡單微分方程的解析解.2、求微分方程的數(shù)值解.58求微分方程的數(shù)值解(一)常微分方程數(shù)值解的定義(二)建立數(shù)值解法的一些途徑(三)用Matlab軟件求常微分方程的數(shù)值解返回59微分方程的解析解求微分方程(組)的解析解命令:dsolve(‘方程1’,‘方程2’,…‘方程n’,‘初始條件’,‘自變量’)結(jié)果:u=tg(t-c)60解輸入命令:y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')結(jié)果為:y=3e-2xsin(5x)61解輸入命令:

[x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z','t');

x=simple(x)%將x化簡

y=simple(y)z=simple(z)結(jié)果為:x=(c1-c2+c3+c2e-3t-c3e-3t)e2t

y=-c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2tz=(-c1e-4t+c2e-4t+c1-c2+c3)e2t

62微分方程的數(shù)值解(一)常微分方程數(shù)值解的定義

在生產(chǎn)和科研中所處理的微分方程往往很復雜且大多得不出一般解。而在實際上對初值問題,一般是要求得到解在若干個點上滿足規(guī)定精確度的近似值,或者得到一個滿足精確度要求的便于計算的表達式。因此,研究常微分方程的數(shù)值解法是十分必要的。63(二)建立數(shù)值解法的一些途徑1、用差商代替導數(shù)

若步長h較小,則有故有公式:此即歐拉法。642、使用數(shù)值積分對方程y’=f(x,y),兩邊由xi到xi+1積分,并利用梯形公式,有:實際應用時,與歐拉公式結(jié)合使用:此即改進的歐拉法。故有公式:653、使用泰勒公式

以此方法為基礎(chǔ),有龍格-庫塔法、線性多步法等方法。4、數(shù)值公式的精度當一個數(shù)值公式的截斷誤差可表示為O(hk+1)時(k為正整數(shù),h為步長),稱它是一個k階公式。k越大,則數(shù)值公式的精度越高。歐拉法是一階公式,改進的歐拉法是二階公式。龍格-庫塔法有二階公式和四階公式。線性多步法有四階阿達姆斯外插公式和內(nèi)插公式。66(三)用Matlab軟件求常微分方程的數(shù)值解[t,x]=solver(’f’,ts,x0,options)ode45ode23ode113ode15sode23s由待解方程寫成的m-文件名ts=[t0,tf],t0、tf為自變量的初值和終值函數(shù)的初值ode23:組合的2/3階龍格-庫塔-芬爾格算法ode45:運用組合的4/5階龍格-庫塔-芬爾格算法自變量值函數(shù)值用于設定誤差限(缺省時設定相對誤差10-3,絕對誤差10-6),命令為:options=odeset(’reltol’,rt,’abstol’,at),rt,at:分別為設定的相對誤差和絕對誤差.67

1、在解n個未知函數(shù)的方程組時,x0和x均為n維向量,m-文件中的待解方程組應以x的分量形式寫成.

2、使用Matlab軟件求數(shù)值解時,高階微分方程必須等價地變換成一階微分方程組.注意:68解:令y1=x,y2=y1’1、建立m-文件vdp1000.m如下:

function

dy=vdp1000(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1000*(1-y(1)^2)*y(2)-y(1);

2、取t0=0,tf=3000,輸入命令:

[T,Y]=ode15s('vdp1000',[03000],[20]);plot(T,Y(:,1),'-')3、結(jié)果如圖69解

1、建立m-文件rigid.m如下:

function

dy=rigid(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);2、取t0=0,tf=12,輸入命令:

[T,Y]=ode45('rigid',[012],[011]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')3、結(jié)果如圖圖中,y1的圖形為實線,y2的圖形為“*”線,y3的圖形為“+”線.70Matlab的應用

-----------插值和擬合71拉格朗日插值分段線性插值三次樣條插值一維插值一、插值的定義二、插值的方法三、用Matlab解插值問題72一維插值的定義已知n+1個節(jié)點其中互不相同,不妨設求任一插值點處的插值節(jié)點可視為由產(chǎn)生,,表達式復雜,,或無封閉形式或未知73構(gòu)造一個(相對簡單的)函數(shù)通過全部節(jié)點,即再用計算插值,即74用MATLAB作插值計算一維插值函數(shù):yi=interp1(x,y,xi,'method')插值方法被插值點插值節(jié)點xi處的插值結(jié)果‘nearest’

:最鄰近插值‘linear’

:線性插值;‘spline’

:三次樣條插值;‘cubic’

:立方插值。缺省時:分段線性插值。注意:所有的插值方法都要求x是單調(diào)的,并且xi不能夠超過x的范圍。75例:在1-12的11小時內(nèi),每隔1小時測量一次溫度,測得的溫度依次為:5,8,9,15,25,29,31,30,22,25,27,24。試估計每隔1/10小時的溫度值。hours=1:12;temps=[589152529313022252724];h=1:0.1:12;t=interp1(hours,temps,h,'spline');(直接輸出數(shù)據(jù)將是很多的)plot(hours,temps,'+',h,t,hours,temps,'r:')%作圖xlabel('Hour'),ylabel('DegreesCelsius’)76擬合2.擬合的基本原理1.擬合問題引例77擬合問題引例1溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032已知熱敏電阻數(shù)據(jù):求600C時的電阻R。

R=at+ba,b為待定系數(shù)78擬合問題引例2

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01已知一室模型快速靜脈注射下的血藥濃度數(shù)據(jù)(t=0注射300mg)求血藥濃度隨時間的變化規(guī)律c(t).作半對數(shù)坐標系(semilogy)下的圖形79曲線擬合問題的提法已知一組(二維)數(shù)據(jù),即平面上n個點(xi,yi)i=1,…n,

尋求一個函數(shù)(曲線)y=f(x),

使f(x)

在某種準則下與所有數(shù)據(jù)點最為接近,即曲線擬合得最好。

+++++++++xyy=f(x)(xi,yi)ii為點(xi,yi)與曲線y=f(x)的距離80擬合與插值的關(guān)系函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個函數(shù)作為近似,由于近似的要求不同,二者的數(shù)學方法上完全不同。實例:下面數(shù)據(jù)是某次實驗所得,希望得到x和

f之間的關(guān)系?問題:給定一批數(shù)據(jù)點,需確定滿足特定要求的曲線或曲面解決方案:若不要求曲線(面)通過所有數(shù)據(jù)點,而是要求它反映對象整體的變化趨勢,就是數(shù)據(jù)擬合,又稱曲線擬合或曲面擬合。若要求所求曲線(面)通過所給所有數(shù)據(jù)點,就是插值問題;81曲線擬合問題最常用的解法——線性最小二乘法的基本思路第一步:先選定一組函數(shù)

r1(x),r2(x),…rm(x),m<n,

f(x)=a1r1(x)+a2r2(x)+…+amrm(x)(1)其中

a1,a2,…am

為待定系數(shù)。

第二步:確定a1,a2,…am

的準則(最小二乘準則):使n個點(xi,yi)與曲線y=f(x)的距離i的平方和最小

。記

問題歸結(jié)為,求

a1,a2,…am

使

J(a1,a2,…am)

最小。82線性最小二乘法的求解:預備知識超定方程組:方程個數(shù)大于未知量個數(shù)的方程組即Ra=y其中超定方程一般是不存在解的矛盾方程組。如果有向量a使得達到最小,則稱a為上述超定方程的最小二乘解。

83線性最小二乘法的求解定理:當RTR可逆時,超定方程組(3)存在最小二乘解,

且即為方程組

RTRa=RTy------正則(正規(guī))方程組的解:a=(RTR)-1RTy所以,曲線擬合的最小二乘法要解決的問題,實際上就是求以下超定方程組的最小二乘解的問題。其中Ra=y(3)84用MATLAB解擬合問題1、線性最小二乘擬合2、非線性最小二乘擬合85用MATLAB作線性最小二乘擬合1.作多項式f(x)=a1xm+…+amx+am+1擬合,可利用已有程序:a=polyfit(x,y,m)2.對超定方程組可得最小二乘意義下的解。,用3.多項式在x處的值y的計算命令:y=polyval(a,x)輸出擬合多項式系數(shù)a=[a1,…,am,am+1]’

(數(shù)組)輸入同長度數(shù)組X,Y擬合多項式

次數(shù)86即要求出二次多項式:中的使得:例對下面一組數(shù)據(jù)作二次多項式擬合871)輸入命令:x=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];R=[(x.^2)',x',ones(11,1)];

A=R\y'解法1.解超定方程的方法2)計算結(jié)果:A=[-9.8108,20.1293,-0.0317]882)計算結(jié)果:A=[-9.8108,20.1293,-0.0317]解法2.用多項式擬合的命令MATLAB(zxec2)1)輸入命令:x=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];A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,'k+',x,z,'r')%作出數(shù)據(jù)點和擬合曲線的圖形891.lsqcurvefit已知數(shù)據(jù)點:xdata=(xdata1,xdata2,…,xdatan)

ydata=(ydata1,ydata2,…,ydatan)用MATLAB作非線性最小二乘擬合兩個求非線性最小二乘擬合的函數(shù):lsqcurvefit、lsqnonlin。相同點和不同點:兩個命令都要先建立M-文件fun.m,定義函數(shù)f(x),但定義f(x)的方式不同,請參考例題。

lsqcurvefit用以求含參量x(向量)的向量值函數(shù)F(x,xdata)=(F(x,xdata1),…,F(xiàn)(x,xdatan))T中的參變量x(向量),使得90輸入格式:(1)x=lsqcurvefit(‘fun’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’,x0,xdata,ydata,lb,ub);

(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,lb,ub,options);(4)[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,options,funval]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(6)[x,options,funval,Jacob]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);fun是一個事先建立的定義函數(shù)F(x,xdata)

的M-文件,自變量為x和xdata說明:x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知數(shù)據(jù)點選項見無約束優(yōu)化91

lsqnonlin用以求含參量x(向量)的向量值函數(shù)

f(x)=(f1(x),f2(x),…,fn(x))T

中的參量x,使得

最小。其中fi(x)=f(x,xdatai,ydatai)

溫馨提示

  • 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
  • 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯(lián)系上傳者。文件的所有權(quán)益歸上傳用戶所有。
  • 3. 本站RAR壓縮包中若帶圖紙,網(wǎng)頁內(nèi)容里面會有圖紙預覽,若沒有圖紙預覽就沒有圖紙。
  • 4. 未經(jīng)權(quán)益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
  • 5. 人人文庫網(wǎng)僅提供信息存儲空間,僅對用戶上傳內(nèi)容的表現(xiàn)方式做保護處理,對用戶上傳分享的文檔內(nèi)容本身不做任何修改或編輯,并不能對任何下載內(nèi)容負責。
  • 6. 下載文件中如有侵權(quán)或不適當內(nèi)容,請與我們聯(lián)系,我們立即糾正。
  • 7. 本站不保證下載資源的準確性、安全性和完整性, 同時也不承擔用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。

評論

0/150

提交評論