1、基本科學運算基本科學運算MATLAB方程求根方程求根n多項式的求根x2-x-1=0 roots(1 -1 -1)ans = -0.6180 1.6180求下列三階線性代數(shù)方程組的近似解5426255452321321321xxxxxxxxxMATLAB程序為：A=2 -5 4;1 5 -2;-1 2 4；b=5;6;5；x=Ab 或A=2 -5 4;1 5 -2;-1 2 4A=inv(A)b=5;6;5x=A*b線性方程組的數(shù)值解法 常微分方程的數(shù)值解法 Euler算法演示 微分方程 時刻的值為 迭代公式 這樣就能把數(shù)值解一步一步計算出來MATLAB下的常微分方程求解函數(shù) 主求解函數(shù) 其他可

2、用的求解函數(shù)，調(diào)用格式一致 ode23、ode15s()、ode113() options可以由odeget和odeset函數(shù)處理 微分方程可以用下面的方式之一描述 匿名函數(shù) M-函數(shù) inline函數(shù)：不建議使用常微分方程的數(shù)值解法n數(shù)值解指令ode45ode23n格式t,y=ode23(f,t0,tf,y0,options)t,y=ode45(f,t0,tf,y0,options)f :微分方程的形式，需要先編寫微分方程的形式，需要先編寫m文件文件t0,tf：自變量的初值和終值y0：函數(shù)的初值options：積分參數(shù)，指定誤差例一 y=-y+x+1, y(0)=1首先編寫方程的M文件fun

3、ction z=ff1(x,y)z=-y+x+1ode23(ff1,0 1,1);grid00.10.20.30.40.50.60.70.80.9111.051.11.151.21.251.31.351.4例二：已知一個剛體在無外力作用下運動方程及初始條件為 仿真時間為0,120(0)y1,(0)y(0)y 2 y123213312321yyyyyyyy編寫m文件： f=(t,y)y(2)*y(3);-y(2)*y(3);-2*y(1)*y(2); tspan=0 12 y0=0 1 1 t,y=ode45(f,tspan,y0) plot(t,y(:,1),r,t,y(:,2),b*,t,y

4、(:,3),k-)024681012-1-0.500.511.5常用options選項常微分方程舉例 Lorenz方程 初值 描述微分方程，然后求解繪圖f=(t,x)-8*x(1)/3+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3)t,x=ode45(f,0,100,0,0,1e-10)plot(t,x(:,1),r,t,x(:,2),k-,t,x(:,3),b*);gridxlabel(t);ylabel(x)0102030405060708090100-30-20-1001020304050txVan der Pol方程求解 Van d

5、er Pol方程 選擇狀態(tài)變量 帶有附加參數(shù)m，可以求解f=(t,x,mu)x(2);-mu*(x(1)2-1)*x(2)-x(1)h_opt=odeset;x0=-0.2;-0.7;t_final=20;mu=1;t1,y1=ode45(f,0 t_final,x0,h_opt,mu)mu=6;t2,y2=ode45(f,0 t_final,x0,h_opt,mu)plot(t1,y1,r,t2,y2,k-);grid02468101214161820-10-8-6-4-20246810隱式微分方程求解舉例 隱式微分方程 可以寫出真正的剛性微分方程求解 常微分方程 初值微分方程組的變換和技巧

6、 解決其他形式微分方程到標準型 的變換方法 單個高階常微分方程處理方法 選擇一組狀態(tài)變量 原微分方程可以變換為高階常微分方程組的變換方法 一般高階微分方程組 選擇狀態(tài)變量Apollo衛(wèi)星軌跡 衛(wèi)星方程 選擇狀態(tài)變量 微分方程變換成 直接求解 直接求解結(jié)果是錯誤的，應該檢驗 微分代數(shù)方程的標準化 求解 所謂符號計算是指在運算時,無須事先對變量賦值,而將所得到結(jié)果以標準的符號形式來表示。 MathWorks公司以Maple的內(nèi)核作為符號計算引擎（Engine），依賴Maple已有的函數(shù)庫，開發(fā)了實現(xiàn)符號計算的兩個工具箱：基本符號工具箱和擴展符號工具箱。MATLAB的符號計算 參與符號運算的對象可以

