matlab符號(hào)計(jì)算及其應(yīng)用

MATLAB

軟件及其應(yīng)用

ApplicationofMatlabLanguage1

MATLAB符號(hào)運(yùn)算

（Symbolic）2Matlab

符號(hào)運(yùn)算介紹

Matlab符號(hào)運(yùn)算是通過(guò)符號(hào)數(shù)學(xué)工具箱（SymbolicMathToolbox）來(lái)實(shí)現(xiàn)的。Matlab符號(hào)數(shù)學(xué)工具箱是建立在功能強(qiáng)大的Maple軟件的基礎(chǔ)上的，當(dāng)Matlab進(jìn)行符號(hào)運(yùn)算時(shí)，它就請(qǐng)求Maple軟件去計(jì)算并將結(jié)果返回給Matlab。

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

符號(hào)運(yùn)算特點(diǎn)計(jì)算以推理方式進(jìn)行，因此不受計(jì)算誤差累積所帶來(lái)的困擾。

符號(hào)計(jì)算指令的調(diào)用比較簡(jiǎn)單，與數(shù)學(xué)教科書(shū)上的公式相近。

符號(hào)計(jì)算可以給出完全正確的封閉解，或任意精度的數(shù)值解（封閉解不存在時(shí)）。

符號(hào)計(jì)算所需的運(yùn)行時(shí)間相對(duì)較長(zhǎng)。4Matlab

符號(hào)運(yùn)算舉例求一元二次方程ax2+bx+c=0

的根>>

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

x)2

的一次導(dǎo)數(shù)>>

x=sym('x');>>

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

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

上的定積分>>

symsabx;>>

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

sym

函數(shù)用來(lái)建立單個(gè)符號(hào)變量，一般調(diào)用格式為：符號(hào)對(duì)象的建立：sym

和syms符號(hào)對(duì)象的建立例：>>

a=sym('a')符號(hào)變量

=

sym(A)參數(shù)A

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

是符號(hào)變量b

是符號(hào)常量>>

b=sym(1/3)C

是符號(hào)矩陣>>

C=sym('[1ab;cd]')7符號(hào)對(duì)象的建立：sym

和syms符號(hào)對(duì)象的建立

syms

命令用來(lái)建立多個(gè)符號(hào)變量，一般調(diào)用格式為：syms符號(hào)變量1符號(hào)變量2...符號(hào)變量n例：>>

symsabc>>

a=sym('a');>>

b=sym('b');>>

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

函數(shù)直接建立符號(hào)表達(dá)式。

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

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

y=sin(x)+cos(x)符號(hào)表達(dá)式的建立>>symsx;>>

y=sin(x)+cos(x)9Matlab符號(hào)運(yùn)算采用的運(yùn)算符和基本函數(shù)，在形狀、名稱(chēng)和使用上，都與數(shù)值計(jì)算中的運(yùn)算符和基本函數(shù)完全相同符號(hào)對(duì)象的基本運(yùn)算基本運(yùn)算符普通運(yùn)算：+

-

*

\

/

^數(shù)組運(yùn)算：.*

.\

./

.^矩陣轉(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符號(hào)對(duì)象的基本運(yùn)算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ù)、對(duì)數(shù)函數(shù)等基本函數(shù)11查找符號(hào)表達(dá)式中的符號(hào)變量若表達(dá)式中有兩個(gè)符號(hào)變量與x

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

按字母順序列出符號(hào)表達(dá)式expr

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

按順序列出expr

中離x

最近的N

個(gè)符號(hào)變量常量pi,

i,

j

不作為符號(hào)變量12例：>>

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

findsym(f)>>

findsym(f,3)>>

findsym(f,1)findsym

舉例13符號(hào)表達(dá)式的替換subs(f,x,a)用a

替換字符函數(shù)

f

中的字符變量x

a

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

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

應(yīng)該具有與x

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

用給定的數(shù)據(jù)替換符號(hào)表達(dá)式中的指定的符號(hào)變量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例：指出下面各條語(yǔ)句的輸出結(jié)果f=2*u15符號(hào)矩陣>>

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

使用sym

函數(shù)直接生成

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

符號(hào)矩陣中元素的引用和修改>>

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六類(lèi)常見(jiàn)符號(hào)運(yùn)算因式分解、展開(kāi)、合并、簡(jiǎn)化及通分等計(jì)算極限計(jì)算導(dǎo)數(shù)計(jì)算積分符號(hào)求和代數(shù)方程和微分方程求解17因式分解因式分解factor(f)>>

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

