MATLAB程序設(shè)計(jì)與應(yīng)用下

1、matlab程序設(shè)計(jì)與實(shí)例應(yīng)用（下）第六章matlab的數(shù)值計(jì)算錯(cuò)誤!未定義書簽。6-1多項(xiàng)式的運(yùn)算846-1-1多項(xiàng)式的生成和表達(dá)846-1-2多項(xiàng)式的乘除856-1-3多項(xiàng)式的求導(dǎo)866-1-4多項(xiàng)式的求根876-2數(shù)據(jù)分析896-1-1極值、均值、標(biāo)準(zhǔn)差和中位值的計(jì)算896-2-2曲線擬合906-3數(shù)值積分和數(shù)值微分926-3-1微分和積分的數(shù)學(xué)表達(dá)式926-3-2函數(shù)數(shù)值積分926-3-3數(shù)值微分976-3 -般非線性方程組的數(shù)值解996- 4微分方程的數(shù)值解1006-4-1微分方程的基本形式1006-4-2 一階常微分方程的求解1006- 4-3二階常微分方程的求解103第七章 ma

2、tlab的符號(hào)計(jì)算1497- 1符號(hào)表達(dá)式的生成1517- 1-1 符號(hào)表達(dá)式白勺建立1517-1-2符號(hào)矩陣的修改1527-2符號(hào)表達(dá)式的基本運(yùn)算1537-2-1基本代數(shù)運(yùn)算1537-2-2因式分解、展開、化簡1547-2-3符號(hào)與數(shù)值間的轉(zhuǎn)化1557-3符號(hào)矩陣1567-3-1符號(hào)矩陣的創(chuàng)立1577-3-2符號(hào)矩陣的基本運(yùn)算1587-3-3符號(hào)矩陣的簡化1617-3-4符號(hào)矩陣的特征值和特征向量1637-4符號(hào)求極限1657-4-1 單變量函數(shù)的極限1657-4-2多變量函數(shù)的極限1667-5 微分1677-6 積分1687-7符號(hào)求解方程1697-1-1線性方程1697-1-2非線性方程

3、1717-1-3常微分方程1727-1-4符號(hào)函數(shù)的圖形顯示1747-1-5符號(hào)函數(shù)的顯示177第六章matlab數(shù)值計(jì)算6-1多項(xiàng)式的運(yùn)算6-1-1多項(xiàng)式的生成和表達(dá)1. 多項(xiàng)式的表達(dá)在matlab環(huán)境下多項(xiàng)式是用向量的形式表達(dá)的.向量最右邊的元素表示多項(xiàng)式的0階，向左數(shù)依 次表示多項(xiàng)式的第1階、第2階、第3階。例如多項(xiàng)式5兀4+3 + 2x+1表示為:5032 102. 多項(xiàng)式的生成語法：p=ploy(ma)說明：1. 若ma為方陣，則生成的多項(xiàng)式為p為方陣ma的特征多項(xiàng)式。2. 若ma為向量，則向量和多項(xiàng)式滿足這樣一種關(guān)系：ma = r, r2 rn ,生成的多項(xiàng)式為: (x-r1)(x

4、-r2)(x-r3)-«(x-rn) = a0xn+a1xn-,+a2xn-2 +3. 直接輸入的方式生成多項(xiàng)式。例6-1利用方陣m二5 6 7;8 9 1;11 12 13生成一個(gè)多項(xiàng)式(為方陣m的特征多項(xiàng)式)。程序設(shè)計(jì):»clearm=5 6 7;8 9 1;11 12 13;p=poly(m);%產(chǎn)生多項(xiàng)式的向量表達(dá)式px=poly2str(p, 'x');%生成常見的多項(xiàng)式表示形式p,px運(yùn)行結(jié)果：p =1.0000 -27.000090.000054.0000px =xa3 27 xa2 + 90 x + 54例6-2利用向量a=2 3 4 5生成

5、一個(gè)多項(xiàng)式。程序設(shè)計(jì)：»cleara二2 3 4 5;%產(chǎn)生多項(xiàng)式的向量表達(dá)式%生成常見的多項(xiàng)式表示形式p=poly(a);px=poly2str(p, 乂)；p,px運(yùn)行結(jié)果：p =1-1471154120px =xa4 14 xa3 + 71 xa2 154 x+ 1206-1-2多項(xiàng)式的乘除語法：a. c=conv(a5b)b. q,r=decony(c,a)說明：1. a、b和c分別是多項(xiàng)式的向量表示形式。a表示兩個(gè)多項(xiàng)式的乘積運(yùn)算，b表示兩個(gè)多項(xiàng)式的 除法運(yùn)算。2. q表示除運(yùn)算的商，r表示除運(yùn)算的余數(shù)。例6-3求多項(xiàng)式f(x) = x2 +5x和g(x) = 2x+1的乘

6、積m(x)o 程序設(shè)計(jì)：%第一個(gè)多項(xiàng)式f(x)%第二個(gè)多項(xiàng)式g(x)%求兩個(gè)多項(xiàng)式的乘積%用常用的方式表示多項(xiàng)式的積»cleara=l 5 0; b=2 1; c=conv(a,b);mx=poly2str(c, 'x'); c,mx%end運(yùn)行結(jié)果：2mx =112 xa3 + 11 xa2 + 5 x例6-4求多項(xiàng)式f(x) = x2+5x和g(x) = 2x + 1的除運(yùn)算d(x)。 程序設(shè)計(jì)：»clearc=l 5 0;a=2 1;q,r=deconv(c,a);%第一個(gè)多項(xiàng)式f(x)%第二個(gè)多項(xiàng)式g(x)%求 f(x)/ g(x)dx=poly2s