7、是符號變量、符號表達式或符號矩陣。符號變量要先定義，后引用?？梢杂胹ym函數(shù)、syms函數(shù)將運算量定義為符號型數(shù)據(jù)。引用符號運算函數(shù)時，用戶可以指定函數(shù)執(zhí)行過程中的變量參數(shù)；若用戶沒有指定變量參數(shù)，則使用findsym函數(shù)默認的變量作為函數(shù)的變量參數(shù)。（一）（一） 定義符號變量定義符號變量1、sym函數(shù) sym函數(shù)的主要功能是創(chuàng)建符號變量，以便進行符號運算，也可以用于創(chuàng)建符號表達式或符號矩陣。用sym函數(shù)創(chuàng)建符號變量的一般格式為： x = sym(x)其目的是將x創(chuàng)建為符號變量，以x作為輸出變量名。每次調(diào)用該函數(shù),可以定義一個符號變量。（一） 定義符號變量【例1】作符號計算：a,b,x,y均為

8、符號運算量。在符號運算前，應先將a,b,x,y定義為符號運算量15byaxbyax（一） 定義符號變量a=sym(a); %定義a為符號運算量，輸出變量名為ay =2/bb=sym(b);x=sym(x);y=sym(y”); x,y=solve(a*x-b*y-1,a*x+b*y-5,x,y) %以a,b為符號常數(shù)，x,y為符號變量即可得到方程組的解：x =3/ay =2/b（一） 定義符號變量2、syms函數(shù)syms函數(shù)的功能與sym函數(shù)類似。syms函數(shù)可以在一個語句中同時定義多個符號變量，其一般格式為： syms arg1 arg2 argN 用于將rg1, arg2,argN等符號創(chuàng)

9、建為符號型數(shù)據(jù)。（一） 定義符號變量 在數(shù)學表達式中，一般習慣于使用排在字母表中前面的字母作為變量的系數(shù)，而用排在后面的字母表示變量。例如： f=ax2+bx+c 表達式中的a,b,c通常被認為是常數(shù)，用作變量的系數(shù)；而將x看作自變量。 （二）默認符號變量（二）默認符號變量例如，數(shù)學表達式 f=xn g=sin(at+b)根據(jù)數(shù)學式中表示自變量的習慣，默認a,b,c為符號常數(shù)，x為符號變量。若在MATLAB中表示上述表達式，首先用syms 函數(shù)定義a，b，n，t，x為符號對象。在進行導數(shù)運算時，由于沒有指定符號變量，則系統(tǒng)采用數(shù)學習慣來確定表達式中的自變量，默認a,b,c為符號常數(shù)，x，t為符