factor(f)

factor

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

s=factor(100)>>

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

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

expand(f)

多項(xiàng)式展開(kāi)

三角函數(shù)展開(kāi)>>

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

expand(f)19合并同類(lèi)項(xiàng)合并同類(lèi)項(xiàng)collect(f,v):按指定變量v

進(jìn)行合并collect(f):按默認(rèn)變量進(jìn)行合并>>

symsxy;

>>

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

collect(f)>>

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

對(duì)f

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

y

為f

的最簡(jiǎn)短形式，How中記錄的為簡(jiǎn)化過(guò)程中使用的方法。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ù)簡(jiǎn)化函數(shù)簡(jiǎn)化y=simplify(f):

對(duì)f

進(jìn)行簡(jiǎn)化>>

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

simplify(f)>>

symscalphabeta;

>>

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

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

symsx;

>>

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

y1=simplify(f)>>

g1=simple(f)>>

g2=simple(g1)多次使用simple

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

N

為通分后的分子，D

為通分后的分母>>

symsxy;

>>

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

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

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

多項(xiàng)式

horner

多項(xiàng)式：嵌套形式的多項(xiàng)式>>

symsx;

>>

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

g=horner(f)例：25計(jì)算極限limit(f,x,a):計(jì)算limit(f,a):當(dāng)默認(rèn)變量趨向于a時(shí)的極限limit(f):計(jì)算a=0時(shí)的極限limit(f,x,a,'right'):計(jì)算右極限limit(f,x,a,'left'):計(jì)算左極限例：計(jì)算，>>

symsxhn;

>>

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

M=limit((1-x/n)^n,n,inf)26計(jì)算導(dǎo)數(shù)g=diff(f,v)：求符號(hào)表達(dá)式

f

關(guān)于

v

的導(dǎo)數(shù)g=diff(f)：求符號(hào)表達(dá)式

f

關(guān)于默認(rèn)變量的導(dǎo)數(shù)g=diff(f,v,n)：求

f

關(guān)于

v

的

n階導(dǎo)數(shù)

diff>>

symsx;>>

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

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

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

I=int(f,x)>>

K=int(exp(-x^2),x,0,inf)例：計(jì)算和28符號(hào)求和>>

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)于默認(rèn)變量求和

例：計(jì)算級(jí)數(shù)及其前100項(xiàng)的部分和例：計(jì)算函數(shù)級(jí)數(shù)>>

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

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

可以是用字符串表示的方程、符號(hào)表達(dá)式或符號(hào)方程；

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

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

y

為輸出的解，

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

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

為自變量例1：求微分方程的通解，并驗(yàn)證。>>

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

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

幾點(diǎn)說(shuō)明如果省略初值條件，則表示求通解；如果省略自變量，則默認(rèn)自變量為t

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

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

表示對(duì)自變量的導(dǎo)數(shù)，如：Dyy'；

D2yy''；

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

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

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

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

ezplot(y);33微分方程求解例3：求微分方程組在初值條件下的特解，并畫(huà)出解函數(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í)，如果所給的輸出個(gè)數(shù)與方程個(gè)數(shù)相同，則方程組的解按詞典順序輸出。34其它運(yùn)算反函數(shù)finverse(f,v)：求f

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

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

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

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

g1=finverse(f,x)>>

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

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

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

1、模型：命令：x=linprog（c，A，b）

2、模型：minz=cX

命令：x=linprog（c，A，b，Aeq,beq）注意：若沒(méi)有不等式：存在，則令A(yù)=[]，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]若沒(méi)有等式約束:,則令A(yù)eq=[],beq=[].[2]其中X0表示初始點(diǎn)4、命令：[x,fval]=linprog(…)返回最優(yōu)解ｘ及ｘ處的目標(biāo)函數(shù)值fval.38解編寫(xiě)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解:編寫(xiě)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應(yīng)用—

解非線性規(guī)劃

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

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

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

（1）其中，是定義在Rn

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

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

0

m;1,2,...,

0..

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

其中X為n維變?cè)蛄浚珿(X)與Ceq(X)均為非線性函數(shù)組成的向量。48