7、tr(q, 'x');q,r,dx%end運(yùn)行結(jié)果：q =0.50002.2500%用常用的方式表示多項(xiàng)式的積00-2.2500dx =0.5 x + 2.25程序說明：1. 在運(yùn)行結(jié)果中變量q是f(x)除以g(x)的商，而r則是除不盡的余數(shù)。2. 運(yùn)行結(jié)果變量dx表示的商沒有加上余數(shù)。6-1-3多項(xiàng)式的求導(dǎo)語法：dp=polyder(p)說明：p為向量表示的多項(xiàng)式。例6-5求多項(xiàng)式f(x) = x2+5x和g(x) = 2x + 1的一階導(dǎo)數(shù)。我們?nèi)菀字酪陨蟽蓚€(gè)方程的導(dǎo)數(shù)手工驗(yàn)算結(jié)果為：(x) =2x+5和g * (x)=2我們看matlab的計(jì)算結(jié)果。 程序設(shè)計(jì)：

8、7;clearf=l 5 0;g 二2 1；df=polyder(f);dg= polyder(g);dfx=poly2str(df, 'x');dgx=poly2str(dg, 'x');df,dg,dfx,dgx%end運(yùn)行結(jié)果：df =25%第一個(gè)多項(xiàng)式f(x) %第二個(gè)多項(xiàng)式g(x) %求f(x)的導(dǎo)數(shù) %求6&)的導(dǎo)數(shù)dg =dfx 二2 x + 5 dgx =96-1-4多項(xiàng)式的求根語法：a.r=roots(p)說明：另外還有一種通過先求多項(xiàng)式的伴隨矩陣，然后再求特征值的方法也可以求得多項(xiàng)式的根。例66求多項(xiàng)式f(x) = x2+5x和g(x

9、) = 2x +1的根。我們?nèi)菀字婪匠痰母鶠椋篺(x)為：xl=o； x2=-5g(x)為：x2=-l/2我們看matlab的計(jì)算結(jié)果。程序設(shè)計(jì)：»clear仁15 0;%第一個(gè)多項(xiàng)式f(x)g 二2 1;%第二個(gè)多項(xiàng)式g(x)rf=roots(f);%求f(x)的根rg= roots (g)；%求6匕)的根rf, rg,%end運(yùn)行結(jié)果：rf =0-5rg =-0.5000程序說明:從運(yùn)行結(jié)果我們可以看到f(x)的根為0和-5, g(x)的根為-0.50,這和我們手工驗(yàn)算的結(jié)果完全一致。例6-7求多項(xiàng)式f(x) = x2+5x-3 = 0和g(x) = 2x3 + x2+1 =

10、0 的根。程序設(shè)計(jì)：%第一個(gè)多項(xiàng)式f(x)»clearhl 5-3;g=l2 1 0 1j； rf=roots(f); rg= roots (g);rf, rg,%end 運(yùn)行結(jié)果：rf =-5.54140.5414%第二個(gè)多項(xiàng)式g(x)%求f(x)的根%求6仇)的根-1.00000.2500 +0.6614i0.2500 0.6614i程序說明：1. 我們可以看到例6-6和例6-7求得的根和我們手工計(jì)算的結(jié)果是一致的，其中例6-7多項(xiàng)式的 根出現(xiàn)了虛根。2. 我們用求伴隨矩陣的方法對(duì)例6-7再做一次求根運(yùn)算：»clearf=15 3; g=2 10 1;cf=compan

11、(f); cg=compan(g);cf, eg,%end運(yùn)行結(jié)果：cf =%第一個(gè)多項(xiàng)式f(x) %第二個(gè)多項(xiàng)式g(x) %求f(x)的伴隨矩陣 %求6匕)的伴隨矩陣-5-3cg =-0.5000-0.50001.00001.0000ccf=eig(cf) ccf=eig(cg) 運(yùn)行結(jié)果：ccf =%求特征根%求特征根-4.3028-0.6972ccf =0.2500 + 0.661410.2500-0.6614i我們看到兩種方法的結(jié)果是一致的。6-2數(shù)據(jù)分析6-1-1極值、均值、標(biāo)準(zhǔn)差和中位值的計(jì)算語法：a. pmax=max(x)b. pmin=min(x)c. pmean=mean(

12、x)d. pstd=std(x)e. pmedian=median (x)說明：1. max(x)、min(x)、mean(x)、std(x)、median (x)分別是用來求數(shù)組或矩陣的最大值、最小值、 均值、標(biāo)準(zhǔn)差、中位值。2. 這里x可以是數(shù)組也可以是矩陣。如果是數(shù)組則對(duì)整個(gè)數(shù)組進(jìn)行計(jì)算，如果是矩陣則是分別 對(duì)矩陣的每個(gè)列向量進(jìn)行運(yùn)算。例6-8對(duì)一組隨機(jī)數(shù)組進(jìn)行均值、方差和屮位值的計(jì)算。程序設(shè)計(jì)：»clearx=randn(l jo);%生成一個(gè)10元素的數(shù)組pmax=max(x);%計(jì)算數(shù)組x的最大值pmin=min(x);%計(jì)算數(shù)組x的最小值pmean=mean(x);ps

