matlab與科學(xué)計算簡單易學(xué),老師的ppt,很不錯.ppt

matlab科學(xué)工程計算,同濟大學(xué)數(shù)學(xué)系 陳雄達(dá) ,2,teapotdemo,3,目錄,matlab入門 matlab基礎(chǔ)編程 matlab作圖 matlab工具箱 matlab上機操作 matlab綜合訓(xùn)練 matlab命令速查,4,matlab入門,matlab是什么 matlab不是什么 matlab集成環(huán)境 matlab ABC,5,matlab入門,沒有matlab就沒有樂趣。 M.N. Nachtigal S.C. Reddy L.N. Trefethen 關(guān)于迭代法的Copper Mountain 論文集 1990,6,matlab是什么,一個可視化的計算程序, 廣泛使用于從個人計算機到超級計算機范圍內(nèi)的各種計算機上 包括命令控制、可編程，上百個預(yù)先定義命令和函數(shù) 有許多強有力的命令, 能完成大量的高級矩陣處理 強有力的二維、三維圖形工具 能與其他程序一起使用 25個(不斷增加中)不同的工具箱應(yīng)用于特殊的應(yīng)用領(lǐng)域 工業(yè)研究與開發(fā)的有力工具 數(shù)學(xué)教學(xué), 尤其線代, 數(shù)值分析, 科學(xué)計算方面的教研工具 電子學(xué), 控制理論, 物理學(xué)等工程科學(xué)方面的教研工具 經(jīng)濟學(xué), 化學(xué)和生物學(xué)等有計算問題的所有領(lǐng)域中的教學(xué)與研究 名字取自矩陣實驗室(matrix laboratory),7,matlab不是什么,不是萬能的解決工具 不是最高性能的編程語言 受計算條件限制, 不能解決超大型實際問題 不能解決工具箱之外的問題種類 - 需要編寫接口、算法甚至工具箱,8,集成環(huán)境,窗口系統(tǒng) View Desktop Layout Five Panel 歷史命令 / 變量和文件管理 / 命令窗口 菜單系統(tǒng) File / Edit / View / Web / Window / Help 按鈕,9,matlab ABC,提示符 注釋符 % 續(xù)行 . 123 + 45 - 67 + 8 9 x = 3 * sin(pi/4) 2 向量(數(shù)組) v =1 3 5 2 4 6; 矩陣(二維數(shù)組) A =1 3; 5,2 4 6;,10,matlab基礎(chǔ)編程,數(shù)據(jù)結(jié)構(gòu) 冒號(:) 矩陣操作入門 標(biāo)識符 初等函數(shù)(elfun) 初等函數(shù)(exp) 函數(shù)fix,round, ceil,floor 邏輯判斷 關(guān)系運算,運算的級別 matlab幫助 結(jié)構(gòu)化編程 腳本文件 函數(shù)文件,11,數(shù)據(jù)結(jié)構(gòu),最基本的數(shù)據(jù)結(jié)構(gòu)- 矩陣 數(shù)和向量看成為特殊的矩陣 矩陣以 為定界符，與多維數(shù)組等同 字符串看成為每個元素都是單個字符的向量，也可以有字符矩陣 高維數(shù)組 細(xì)胞 (cell) 結(jié)構(gòu)體 (struct),12,冒號(:),冒號 a:s:b 從a開始, 步長為s, 界為b 例如 1:2:10 1 3 5 7 9 7:-2:0 7 5 3 1 1:6 1 2 3 4 5 6 6:2 (空矩陣),13,標(biāo)識符,文件名、變量名、函數(shù)名 標(biāo)識符的規(guī)定 最長不超過19個字符 為配合matlab的風(fēng)格，采用見名知意的小寫短名稱 以字母開頭，包含數(shù)字、大小寫字母、下劃線 變量命名規(guī)則 常用從簡，專用從繁 矩陣大寫 A, B, C 向量小寫 u, v, w, x, y, z 函數(shù)小寫 f, g, h, fun, f1, f2 常量小寫 alpha, beta, a, b, c,14,內(nèi)部變量,pi 圓周率=3.1415926 i,j 虛根 inf 無限大 eps 2.2204460e-16 NaN 不定型(Not a Number) 內(nèi)部變量的運算規(guī)則 inf 參與的運算 NaN 參與的運算,15,點運算(向量運算), a = 1 2 3 4; b = 0 1 2 3; a + b ans = 1 3 5 7 b./a ans = 0 0.5000 0.6667 0.7500 b.*a ans = 0 2 6 12 a.b ans = 1 2 9 64,16,矩陣操作入門,矩陣的抽取 A = 8 1 6; 3 5 7; 4 9 2 ; A(2,3) 得到 7 A(2, 3 2 1) 得到 7 5 3 A(3 2,2:-1:1) 得到 9 4; 5 3 A(2,:) 得到 3 5 7 A(:,end) 得到 6; 7; 2,17,賦值語句,A(2)= A(n-1,:)= A(2+3, 2*3)= x= A(sqrt(4),x-3)= A(:)= v(end)=,A(2) + 2 = x2 + x = A(sqrt(4),x)-3 = A = inv(A) = (x=0) =,左值表達(dá)式: 指向某一存儲空間的表達(dá)式，放在賦值語句左邊。,18,初等函數(shù)(elfun),三角函數(shù)(sin,sinh,asin,asinh) 指數(shù)對數(shù)(exp,log,log10,sqrt) 復(fù)變函數(shù)(abs,conj,real,imag) 取整取余(floor,ceil,round,mod) 全部都具有向量功能: 輸入向量, 返回對應(yīng)函數(shù)值組成的向量,19,初等函數(shù)(exp), x = 0:0.2:1 x = 0 0.2 0.4 0.6 0.8 1.0 y = exp(x) y = 1.00 1.22 1.49 1.82 2.22 2.71,向量功能,20,函數(shù)fix,round,ceil,floor,21,簡單矩陣操作,rand(n), eye(n), zeros(n) 生成 n 階隨機、單位、零矩陣, rand(2) ans = 0.9501 0.6068 0.2311 0.4860 rand(3,2) ans = 0.8913 0.0185 0.7621 0.8214 0.4564 0.4447, eye(2,3) ans = 1 0 0 0 1 0 zeros(3) ans = 0 0 0 0 0 0 0 0 0,22,diag, diag(1 2 3; 4 5 6 7 8 9) ans = 1 5 9 diag(1 2 3; 4 5 6 7 8 9,1) ans = 2 6, diag(1 5 9) ans = 1 0 0 0 5 0 0 0 9 diag(4 8,-1) ans = 0 0 0 4 0 0 0 8 0,23,triu/tril, A = magic(3) A = 8 1 6 3 5 7 4 9 2 U = triu(A,-1) U = 8 1 6 3 5 7 0 9 2, L1 = tril(A,-1) L1 = 0 0 0 3 0 0 4 9 0 L2 = tril(A) L2 = 8 0 0 3 5 0 4 9 2,24,字符運算, char(4,1,8,8,11+100) ans = hello abs(matlab) ans = 109 97 116 108 97 98 eval(exp(-pi*i) ans = -1 + 0i, feval(sin,pi/2) ans = 1 fun = cos; feval(fun,pi/2) ans = 0,字符串,25,字符運算,for k = 1:9, st = fun . num2str(k); x(k) = . feval(st, x(k); end function y = fun1(x) function y = fun2(x),str2num(2*pi) ans = 6.2832 str2num(cos(2*pi) ans = 1 eval(cos(2*pi) ans = 1,26,邏輯判斷,真值 邏輯真(true,1)或假(false,0) 邏輯運算 , =, =, =, = 例如 又如, pi = sqrt(10) ans = 0, 1:5 = pi ans = 0 0 0 1 1,向量功能,27,關(guān)系運算,關(guān)系運算有： 與(and) and(a,b) a & b 或(or) or(a,b) a | b 非(not) not(a) a 異或(xor) xor(a,b),28,兩個特殊關(guān)系運算(any,all),僅當(dāng)輸入向量 v 每個分量為真(非零值), all(v)返回1；否則返回0。 僅當(dāng)輸入向量 v 每個分量為假(零), any(v)返回0；否則返回1。 any(0 1 2 3) ans = 1 all(0 1 2 3) ans = 0,29,運算的級別,函數(shù)、括號 算術(shù)運算 邏輯運算 關(guān)系運算 每一個運算級別內(nèi)仍有優(yōu)先級的高低,30,matlab幫助,help help 命令名/工具箱名 help help doc doc eig what 列出當(dāng)前matlab文件 who/whos 列出當(dāng)前變量 which 尋找matlab文件或命令所在目錄 lookfor 比較 lookfor curve / help curve,31,結(jié)構(gòu)化編程,分支結(jié)構(gòu)(if) for循環(huán) while循環(huán) 嵌套結(jié)構(gòu) 程序規(guī)范,32,分支結(jié)構(gòu)(if),if 判斷條件1 執(zhí)行語句1 elseif 判斷條件2 執(zhí)行語句2 elseif 判斷條件3 執(zhí)行語句3 else 執(zhí)行語句n+1 end,if temp 10, disp(cold!); elseif temp 20, disp(cool); elseif temp 27, disp(warm); elseif temp 31, disp(hot); else disp(too hot); end,33,for循環(huán),for 循環(huán)變量 = 循環(huán)值 循環(huán)操作 end,s = 0; for k = 1:100, s = s + 1/k2; end,for k = hello world, fprintf(%c,k); pause(0.1); end,for k = magic(3), k+1 end,向量功能,34,while循環(huán),while 循環(huán)條件 循環(huán)操作 end while 1, if 循環(huán)條件, break; end 循環(huán)操作 end,n = input(n= ); while n=1, if mod(n,2)=1, n = n*3 + 1; else n = n / 2; end end,Collatz猜想,35,嵌套結(jié)構(gòu),for k = 1:9, for j = 1:9, fprintf(%d*%d=%3d , k,j,k*j); end fprintf(n); end,九九乘法表,36,嵌套結(jié)構(gòu),n = input(n=0: ); flag = 0; while flag, nn = n; rn = 0; while nn=0, d = mod(nn,10); rn = rn * 10 + d; nn = floor(nn/10); end disp(rn) if rn=n, flag = 1; else n = rn + n; end end,輸入n 設(shè)置迭代指標(biāo)為0 當(dāng)?shù)怀晒r(=0), 保留當(dāng)前n 計算n的逆序為rn 顯示rn 重新設(shè)置迭代指標(biāo) 或 n = rn + n end %迭代不成功,偽代碼,源代碼,對一個整數(shù)n, 把它反過 來讀的數(shù)稱為逆序數(shù)。 如132的逆序數(shù)為231。 寫一個算法，輸入n, 若 它與自己的逆序數(shù)相等, 則停止；否則把它與自 己的逆序數(shù)相加，重新 判斷。,37,程序規(guī)范,程序按照結(jié)構(gòu)縮進(jìn)(matlab會自行處理) 按照邏輯關(guān)系分段 添加適當(dāng)?shù)淖⑨?5%-30%) 規(guī)范所有的標(biāo)識符(文件名、變量名) 適當(dāng)增加空行、空格以保證可讀性 辦法: 參考Matlab內(nèi)部程序風(fēng)格,38,腳本文件,n = input(years: ); for k = n, if mod(k,400)=0, disp(num2str(k) is a leap year); elseif mod(k,100)=0, disp(num2str(k) isnt a leap year); elseif mod(k,4)=0, disp(num2str(k) is a leap year); else disp(num2str(k) isnt a leap year); end end,判別閏年不用 這么煩吧？,n = input(years: ); for k = n, if mod(k,400)=0 | . (mod(k,4)=0 end end,好像簡單 多了.,39,函數(shù)文件,function p,t = fun(n) % n: input integer % p: sum of factors % t: number of factors p = 0; t = 0; for k = 1:n, if mod(n,k)=0, p = p + k; t = t + 1; end end,function p,t= fun(n) % n: input integer % p: sum of factors % t: number of factors k = 1:n; z = mod(n,k); t = sum(z=0); p = sum(k(z=0);,保存為fun.m,40,函數(shù)的調(diào)用, z,p = fun(15) z = 24 p = 4 x = fun(6) x = 12 fun(28) ans = 56, p,z = fun Input argument n is undefined. p,z = fun(20,18) Too many input arguments. p,t,b = fun(15) Too many output arguments.,function p,t = fun(n),41,函數(shù)的相互調(diào)用,尋找親和數(shù)對, 即這樣兩個數(shù), 它們的真因子之和等于對方。,function X = amic(n) s = 1; for k = 2:n, if k=fun(fun(k), X(s,1 2) = k, fun(k); s = s + 1; end end,function p = fun(n) k = 1:n-1; z = mod(n,k); p = sum(k(z=0);,調(diào)用,子函數(shù),42,遞歸調(diào)用,函數(shù)調(diào)用自己稱為遞歸因此必須有停止自調(diào)用的時刻。,例：Frayer級數(shù), 即0,1上分母不超過n的既約真分?jǐn)?shù)從小到大排列。,每一次把新出現(xiàn)的分?jǐn)?shù), 放在兩個已有分?jǐn)?shù)中, 這兩個 分?jǐn)?shù)的分子和分母各自的和就是新分?jǐn)?shù)的分子和分母。,43,遞歸調(diào)用,function p,q = frayer(n) if n1, p,q = frayer(n-1); ind = find(q(1:end-1)+q(2:end)=n); for k = ind(end:-1:1), q = q(1:k) n q(k+1:end); p = p(1:k) p(k)+p(k+1) p(k+1:end); end elseif n=1, p = 0,1; q = 1,1; end,44,遞歸調(diào)用,function v = fact(n) if n=0|n=1, v = 1; else v = n * fact(n-1); end, v = prod(1:n);,階乘,45,算法的復(fù)雜度,s = 0; for k = 1:n, p = 1; for j = 1:k, p = p * j; end s = s + 1/(p+1); end,op = 0; for k = 1:n, op = op + 1; for j = 1:k, op = op + 1; end op = op + 3; end,op=4n+n(n+1)/2,46,算法的復(fù)雜度,s = 0; p = 1; for k = 1:n, p = p*k; s = s + 1/(p+1); end,op = 2; for k = 1:n, op = op + 1; op = op + 3; end,op=4n+2,47,matlab作圖,二維作圖 三維作圖 其它作圖,48,二維作圖(plot),plot(x,y,s) 例如 x = linspace(0,2*pi); y = sin(x); plot(x,y,r-);,49,二維作圖(plot),plot(x,y,s) - s的規(guī)定,50,二維作圖(plot),xlabel(x); ylabel(y); title(y=sin(x);,51,二維作圖(polar), theta = linspace(0,6*pi); rho = theta; polar(theta,rho,r-);,52,三維作圖(plot3),plot3(x,y,z,s) t = 0:pi/50:10*pi; plot3(sin(t),cos(t),t);,53,三維作圖(comet3), t = 0:pi/500:20*pi; comet3(sin(t),cos(t),t);,54,三維作圖(meshgrid),X,Y = meshgrid(x,y) 例如 X，Y = meshgrid(1:3, 4:7) 則 X =1 2 3 且 Y = 4 4 4 1 2 3 5 5 5 1 2 3 6 6 6 1 2 3 7 7 7,55,三維作圖(mesh),x = linspace(-5,5); X,Y = meshgrid(x); Z = X.2 Y.2; mesh(X,Y,Z); view(-60,28);,如果X,Y由meshgrid(x,y) 生成，可寫 mesh(x,y,Z) surf有同樣的格式,向量功能,56,其它作圖(bar3), A = 82 82 63 81 82 81 84 61 84 85 80 79 65 84 86; bar3(A,grouped);,57,其它作圖(pie3),pie3(2 4 3 5,0 1 1 0,. North,South,East,West),58,surf, contour,例: 山區(qū)地貌- 在某山區(qū)測得一些地點的高程如下表。平面區(qū)域為 1200=x=4000,1200=y=3600)。試作出該山區(qū)的地貌圖和等高線圖。,59,surf, contour, A = 1130 1250 1280 1230 1040 900 500 700 1320 1450 1420 1400 1300 700 900 850 1390 1500 1500 1400 900 1100 1060 950 1500 1200 1100 1350 1450 1200 1150 1010 1500 1200 1100 1550 1600 1550 1380 1070 1500 1550 1600 1550 1600 1600 1600 1550 1480 1500 1550 1510 1430 1300 1200 980; x = 1200:400:4000; y = 1200:400:3600; surf(x,y,A); contour(x,y,A,100);,60,matlab工具箱(toolbox),矩陣計算(help matfun) 最優(yōu)化工具箱(help optim) 微分方程工具箱(help ode45,help pde) 多項式計算(help polyfun) 統(tǒng)計工具箱(help stats),61,矩陣計算(inv),inv(A) 求矩陣A的逆 例 A = magic(3) A = 8 1 6 3 5 7 4 9 2 inv(A) ans = 0.1472 -0.1444 0.0639 -0.0611 0.0222 0.1056 -0.0194 0.1889 -0.1028,62,矩陣計算(), A = hilb(3) A = 1.0000 0.5000 0.3333 0.5000 0.3333 0.2500 0.3333 0.2500 0.2000 b = ones(3,1); % 3行1列全是1的向量 x = inv(A) * b x = 3.0000 -24.0000 30.0000, A = hilb(3) A = 1.0000 0.5000 0.3333 0.5000 0.3333 0.2500 0.3333 0.2500 0.2000 b = ones(3,1); % 3行1列全是1的向量 x = A b x = 3.0000 -24.0000 30.0000,63,矩陣的秩(rank), A = magic(4) A = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 rank(A) ans = 3,64,矩陣的行列式(det), A = pascal(4) A = 1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20 det(A) ans = 1,65,矩陣的跡(trace), A = 1:4 * 1:4 A = 1 2 3 4 2 4 6 8 3 6 9 12 4 8 12 16 trace(A) ans = 30,66,矩陣的Cholesky分解(chol), A = pascal(4) A = 1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20 B = chol(A) B = 1 1 1 1 0 1 2 3 0 0 1 3 0 0 0 1,67,矩陣的特征值(eig), A = gallery(3) A = -149 -50 -154 537 180 546 -27 -9 -25 eig(A) ans = 1.00000000001122 1.99999999999162 2.99999999999700, V,D = eig(A) V = 0.3162 -0.4041 -0.1391 -0.9487 0.9091 0.9740 -0.0000 0.1010 -0.1789 D = 1.0000 0 0 0 2.0000 0 0 0 3.0000,68,最優(yōu)化工具箱(linprog),什么是線性規(guī)劃？, f = 1 1; A = -2 -1; -5 -8; b = -10 -40; Aeq = ; beq = ; lb = 0 4; ub = 7 inf; linprog(f,A,b,Aeq,beq,lb,ub) Optimization terminated successfully. ans = 3.0000 4.0000,69,最優(yōu)化工具箱(fminu),x,opts=fminu(f,x0,options,g,p1,p2,),最優(yōu)解 向量,求解 情況,工具箱 函數(shù),目標(biāo)函數(shù) 字符串,初值,算法參數(shù),目標(biāo)函數(shù) 梯度 字符串,目標(biāo)函數(shù) 其他參數(shù),help foptions,70,最優(yōu)化工具箱(fminu),function v=f(x) v = 100*( x(2)-x(1)2 )2 + . (x(1)-1)2; x = fminu(f,0 0) x = 1.0000 1.0000,71,最優(yōu)化工具箱(fzero),x,fv=fzero(f,x0,options,p1,p2,),零點向量,函數(shù)值,工具箱 函數(shù),目標(biāo)函數(shù) 字符串,初值,算法參數(shù),目標(biāo)函數(shù) 其他參數(shù),help foptions,72,最優(yōu)化工具箱(fzero), f = inline(x+exp(x)-3); x,fv=fzero(f,2) x = 0.79205996843068 fv = 0 只能求單變量函數(shù)的根！,73,最優(yōu)化工具箱(fmincon),x = fmincon(fun,x0,. 1 1 1,8,. 1 2 0,6,. 0 0 0,. cons. ); function v = fun(x) v =-3*x(1)-x(2)2+x(3); function c,ceq = cons(x) c(1)= x*x-16; c(2)= x(2)*(x(1)+x(3)-12; ceq = ;,最優(yōu)解 6 0 0,74,微分方程工具箱(ode45),t,y = ode45(odefun,tspan,y0,options,p1,p2.),求解函數(shù),工具箱 函數(shù),算法參數(shù),help odeset,求解區(qū)間,初值,解的向量 可用plot(x,y)畫出解,目標(biāo)函數(shù) 其他參數(shù),75,微分方程工具箱(ode45),function dydx =f(x,y) dydx = y x + 1; x,y = . ode45(f,0,2,1); plot(x,y,r-);,76,多項式計算,多項式表示：向量,p = 1 0 2 表示,多項式求值 polyval(p,1:5) 多項式求根 r = roots(p) 用根生成多項式 p = poly(r),多項式相乘 p = 1 0 2; q = 2 1; s = conv(p,q); 多項式相除 q,r = deconv(s,p); 多項式求導(dǎo) d = polyder(p),77,多項式計算,多項式擬合(polyfit) x = linspace(0,pi,8) x = 0 0.4488 0.8976 1.3464 1.7952 2.2440 2.6928 3.1416 y = sin(x) y = 0 0.4339 0.7818 0.9749 0.9749 0.7818 0.4339 0.0000 p = polyfit(x,y,2) p = -0.4011 1.2601 -0.0179 xx = linspace(0,pi,200); yy = polyval(p,xx); plot(xx,yy,r-,x,y,b-,linewidth,2),78,插值(interp1),yi=interp1(x，y，xi，method),nearest: 最鄰近插值linear: 線性插值 spline: 三次樣條插值 cubic: 立方插值 缺省時: 分段線性插值,79,插值(interp2),z=interp2(x0,y0,z0,x,y,method),nearest 最鄰近插值 linear 雙線性插值 cubic 雙三次插值 缺省時, 雙線性插值,80,統(tǒng)計工具箱,參數(shù)估計 betafit, binofit, mle, 概率密度函數(shù) betapdf, binopdf, hygepdf, 各種分布的均值、方差 betastat, binostat, chi2stat, 線性、非線性的各種回歸模型 anova1, anova2, polyf