10、號變量。即 ： 對函數(shù)f求導為：df/dx 對函數(shù)g求導為：dg/dt（二）默認符號變量 符號表達式由符號變量、函數(shù)、算術(shù)運算符等組成。符號表達式的書寫格式與數(shù)值表達式相同。例如,數(shù)學表達式 其符號表達式為： 1+sqr(5*x)/2 注意，在定義表達式前應先將表達式中的字符x定義為符號變量。251x（三） 符號表達式（四） 生成符號函數(shù) 將表達式中的自變量定義為符號變量后，賦值給符號函數(shù)名，即可生成符號函數(shù)。例如有一數(shù)學表達式：222),(cbyaxyxf 其用符號表達式生成符號函數(shù)fxy的過程為： syms a b c x y %定義符號運算量 fxy=(a*x2+b*y2)/c2 %生成

11、符號函數(shù) 生成符號函數(shù)fxy后，即可用于微積分等符號計算。（四） 生成符號函數(shù)【例4】定義一個符號函數(shù) fxy=(a*x2+b*y2)/c2 ，分別求該函數(shù)對x、y的導數(shù)和對x的積分。syms a b c x y %定義符號變量fxy=(a*x2+b*y2)/c2； %生成符號函數(shù) diff(fxy，x) %符號函數(shù)fxy對x求導數(shù)ans =2*a*x/c2diff(fxy， y) %符號函數(shù)fxy對y求導數(shù) ans =2*b*y/c2 %符號函數(shù)fxy對x求積分int(fxy， x) ans =1/c2*(1/3*a*x3+b*y2*x)（四） 生成符號函數(shù)（一） 微積分函數(shù)1.求極限 函數(shù)

12、limit用于求符號函數(shù)f的極限。系統(tǒng)可以根據(jù)用戶要求，計算變量從不同方向趨近于指定值的極限值。該函數(shù)的格式及功能：微積分 limit(f,x,a)：求符號函數(shù)f（x）的極限值。即計算當變量x趨近于常數(shù)a時，f（x）函數(shù)的極限值。 limit(f,a)：求符號函數(shù)f（x）的極限值。由于沒有指定符號函數(shù)f（x）的自變量，則使用該格式時，符號函數(shù)f（x）的變量為函數(shù)findsym(f)確定的默認自變量，既變量x趨近于a。 limit(f)：求符號函數(shù)f（x）的極限值。符號函數(shù)f（x）的變量為函數(shù)findsym(f)確定的默認變量；沒有指定變量的目標值時，系統(tǒng)默認變量趨近于0，即a=0的情況。 li

13、mit(f,x,a,right)：求符號函數(shù)f的極限值。right表示變量x從右邊趨近于a。 limit(f,x,a,left)：求符號函數(shù)f的極限值。left表示變量x從左邊趨近于a。二、微積分syms x； %定義符號變量f=(x*(exp(sin(x)+1)-2*(exp(tan(x)-1)/sin(x)3； %確定符號表達式w=limit(f) %求函數(shù)的極限w = -1/2xeextgxxx3sin0sin) 1(2) 1(lim例 求極限2. 微分函數(shù)diff函數(shù)用于對符號表達式s求微分。該函數(shù)的一般引用格式為： diff(s,v,n)說明： 應用diff（s）沒有指定微分變量和微

14、分階數(shù)，則系統(tǒng)按findsym函數(shù)指示的默認變量對符號表達式s求一階微分。 應用diff（s，v）或diff（s，sym（v） 格式，表示以v為自變量，對符號表達式s求一階微分。 應用diff（s，n）格式，表示對符號表達式s求n階微分，n為正整數(shù)。 應用diff（s，v，n）diff（s，n，v） 格式，表示以v為自變量，對符號表達式s求n階微分。x = sym(x); %定義符號變量t = sym(t);diff(sin(x2) %求導運算ans =2*cos(x2)*xdxxd2sin例 求導數(shù)3積分函數(shù)積分函數(shù)int（s ，v，a，b)可以對被積函數(shù)或符號表達式s求積分。其引用格式為：

15、 int（s ，v，a，b)應用int（s）格式，表示沒有指定積分變量和積分階數(shù)時，系統(tǒng)按findsym函數(shù)指示的默認變量對被積函數(shù)或符號表達式s求一階積分。應用int（s，v）格式，表示以v為自變量，對被積函數(shù)或符號表達式s求一階不定積分。應用積分函數(shù)時，如果給定 a、b兩項，表示是進行定積分運算。a、b分別表示定積分的下限和上限。不指定積分的下限和上限表示求不定積分。syms xint(1/(1+x2) ans =atan(x)dxx211求積分求積分4. 級數(shù)(級數(shù)求和)級數(shù)求和運算是數(shù)學中常見的一種運算。例如： f(x)=a0+a1x+a2x2+a3x3+anxn函數(shù)symsum可以用

16、于此類對符號函數(shù)f的求和運算。該函數(shù)的引用時，應確定級數(shù)的通項式s，變量的變化范圍a和b。該函數(shù)的引用格式為： symsum(s, a,b)求級數(shù)的和:鍵入：1/12+1/22+1/32+1/42+ syms k symsum(1/k2,1,Inf) %k值為1到無窮大ans =1/6*pi2其結(jié)果為：1/12+1/22+1/32+1/42+ =2/6dxdtydydtx ,S=dsolve(Dx=y,Dy=-x);disp(blanks(12),x,blanks(21),y),disp(S.x,S.y) x ycos(t)*C1+sin(t)*C2, -sin(t)*C1+cos(t)*C2

17、 微分方程符號解微分方程符號解dsolve(eqn1,eqn2)MATLAB解插值擬合問題解插值擬合問題 插插 值值用用MATLABMATLAB作插值計算作插值計算一維插值函數(shù)一維插值函數(shù)：yi=interp1(x，y，xi，method)插值方法插值方法被插值點被插值點插值節(jié)點插值節(jié)點xixi處的插處的插值結(jié)果值結(jié)果nearest ：最鄰近插值：最鄰近插值linear ： 線性插值；線性插值；spline ： 三次樣條插值；三次樣條插值； cubic ： 立方插值。立方插值。缺省時：缺省時： 分段線性插值。分段線性插值。注意：所有的插值方法都要求注意：所有的插值方法都要求x x是單調(diào)的，是單

18、調(diào)的， 并且并且xi不能夠不能夠 超過超過x的范圍。的范圍。splinespline 函數(shù)函數(shù)：yi=spline(x，y，xi)被插值點被插值點插值節(jié)點插值節(jié)點xixi處的插處的插值結(jié)果值結(jié)果pp=spline(x，y)通過插值節(jié)點通過插值節(jié)點的分段三次多的分段三次多項式項式xy機翼下輪廓線例例 已知飛機下輪廓線上數(shù)據(jù)如下，求已知飛機下輪廓線上數(shù)據(jù)如下，求x每改變每改變0.1時的時的y值。值。x=0,3,5,7,9,11,12,13,14,15;x=0,3,5,7,9,11,12,13,14,15;y=0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6;y=0,1.2

19、,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6;xi=0:0.1:15;xi=0:0.1:15; yi=interp1(x,y,xi,spline);plot(xi,yi,k);grid05101500.511.522.5 要求要求x0,y0 x0,y0單調(diào)；單調(diào)；x x，y y可取可取為矩陣，或為矩陣，或x x取行向量，取行向量，y y取為列向量，取為列向量，x,yx,y的值分別不能超出的值分別不能超出x0,y0 x0,y0的范圍。的范圍。z=interp2(x0,y0,z0,x,y,method)被插值點插值方法用用MATLAB作網(wǎng)格節(jié)點數(shù)據(jù)的插值作網(wǎng)格節(jié)點數(shù)據(jù)的插值插值

20、節(jié)點被插值點的函數(shù)值nearest nearest 最鄰近插值最鄰近插值linear linear 雙線性插值雙線性插值cubic cubic 雙三次插值雙三次插值缺省時缺省時, , 雙線性插值雙線性插值例：測得平板表面例：測得平板表面3 3* *5 5網(wǎng)格點處的溫度分別為：網(wǎng)格點處的溫度分別為： 82 81 80 82 84 82 81 80 82 84 79 63 61 65 81 79 63 61 65 81 84 8484 84 82 85 86 82 85 86 試作出平板表面的溫度分布曲面試作出平板表面的溫度分布曲面z=f(x,y)z=f(x,y)的圖形。的圖形。輸入以下命令：x=

21、1:5;y=1:3;temps=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;mesh(x,y,temps)1.先在三維坐標畫出原始數(shù)據(jù)，畫出粗糙的溫度分布曲圖.2以平滑數(shù)據(jù),在x、y方向上每隔0.2個單位的地方進行插值.再輸入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)畫出插值后的溫度分布曲面圖. 1234511.522.53606570758085901234511.522.5360657075808590例例 山區(qū)地貌：山區(qū)地貌： 在某山區(qū)

22、測得一些地點的高程如下表。平面區(qū)域為在某山區(qū)測得一些地點的高程如下表。平面區(qū)域為 1200=x=4000,1200=y=3600)試作出該山區(qū)的地貌圖和等高線圖，并對幾種插值方法進行比較。試作出該山區(qū)的地貌圖和等高線圖，并對幾種插值方法進行比較。 X Y12001600200024002800320036004000120011301250128012301040900500700160013201450142014001300700900850200013901500150014009001100106095024001500120011001350145012001150101028001

23、500120011001550160015501380107032001500155016001550160016001600155036001480150015501510143013001200980 通過此例對最近鄰點插值、雙線性插值方法和雙三次插值方法的插值效果進行比較。To MATLAB (moutain)編寫m文件如下：clcx=1200:400:4000y=1200:400:3600z=1130 1250 1280 1230 1040 900 500 700; 1320 1450 1420 1400 1300 700 900 850; 1390 1500 1500 1400 90

24、0 1100 1060 950; 1500 1200 1100 1350 1450 1200 1150 1010; 1500 1200 1200 1550 1600 1550 1380 1070; 1500 1550 1620 1550 1600 1600 1600 1550; 1480 1500 1550 1510 1430 1300 1200 980mesh(x,y,z)xi=1200:10:4000yi=1200:10:3600zi=interp2(x,y,z,xi,yi,nearest);020004000600002000400060000500100015002000XYZ 擬擬

25、合合 用用MATLAB解擬合問題解擬合問題p = polyfit (x，y，n）多項式次多項式次數(shù)數(shù)已知節(jié)點已知節(jié)點擬合多項擬合多項式系數(shù)向式系數(shù)向量量多項式曲線擬合多項式曲線擬合即要求即要求 出二次多項式出二次多項式:3221)(axaxaxf中中 的的),(321aaaA 使得使得:最小 )(1112iiiyxf例例 對下面一組數(shù)據(jù)作二次多項式擬合對下面一組數(shù)據(jù)作二次多項式擬合xi0.10.20.40.50.60.70.80.91yi1.9783.286.167.347.669.589.489.3011.200.20.40.60.81-20246810121）輸入以下命令：）輸入以下命令：

26、 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ù)點和擬合曲線的圖形作出數(shù)據(jù)點和擬合曲線的圖形2）計算結(jié)果：）計算結(jié)果： = -9.8108 20.1293 -0.0317解解: 用多項式擬合的命令用多項式擬合的命令00.20.40.60.81-20246810120317.01293.208108.9)(2xxxf已知數(shù)據(jù)點數(shù)據(jù)點： xdataxdata= =（xdata1，xdata

27、2，xdataxdatan n），）， ydataydata= =（ydataydata1 1，ydataydata2 2，ydataydatan n） 用用MATLAB作作非線性非線性最小二乘擬合最小二乘擬合 MatlabMatlab的提供了求非線性最小二乘擬合的函數(shù)：的提供了求非線性最小二乘擬合的函數(shù)：最小 ),(21niiiydataxdataxF lsqcurvefitlsqcurvefit用以求含參量用以求含參量x x（向量）的向量值函數(shù)（向量）的向量值函數(shù)F(x,xdataF(x,xdata)=)=（F F（x x，xdataxdata1 1），），F(xiàn) F（x x，xdataxda

28、tan n）T T中的參變量中的參變量x(x(向量向量),),使得使得 lsqcurvefitx = lsqcurvefit (fun,x0,xdata,ydata); fun是一個事先建立的是一個事先建立的定義函數(shù)定義函數(shù)F(x,xdata), 自自變量為變量為x和和xdata x = lsqcurvefit (fun,x0,xdata,ydata);迭代初值迭代初值已知數(shù)據(jù)點已知數(shù)據(jù)點 輸入格式為輸入格式為:擬合函數(shù)的參數(shù)擬合函數(shù)的參數(shù)100200 30040050060070080090010004.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6

29、.59jt310jc210102. 0),(minjjktcbeakbaFj 例例 用下面一組數(shù)據(jù)擬合用下面一組數(shù)據(jù)擬合 中的參數(shù)中的參數(shù)a，b，kktbeatc2 . 0 . 0)(該問題即解最優(yōu)化問題：該問題即解最優(yōu)化問題： 1 1）編寫）編寫M-M-文件文件 curvefun1.mcurvefun1.m function f=curvefun1(x,tdata) f=x(1)+x(2)*exp(-0.02*x(3)*tdata) %其中其中 x(1)=a; x(2)=b；x(3)=k;2）輸入命令）輸入命令tdata=100:100:1000tdata=100:100:1000cdata

30、=cdata=1e-03* *4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59; x0=0.2,0.05,0.05; x0=0.2,0.05,0.05; x=lsqcurvefitx=lsqcurvefit (curvefun1,x0,tdata,cdata) (curvefun1,x0,tdata,cdata) f= f= curvefun1(x,tdata) F(x，tdata)= ，x=(a，b，k)Tktktbeabea),(10102.

31、 002. 0解法解法1 1. 用命令用命令lsqcurvefitlsqcurvefit3 3）運算結(jié)果為）運算結(jié)果為：f =0.0043 0.0051 0.0056 0.0059 0.0061 f =0.0043 0.0051 0.0056 0.0059 0.0061 0.0062 0.0062 0.0063 0.0063 0.0063 0.0062 0.0062 0.0063 0.0063 0.0063 x = 0.0063 -0.0034 0.2542 x = 0.0063 -0.0034 0.25424）結(jié)論）結(jié)論：a=0.0063, b=-0.0034, k=0.2542MATLAB(fzxec