13、td=std(x);%計(jì)算數(shù)組x的平均值%計(jì)算數(shù)組x的標(biāo)準(zhǔn)差psqu= pstda2;%計(jì)算數(shù)組x的方差pmed=median (x);%計(jì)算數(shù)組x的中位值answers= pmax pmin pmean pstd psqu pmed;x,answers運(yùn)行結(jié)果：columns 1 through 6-0.4326-1.66560.12530.2877 1.14651.19091.1892-0.03760.32730.1746answers =1.1909-1.66560.00130.90340.81620.1500程序說明：1. 注意方差為標(biāo)準(zhǔn)差的平方。2. 需要注意的是matlab的變量定

14、義是區(qū)分大小寫的。例6-9均值、標(biāo)準(zhǔn)差、中位值的計(jì)算。%生成一個(gè)6x 3的隨機(jī)矩陣%計(jì)算矩陣x的最大值%計(jì)算矩陣x的最小值%計(jì)算矩陣x的平均值%計(jì)算矩陣x的標(biāo)準(zhǔn)差%計(jì)算矩陣x的中位值對(duì)一個(gè)隨機(jī)矩陣進(jìn)行最大值、最小值、 程序設(shè)計(jì)：»clearx=randn(6,3)；pmax=max(x);pmin=min(x);pmean=mean(x);pstd=std(x);pmed=median (x);x,pmax,pmin,pmean,pstd,pmed;運(yùn)行結(jié)果:-0.18671.06680.71430.72580.05931.6236-0.5883-0.0956-0.69182.183

15、2-0.83230.8580-0.13640.29441.25400.1139-1.3362-1.5937pmax =2.18321.06681.6236pmin =-0.5883-1.3362-1.5937pmean =0.3519-0.14060.3607pstd =0.99620.84821.2404a.b.622語法:p 二 polyfit(xyn) yval 二 polyval(rx)說明：1. polyfit(x,y,n)根據(jù)輸入數(shù)據(jù)x和y生成一個(gè)n階的擬合多項(xiàng)式。2. polyval(p,x)根據(jù)數(shù)據(jù)x,用擬合多項(xiàng)式p生成擬合好的數(shù)據(jù)。3. 這里x和y構(gòu)成了一系列數(shù)據(jù)點(diǎn)的坐標(biāo)。例

16、6t0下列數(shù)據(jù)為一段時(shí)間的日氣溫平均實(shí)際數(shù)據(jù)，求該數(shù)據(jù)的擬合方程。設(shè)這組氣溫?cái)?shù)據(jù)為:18 19 17 18 20 22 22 23 21 21 20 22 23 23 22 程序設(shè)計(jì)：%時(shí)間%和吋間相對(duì)應(yīng)的日平均氣溫%生成擬合多項(xiàng)式(向量形式)，階次為3%生成擬合多項(xiàng)式(字符形式)%生成擬合值»cleard=l:15;t=18 19 17 18 20 22 22 23 21 21 20 22 23 23 22; p=polyfit(d,t,3)；px=poly2str(p, 'x');pv=polyval(p,d);p，pxplot(d,t, ':*f,d,p

17、v,'-k', 'linewidlh； 2)title('fontsize 18bf數(shù)據(jù)擬合圖像:'color1,y1) xlabel cfontsize 14rm時(shí)間d軸；'color1,'k1)ylabel cfontsize 14rm氣溫t軸；'color','k') gtext('fontsize 16fonlname宋體實(shí)際曲線；'color','r') gtext('fontsize 16fonlname宋體擬合曲線；'color'

18、,'k') %end運(yùn)行結(jié)果：p =0.0010-0.05380.964416.4623px =0.0010474 xa3 0.053823 xa2 + 0.96436 x + 16.4623圖形如6-1所示：23數(shù)據(jù)擬合圖像22擬合曲線18- 夕%仃11051015時(shí)間d軸圖61擬合曲線圖程序說明:1. <1表示時(shí)間第一天、第二天，t表示每天的氣溫。2. 本例我們選取的擬合多項(xiàng)式的階次是3。3. p為向量形式的擬合多項(xiàng)式，px為字符形式的擬合多項(xiàng)式。4. plot為繪圖函數(shù)，它的兩個(gè)參數(shù)分別為圖像的橫坐標(biāo)和縱坐標(biāo)。4個(gè)參數(shù)則分別為兩條曲線的 橫坐標(biāo)和縱坐標(biāo)。5. 圖中的

19、箭頭可在繪圖窗口中打開insert菜單，選擇arrow然后用鼠標(biāo)從起始位置拖動(dòng)到終止位 置即可。例 611根據(jù)上例中的數(shù)據(jù)建立的天氣溫度，擬合方程求5天后的天氣氣溫是多少？程序設(shè)計(jì)及運(yùn)行結(jié)果：»clearx=20;t=0.0010474*xa3 - 0.053823*xa2 + 0.96436檢 + 16.4623t =22.5995程序設(shè)計(jì)及運(yùn)行結(jié)果：1. 依照上例710天的時(shí)間序列d,5天后應(yīng)該是20天，所以令x=20o2. t即為上例中我們求得的擬合公式，結(jié)果求得5天后的日平均氣溫為22.5995。6-3數(shù)值積分和數(shù)值微分6-3-1微分和積分的數(shù)學(xué)表達(dá)式對(duì)函數(shù)f(x)在區(qū)間xi