1．首先建立M文件fun.m,用來(lái)定義目標(biāo)函數(shù)F（X）:functionf=fun(X);f=F(X);MATLAB求解上述問(wèn)題，基本步驟分三步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(…)輸出極值點(diǎn)M文件迭代的初值參數(shù)說(shuō)明變量上下限fmincon函數(shù)可能會(huì)給出局部最優(yōu)解，這與初值X0的選取有關(guān)．501．寫(xiě)成標(biāo)準(zhǔn)形式：

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．運(yùn)算結(jié)果為：

x=0．76471．0588fval=-2．0294521．先建立M文件fun4．m定義目標(biāo)函數(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．運(yùn)算結(jié)果為：

x=-1．22501．2250fval=1．895154例1．先建立M文件fun．m定義目標(biāo)函數(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應(yīng)用—

微分方程解法

57要求目的主要內(nèi)容2、學(xué)會(huì)用Matlab求微分方程的數(shù)值解.1、學(xué)會(huì)用Matlab求簡(jiǎn)單微分方程的解析解.1、求簡(jiǎn)單微分方程的解析解.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化簡(jiǎn)

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)和科研中所處理的微分方程往往很復(fù)雜且大多得不出一般解。而在實(shí)際上對(duì)初值問(wèn)題，一般是要求得到解在若干個(gè)點(diǎn)上滿足規(guī)定精確度的近似值，或者得到一個(gè)滿足精確度要求的便于計(jì)算的表達(dá)式。因此，研究常微分方程的數(shù)值解法是十分必要的。63（二）建立數(shù)值解法的一些途徑1、用差商代替導(dǎo)數(shù)

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

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

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

2、使用Matlab軟件求數(shù)值解時(shí)，高階微分方程必須等價(jià)地變換成一階微分方程組.注意: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的圖形為實(shí)線，y2的圖形為“*”線，y3的圖形為“+”線.70Matlab的應(yīng)用

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

：最鄰近插值‘linear’

：線性插值；‘spline’

：三次樣條插值；‘cubic’

：立方插值。缺省時(shí)：分段線性插值。注意：所有的插值方法都要求x是單調(diào)的，并且xi不能夠超過(guò)x的范圍。75例：在1-12的11小時(shí)內(nèi)，每隔1小時(shí)測(cè)量一次溫度，測(cè)得的溫度依次為：5，8，9，15，25，29，31，30，22，25，27，24。試估計(jì)每隔1/10小時(shí)的溫度值。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.擬合問(wèn)題引例77擬合問(wèn)題引例1溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032已知熱敏電阻數(shù)據(jù)：求600C時(shí)的電阻R。

設(shè)

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

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

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

使f(x)

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

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

f之間的關(guān)系？問(wèn)題：給定一批數(shù)據(jù)點(diǎn)，需確定滿足特定要求的曲線或曲面解決方案：若不要求曲線（面）通過(guò)所有數(shù)據(jù)點(diǎn)，而是要求它反映對(duì)象整體的變化趨勢(shì)，就是數(shù)據(jù)擬合，又稱(chēng)曲線擬合或曲面擬合。若要求所求曲線（面）通過(guò)所給所有數(shù)據(jù)點(diǎn)，就是插值問(wèn)題；81曲線擬合問(wèn)題最常用的解法——線性最小二乘法的基本思路第一步:先選定一組函數(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

的準(zhǔn)則（最小二乘準(zhǔn)則）：使n個(gè)點(diǎn)（xi,yi)與曲線y=f(x)的距離i的平方和最小

。記

問(wèn)題歸結(jié)為，求

a1,a2,…am

使

J(a1,a2,…am)

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

83線性最小二乘法的求解定理：當(dāng)RTR可逆時(shí)，超定方程組（3）存在最小二乘解，

且即為方程組

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

（數(shù)組）輸入同長(zhǎng)度數(shù)組X，Y擬合多項(xiàng)式

次數(shù)86即要求出二次多項(xiàng)式:中的使得:例對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合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）計(jì)算結(jié)果：Ａ=[-9.8108,20.1293,-0.0317]882）計(jì)算結(jié)果：Ａ=[-9.8108,20.1293,-0.0317]解法2．用多項(xiàng)式擬合的命令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ù)點(diǎn)和擬合曲線的圖形891.lsqcurvefit已知數(shù)據(jù)點(diǎn)：xdata=（xdata1，xdata2，…，xdatan）

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

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是一個(gè)事先建立的定義函數(shù)F(x,xdata)

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

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

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

中的參量x，使得

最小。其中fi（x）=f（x，xdatai，ydatai）