20、,x2上求積分在數(shù)學(xué)上常常表示為x2f= jf(x)dxx|對(duì)函數(shù)f(x)求導(dǎo)數(shù)運(yùn)算在數(shù)學(xué)上常常表示為：df(x)dx命令函數(shù)：1. 連續(xù)被積函數(shù) quad(t,a,b,t,trace) quadl(f,a,b,t,trace) quad8(rf5a,b,t,trace)2. 離散被積函數(shù) trapz(x,y) cumsum(y)6-3-2函數(shù)數(shù)值積分自適應(yīng)遞推辛普生(simpso)法求積分(精度較差，建議不用) 自適應(yīng)遞推牛頓一科西(newton otes)法求積分 自適應(yīng)遞推牛頓一科西(newton一cotes)法求積分(舊語法)梯形法求數(shù)值積分歐拉法求數(shù)值積分說明:1.1) 參數(shù)f是被積

21、函數(shù)表達(dá)式字符串或者函數(shù)文件名。2) 參數(shù)a、b定義函數(shù)積分的上下限。3) 參數(shù)t定義積分的精度(默認(rèn)值1.0e-6); trace設(shè)置是否用圖形展示積分過程，1為展示、0 不展示。2.1) trapz積分中給出y相對(duì)于x的積分值。當(dāng)y是(m*n)矩陣時(shí)，積分對(duì)y的列向量分別進(jìn) 行，得到一個(gè)(l*n)矩陣是y的列向量對(duì)應(yīng)于x的積分結(jié)杲。2) cumsum對(duì)y的列向量進(jìn)行積分運(yùn)算，采用等距離單位步長，但精度比較差。積分結(jié)果和 y是同維的，還要注意計(jì)算結(jié)果除以采樣頻率(或乘以步長)才是實(shí)際的枳分序列。例 6-1271設(shè)y(x) = sinx ,用辛普生法求s = jsinxdx oo這個(gè)積分可以用

22、手工算出來是2,我們看matlab的計(jì)算結(jié)果。程序設(shè)計(jì)及運(yùn)行結(jié)果：»clears=quad(,sin(x),o,pi)s =2.0000例 6-13兀設(shè)y(x) = sinx ,用牛頓一科西法求s = jsin xdx。0這個(gè)積分可以用手工算出來是2,我們看matlab的計(jì)算結(jié)果。程序設(shè)計(jì)及運(yùn)行結(jié)果：»clears=quadl(,sin(x),o,pi)s =2.0000例 6-14設(shè)y(x) = sinx,用牛頓一科西法求y(x)在區(qū)間0,2兀內(nèi)每間隔為兀/100小段上的積分值序列， 并畫出這個(gè)序列的圖像。程序設(shè)計(jì)及運(yùn)行結(jié)果：»clearx=ri；t 二 0:p

23、i/100:2*pi;m=max(size(t)-l;for n=l:m;x(n)=quadl(,sin(x)t(n),t(n+1);endplot(x,'k,linewidth； 2)titlecfontsize 18bf積分值序列的圖像;'colof,k)grid onx%end運(yùn)行結(jié)果(如圖6-2所示)：x =columns 1 through 60.00050.00150.00250.00340.00440.0054columns 195 through 200-0.0054-0.0044-0.0034-0.0025-0.00150.0005積分值序列的圖像圖6-2積分

24、值序列圖像程序說明：1. 在運(yùn)行結(jié)果里我們省掉了序列x第7到194個(gè)元素之間的元素。2. 積分序列x和序列的圖像更加直接地說明了積分的物理意義，一個(gè)累積的過程，對(duì)序列x求和, 我們看到其結(jié)果和積分結(jié)果一樣。從圖上可以看出最大值和最小值分別為士 sin(p"2)*pi/10()。對(duì)序列x求和：程序設(shè)計(jì)：»clearx 二;t=0:pi/100:2*pi;m=max(size(t)-l;sn=0;for n=l :m;x(n)=quadl('sin(x)',t(n),t(n+1);sn=sn+x(n);endsn%encl運(yùn)行結(jié)果：sn =-1.6827e-01

25、6對(duì)sinx直接積分 程序設(shè)計(jì)：»clear s=quadl(,sin(x),0,2*pi) %end運(yùn)行結(jié)果：1.7588e-016程序說明：兩者結(jié)果一致，均接近0,實(shí)際可看作0。如果將積分區(qū)間改為0,7t,結(jié)果與我們熟知的結(jié)果一 致等于2。例6t52設(shè)y(x) = (sin2 x + x3)/(x -2cos3x),求2 jy(x)dx «o程序設(shè)計(jì)及運(yùn)行結(jié)果：»cleary=quadl('(sin(x).a2+x.a3)./(x-2.*cos(3.*x)',0,2)%endans =11.1943例 616用離散數(shù)據(jù)表示sinx程序設(shè)計(jì)及運(yùn)行

26、結(jié)果(如圖6-3):»clearx=0:0.001:pi;y=sin(x);z=trapz(x,y);s=cumsum(y)*0.001;plot(x,s,'k；'linewidth', 2)積分區(qū)間為0, tl分別用離散積分方法對(duì)其進(jìn)行積分%離散化積分變量%生成被積函數(shù)的離散序列%梯形法求積分%歐拉法求積分%繪出積分結(jié)果序列s的圖像title('fontsize 18bf歐拉積分;'color'/k')xlabel('fontsize 18)bfx', 'color','k')y

27、label('fontsize 18 bfcumsum*0.001'color','k') grid onzz =2.00002歐拉積分18x圖6-3 cumsum積分結(jié)果序列圖像0.40.200.511 522 5335614 28 6 o.q50.0*e3we3o程序說明：本例中梯形法的結(jié)果是2,歐拉法積分的結(jié)果是一個(gè)和被積函數(shù)y同維的序列s。所以為了表示積 分結(jié)果我們給出了 s序列的圖像。注意歐拉積分的序列值乘以積分步長(或除以間隔頻率)為實(shí)際積 分值。例 6-172設(shè) y(x) = (sin2 x + x')/(x-2cos3x),求s

28、= jy(x)dx。0程序設(shè)計(jì)及運(yùn)行結(jié)果(如圖6-4)：»clearx=0:0.001:2;%離散化積分變量y=(sin(x).a2+x.a3)./(x-2.*cos(3.*x);z=trapz(x,y);s=cumsum(y)*0.001zplot(x,s；-kvlinewidth 2)title(fontsize 18bf歐拉積分'color','k')xlabel(、fontsize18bfx', ,color,k,)ylabel('fontsize 18bfcumsum*0.001 'color','k&

29、#39;)grid on%endz =10.9607s =columns l through 60 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000(中間數(shù)據(jù)省略)columns 1999 through 200110.795610.905311.0161程序說明：1. 輸入公式時(shí)一定要注意算符都是對(duì)數(shù)組進(jìn)行運(yùn)算的，所以是"" ./"。2. 我們可以看到離散的積分計(jì)算結(jié)果和連續(xù)函數(shù)的積分計(jì)算結(jié)果很相近，但不完全相同。歐拉積分0"0020.40 60 81121 41.6182x圖6-4 cumsum積分結(jié)果序列圖像eo.o+

30、e3we3012106-3-3數(shù)值微分命令函數(shù)：d=diff(y)說明：y為一組離散序列，d為與y同維的離散矩陣，則微分運(yùn)算是對(duì)矩陣的每個(gè)列向量微分。例 6-18設(shè)函數(shù)y(x) = sinx和z = 2x2 + 3x2 + x-1進(jìn)行微分運(yùn)算，并畫出計(jì)算結(jié)果的圖像。程序設(shè)計(jì)：»clear%生成一個(gè)從-pi到pi間隔為0.001的序列%對(duì)丫進(jìn)行微分運(yùn)算%對(duì)進(jìn)行求導(dǎo)運(yùn)算%對(duì)乙進(jìn)行微分運(yùn)算%對(duì)乙進(jìn)行求導(dǎo)運(yùn)算%產(chǎn)生兩行一列的繪圖區(qū)間的第一個(gè)繪圖區(qū)間dx=0.001;x=-3:dx:3;xl=-3:dx:3-dx;y=sin(x);z=2.*x.a3+3.*x.a2+x-l;y=diff(y)

31、;dy= diff(y)/dx;z=diff(z);dz=diff(z)/dx;subplot(2,2,l);plot(y；-k', linewidth； 1)title('fontsize 12fontnametime new roman deltay(x) ； 'color；'k')grid onsubplot(2,2,2);%產(chǎn)生兩行一列的繪圖區(qū)間的第二個(gè)繪圖區(qū)間plot(z ,linewidth 1)title('fontsize12fontnametime new roman) delta z (x)', 'color&

32、#39;/k')grid onsubplot(2,2,3)；%產(chǎn)生兩行一列的繪圖區(qū)間的第三個(gè)繪圖區(qū)間plot(xl,dy；-k linewidth； 1)title('fontsize 12)fontname time new roman yh(x)=cosx', 'color'/k')xlabel('fontsize 12 fontname time new roman x', 'color','k')ylabel('fontsize 12fontnametime new roman dy

33、/dx 'color','k')axis (-3 3 1 1)grid onsubplot(2,2,4)；%產(chǎn)生兩行一列的繪圖區(qū)間的第四個(gè)繪圖區(qū)間plot(xl,dz ；-k linewidth； 1)title('fontsize 12fontname time new roman z',(x)=6xa2+6x+1'color','k')xlabel(rfontsize12fontnametime new roman x', 'color','k*)ylabel('font

34、size 12fontnametime new roman dz/dx', 'color','k')axis (-3 3 -20 80)grid on%end運(yùn)行結(jié)果(如圖6-5)：az(x)xx圖65函數(shù)的微分與導(dǎo)數(shù)圖程序說明:1. 輸入公式時(shí)一定要注意算符都是對(duì)數(shù)組進(jìn)行運(yùn)算的，所以是"2. y=diff(y)和z二diff(z)和依此畫出的圖像，縱坐標(biāo)是兩函數(shù)的微分，橫坐標(biāo)是數(shù)據(jù)序列的值并不 是實(shí)際圖像的橫坐標(biāo)。由于縱坐標(biāo)是 fn =f (xn+1)-f (xn),所以縱坐標(biāo)在數(shù)值上是導(dǎo)數(shù) 的(xn+1xj 倍。3. 導(dǎo)數(shù)圖的縱坐標(biāo)是 fn

35、/(xn+i -xn) = f (xn+1)-f (xn)/(xn+1 -xn),橫坐標(biāo)的取值范國應(yīng) 該是xpxn -1/n o這樣畫出的導(dǎo)數(shù)圖像才是正確的。這是由于x取值是x到xn+1,而微分 相應(yīng)取值是到,所以兩序列的數(shù)目相差1,所以畫圖時(shí)不能用plot(x,diff(y) /dx)o6-3 一般非線性方程組的數(shù)值解命令函數(shù)：x=fsolvecfunxo,options)說明：1. 'fun'是我們要求解的方程組，可以直接輸入，不過由于方程組比較復(fù)雜，一般都先生成一個(gè)m 文件，然后再進(jìn)行調(diào)用。2. xo是給出的這個(gè)方程組的初值解，我們可以隨意地給出。3. options是命

36、令函數(shù)fsolve的參數(shù)設(shè)置項(xiàng)。因?yàn)閒solve求解的過程本身就是一個(gè)優(yōu)化的過程， 而options則是設(shè)置優(yōu)化過程的參數(shù)的。這里對(duì)于一般的非線性方程組求解我們常設(shè)置為 optiomset('display', off),意思為不顯示優(yōu)化過程。例：6-19求方程組x+y+z=0v x-y + 3z = 03x + y-z = 4的數(shù)值解。這個(gè)方程組的解經(jīng)過手工演算我們知道解為：x=4/3, y二2/3, z二2/3。程序設(shè)計(jì)與運(yùn)行結(jié)果：1. 建立函數(shù)文件xyz,并保存：function fun=xyz(x)x=x(l);y=x ;z 二 x(3);fun=zeros(3,l);

37、fun( 1 )=x+y+z;fun(2)=x-y+3*z;fun( 3 )=3 *x+y-z-4;2. 在matlab的工作區(qū)里求解：»xo=1 1 1;x=fsolve(xyz,xo, optimset ('display', 'off)1.3333-0.6667-0.6667程序說明:1. 在建立一個(gè)函數(shù)文件時(shí)一定要以function作為文件頭說明。2. 把建立的函數(shù)文件保存到matlab默認(rèn)的文件夾里，如work文件夾。3. zeros(3,l)用來產(chǎn)生一個(gè)3行1列的矩陣，矩陣的所有元素都為0。4. 我們可以看到計(jì)算的結(jié)果和我們手工計(jì)算的結(jié)果非常近似，

38、差別是計(jì)算精度引起的，fsolve的 默認(rèn)求解精度是0.001 o5. 事實(shí)上對(duì)方程組的求解過程是一個(gè)優(yōu)化過程，optimset ('display；'off)設(shè)置不顯示求解過程。例：6-20求非線性方程組x + y+ z = 0 cos x + y3 + z = 54x + 2、+ in z = 1的數(shù)值解。這個(gè)方程組手工計(jì)算是難以完成的。程序設(shè)計(jì)與運(yùn)行結(jié)果：1. 建立函數(shù)文件fun,并保存：function q=fun(t) x=t(l);y=t(2);z=t(3);q=zeros(3,l);q(l)=x+yz;q(2)=cos(x)+ya3+z-5;q(3)=4*x+2a

39、y+log(z)-1;2. 在matlab的工作區(qū)里求解：>>clearx=fsolve(fun, 1 1 1j, optimset ('display', 'off)-0.43921.45481.0156程序說明：fsolve的默認(rèn)求解精度是0.001,方程成立的精度為0.001。6-4微分方程的數(shù)值解6-4-1微分方程的基本形式f(x)的導(dǎo)數(shù)等于g(x),即：6-4-2 一階常微分方程的求解語法：a. t,y=ode45(odefun,tspan,yo)b. t,y=ode23(odefun,tspan,yo)說明：1. ode45()和ode23()是

40、兩個(gè)求解微分方程的命令函數(shù)，ode45()應(yīng)用最廣，在容許誤差大的情況下 ode23()的效率要比ode45()高一些。2. t和y分別表示微分方程的解的兩個(gè)變量的數(shù)值序列。3. odefun是待解微分方程表達(dá)式的函數(shù)文件名。4. tspan表示運(yùn)算的起至?xí)r間，是行向量，例如0 10表示運(yùn)算時(shí)間從0開始到10結(jié)束。1. yo初始狀態(tài)，用列向量表示，其元素分別表示微分方程組各個(gè)表達(dá)式的初始狀態(tài)值。例：6-21求解微分方程dy=cos xdx并畫出解在區(qū)間0 2*pi上的圖像。程序設(shè)計(jì)與運(yùn)行結(jié)果：1. 在程序編輯器中建立待解微分方程表達(dá)式的函數(shù)文件funx.m,并存入work工作目錄： funct

41、ion y=funx(x,y)y=cos(x);2. 在matlab的命令窗口求解，結(jié)果如圖66所示：>>clearx,y=ode45(funx,0 2*pi,0); %可以用 fx,yl=ode45(funx,0:0.01:2*pi,0); plot(x,y；-k linewidlh； 2)title('fontsize 18)fontnametime new roman)微分方程y"(x)=cosx求解圖像,'color','r') xlabel('fontsize 16fontnametime new roman x&

42、#39;, 'color','k')ylabel('fontsize 16fontnametime new roman y=sinx 'color','k')axis(fo 2*pi-l 1)grid on微分方程y'(x)=cosx求解圖像圖65 階微分方程解的圖像1程序說明:2. 注意函數(shù)文件建立的格式，一定要以function作文件頭，函數(shù)名為funx,兩個(gè)變量為x和y。 這里還要注意兩個(gè)變量x和y的順序，否則繪出的圖像橫坐標(biāo)和縱坐標(biāo)就顛倒了。3. 要把建立的函數(shù)文件保存在matlab默認(rèn)的文件夾里，否則在調(diào)

43、用時(shí)要指明路徑。4. ode45(funxj0 2*pi,0),注意函數(shù)文件的調(diào)用格式。本例中計(jì)算時(shí)間(x變量)是從0開始, 到2pi結(jié)束，初始值為0。5. 本例的微分方程是我們熟悉的，它的解和繪出的圖像與手工計(jì)算結(jié)果是吻合的。例：6-22求解微分方程尋一 5y = (x-l)sinx + (x + l) cos x并畫出解在區(qū)間0 5上的圖像。程序設(shè)計(jì)與運(yùn)行結(jié)果：1.在程序編輯器中建立待解微分方程表達(dá)式的函數(shù)文件funxx.m,并存入work工作目錄:function y=funxx(x,y)y=5*y+(x- l)*sin(x)+(x+1 )*cos(x);2.在matlab的命令窗口求解

44、，結(jié)果如圖6-6所示:»clear x,y=ode45(funxx,0 1,1);plot(x,y；-k 'lincwidth； 2)title('fontsize 18fontnametime new roman微分方程yn(x)=5y+(x-1 )sinx+(x+1 )cosx求解圖像; color*,y)xlabel('fontsize 16 fontname time new roman x', 'color1,'k')ylabel('fontsize 16) fontname time new roman y

45、'color1,'k')axis (0 1 0 1.8e+2)grid on微分方程 y'(x)=5y+(x-1 )sinx+(x+1 )cosx 求解圖像圖66微分方程解的圖像2程序說明:1. 在本例中我們起的函數(shù)名為funxx,兩個(gè)變量為x和y。2. 本例就比較復(fù)雜了，手工求解比較困難，但是我們看到用matlab很容易就給出了解的圖像。6-4-3二階常微分方程的求解二階常微分方程的求解使用的仍然是上節(jié)介紹的函數(shù)，只是在微分方程式的表達(dá)上做了一些變換, 轉(zhuǎn)換成一階微分微分方程的表達(dá)形式，然后求解。具體的運(yùn)用，我們?cè)谙旅娴睦渝榻B。例：6-22求解二階微分方程

46、d2x+4cost = 0dt2并畫出解在區(qū)間0 5上的圖像。程序設(shè)計(jì)與運(yùn)行結(jié)果：1. 我們先對(duì)微分方程進(jìn)行變換。設(shè)：x| = x, x2二文，于是微分方程可以表達(dá)成為一個(gè)微分方程組:f嚴(yán)x2x2 =-4 cost2. 在程序編輯器中建立待解微分方程表達(dá)式的函數(shù)文件funerx.m,并存入work工作目錄：function q=funerx(t,x)%x 有兩個(gè)元素 x(l)和 x(2)q二x(2);-4* cos(t);%q 有兩個(gè)列向量3. 在matlab的命令窗口求解，結(jié)果如圖6-7所示：>>cleart,x=ode45(funerx, 0 2*pi,4 0);% t,x =

47、ode45(funerx, 0:0.01:2*pi, 4 0);plot(t,x(:j)；-k t,x(:,2),:k,；linewidth, 2)title('fontsize18fontnametime new roman微分方程xm"+4cost=0求解圖像;'color','k') xlabel('fontsize 16)fontnametime new roman)t 'color；'k') ylabel('fontsize 16)fontname time new roman x',

48、'color','k') legend ('x=4cos(t)','v=-4sin(t)')axis (0 2*pi -4 4)grid on微分方程xn+4cost=0求解圖像圖67二階微分方程解的圖像1程序說明:1. t,x=ode45(funerx,02*pi,40),指定計(jì)算區(qū)間(t的范圍)為02*pi,指定初始值(x的兩 個(gè)元素)分別為4和0,實(shí)際上相當(dāng)于給定了初始位移是4 ,而初始速度為0。 t,x=ode45(funerx,0:0.01:2*pi,4 0),指定計(jì)算區(qū)間(t 的范圍)為0 2*pi,計(jì)算的步長為 0.

49、01 o 而在t,x=ode45(funerx,0 2*pi,4 0)中步長采用默認(rèn)值。2. plot(t,x(:,l)；-k t,x(:,2)；:k7linewidth 2),分別繪出t與x的第一列和第二列數(shù)據(jù)的圖像，第一 列第二列的圖像分別用實(shí)線和點(diǎn)線繪岀。3. 圖67中實(shí)線線表示解函數(shù)x=4cost在區(qū)間02*pi上的圖像，而另一條點(diǎn)線則是x的導(dǎo)數(shù)，即 v=-4sint在區(qū)間02*pi上的圖像。例：6-23求解二階微分方程x +(x2 +l)« + x =0并畫出解在區(qū)間0 5上的圖像。程序設(shè)計(jì)與運(yùn)行結(jié)果：1. 我們先對(duì)微分方程進(jìn)行變換。設(shè)：x| = x, x2 = * ,于是

50、微分方程可以表達(dá)成為一個(gè)微分方程組：x2=-x1-(x + l)x22. 在程序編輯器中建立待解微分方程表達(dá)式的函數(shù)文件funerxx.m,并存入work工作目錄：function q=funerxx(t,x)%x 有兩個(gè)元素 x(l)和 x(2)q=lx(2);- x(l)- (1+ x(l)a2)*x(2)j;%q 有兩個(gè)列向量3. 在matlab的命令窗口求解，結(jié)果如圖6-8所示：>>clearft,x=ode45(funerxxj0 12,1 1); plot(t,x(:,l),-k, t,x(:,2),':k,'linewidth', 2) tit

51、le('fontsize 18fontnametime new roman微分方程x,'"+(xa2+l)x"+x=0求解圖像；'color'/k')xlabel('fontsize 16fontnametime new roman t', 'color','k')ylabel('fontsize 16fontnametime new roman x', 'color','k') legend (,x,xh,)axis (0 12-0.6

52、1.4)微分方程xu+(x2+1)x'+x=0求解圖像圖6-8二階微分方程解的圖像grid on 程序說明：1. t,xj=ode45(funerxx,0 12,1 1),指定計(jì)算區(qū)間(t的范圍)為0 12,指定初始值(x的兩個(gè) 元素)分別為1和1。2. plot(t,x(:,l),k；t,x(:,2),*,linewidth； 2),分別繪出t與x的第一列和第二列數(shù)據(jù)的圖像，第一 列第二列的圖像分別用實(shí)線和點(diǎn)線繪出。3. 圖中實(shí)線線表示解函數(shù)x在區(qū)間0 12上的圖像，而另一條點(diǎn)線則是x的導(dǎo)數(shù)，即v在區(qū) 間0 12上的圖像。第七章matlab的符號(hào)計(jì)算第七章matlab的符號(hào)計(jì)算在上一

53、章我們知道了使用matlab對(duì)各種問題的數(shù)值求解方法。本章介紹matlab的符號(hào)運(yùn)算 功能。應(yīng)該說這一章更加精彩，它給出的不是具體的數(shù)值，而是一個(gè)含有字母變量的表達(dá)式。這是一種 精確計(jì)算，給出的答案是具有普遍意義的-般解。在系統(tǒng)講解matlab的符號(hào)運(yùn)算功能前，讓我們先通過兩個(gè)例子來看一下matlab的符號(hào)運(yùn)算 是怎么工作的。例：7-1求不定積分s = j sin xdx。我們知道這個(gè)不定積分的結(jié)果為x = -sin x。程序設(shè)計(jì)與運(yùn)行結(jié)果：»clearf='sin(x)'%定義被積函數(shù)f=int(f)%求不定積分f =-cos(x)例：7-2求不定積分r x(sin

54、 x)? +(cosx)2s = * dx oj x+sinx這個(gè)比較復(fù)雜，用手工計(jì)算就比較麻煩了。程序設(shè)計(jì)與運(yùn)行結(jié)果：»clearf二'(x*sin(x)人 3+cos(x)八 2)/(xa2+sin(x)'% 定義被積函數(shù)f=int(f)%求不定積分f =l/6*xa6+l/3*xa3+l/4*xa24-l/16*i*(i+2*x)*exp(i*x)a2+(l/2-3*i-3*x+3/2*i*xa2+l/2*xa3)*exp(i*x)+(l/2+3*i-3*x-3/2*i*xa2+l/2*xa3)/exp(i*x)-l/16*i*(-i+2*x)/exp(i*x)

55、a2+int(-2*exp(i*x)*(-l+xa4+xa7)/(2* xa2*exp(i*x)-i*exp(i*x)a2+i),x)第一個(gè)例子我們可以很容易地看出積分的結(jié)果，第二個(gè)例子計(jì)算起來就沒有那么簡單了，但是我們 看到用matlab計(jì)算起來依然和第一個(gè)一樣輕松。從這個(gè)例子中我們還可以看出，matlab的符號(hào)計(jì) 算應(yīng)用是如此的簡單，只需一個(gè)命令就解決問題了。當(dāng)然第二個(gè)例子給出的答案不是我們所熟悉的形式, 在下面的課程中我們將介紹如何給出符合我們書寫習(xí)慣的表達(dá)式。7-1符號(hào)表達(dá)式的生成7-1-1符號(hào)表達(dá)式的建立語法：f=sym(,expression,)說明：1. 在建立符號(hào)表達(dá)式時(shí)一定要注意不要隨意的添加空格