【正文】 %通過圖形來比較 2)0(,100,)1(4 22???? ???? ytt tyy11 ??? ty結(jié)果如下 數(shù)值解為藍(lán)色圓點(diǎn)，精確值為紅色實(shí)線，顯然兩者結(jié)果近似。r39。b.39。 %求數(shù)值解 y1=sqrt(t+1)+1。funt39。 y0=2。 t0=0。 function yp=funt(t,y) yp=(y^2t2)/4/(t+1)。 第六章 MATLAB數(shù)值計(jì)算 第六章 MATLAB數(shù)值計(jì)算 例 631 設(shè)有初值問題 求其數(shù)值解，并與精確解比較 (精確解為 )。 y0是初始狀態(tài)列向量。,tspan,y0) 其中 fname是定義 f(t,y)的函數(shù)文件名，該函數(shù)文件必須返回一個(gè)列向量。,tspan,y0) [t,y]=ode45(39。 第六章 MATLAB數(shù)值計(jì)算 Tttytfy ???? 0),( 00 )( yty ?第六章 MATLAB數(shù)值計(jì)算 龍格－庫(kù)塔法遞推公式 其中 )22(6)( 432110 kkkkhyyihty ii ??????? ?),()2,2()2,2(),(311421131112111hkyhtfkkhyhtfkkhyhtfkytfkiiiiiiii?????????????????? 龍格－庫(kù)塔法的實(shí)現(xiàn) 基于龍格－庫(kù)塔法， MATLAB提供了求常微分方程數(shù)值解的函數(shù)，一般調(diào)用格式為： [t,y]=ode23(39。一般采用等距節(jié)點(diǎn) tn=t0+nh, n=0,1,…,m, h稱為 步長(zhǎng) 。,x0,A,b,[],[],lb,[],[],[]) function f=fop(x) f=*x(2)+x(1)^2+x(2)^2x(1)*x(2)+1/30*x(1)^3 第六章 MATLAB數(shù)值計(jì)算 ? ? 2 2 32 1 2 1 21212121m in 30. 0 , 0f x x x x x x xs t x xxxxx? ? ? ? ???????? 數(shù)據(jù)處理與多項(xiàng)式計(jì)算 ? 數(shù)值微積分 ? 線性方程組求解 ? 非線性方程與最優(yōu)化問題 ?常微分方程的數(shù)值求解 常微分方程的數(shù)值求解 第六章 MATLAB數(shù)值計(jì)算 考慮常微分方程的初值問題 , 所謂其數(shù)值解法，就是求 它的解 y(t)在節(jié)點(diǎn)t0t1…t m處的近似值 y0,y1,…,y m的方法。] lb=[0,0] [x,f]=fmincon(39。] A=[1,。如果某個(gè)約束不存在，則用空矩陣表示。fname39。,0,5) xmin= xmax= fmin(‘mymin39。 xmin=fmin(39。 function fx=mymin(x) fx=x.^32*x5。 第六章 MATLAB數(shù)值計(jì)算 例 629 求 f(x)=x32x5在 [0,5]內(nèi)的最小值點(diǎn)與最大值點(diǎn)。 第六章 MATLAB數(shù)值計(jì)算 注意： 1. 當(dāng)目標(biāo)函數(shù)的階數(shù)大于 2時(shí)，使用 fminunc比fminsearch更有效。,x0) 這兩個(gè)函數(shù)的調(diào)用格式相似。,x1,x2) x=fminunc(39。,x1,x2) 求一元函數(shù)在區(qū)間 (x1,x2)中的最小值點(diǎn) x和最小值 fval. x=fminsearch(39。 第六章 MATLAB數(shù)值計(jì)算 無約束最優(yōu)化問題求解 MATLAB提供了 3個(gè)最小值函數(shù)，它們的調(diào)用格式為： [x,fval]=fminbnd(39。off39。Display39。,[,]39。 x=fsolve(39。 q(2)=*cos(x)+*sin(y)。 y=p(2)。 (1) 建立函數(shù)文件 。off39。Display39。,[,]39。 x=fsolve(39。 q(2)=*cos(x)+*sin(y)。 y=p(2)。 (1) 建立函數(shù)文件 。)將設(shè)定 Display選項(xiàng)為‘ off’。, 39。 optimset(39。如果想改變其中某個(gè)選項(xiàng)，則可以調(diào)用optimset()函數(shù)來完成。,X0,option) 其中 X為返回的解， fun是用于定義需求解的非線性方程組的函數(shù)文件名， X0是求根過程的初值，option為最優(yōu)化工具箱的選項(xiàng)設(shè)定。 fsolve函數(shù)的調(diào)用格式為： X=fsolve(39。funx39。 (2) 調(diào)用 fzero函數(shù)求根。 步驟如下： (1) 建立函數(shù)文件 。 tol控制結(jié)果的相對(duì)精度，缺省時(shí)取 tol=eps， trace 指定迭代信息是否在運(yùn)算中顯示，為 1時(shí)顯示，為 0時(shí)不顯示，缺省時(shí)取 trace=0。,x0,tol,trace) 其中 fname是待求根的函數(shù)文件名， x0為搜索的起點(diǎn)。該函數(shù)的調(diào)用格式為： z=fzero(39。0。0。6]。 b=[9。1,1,1。,) 第六章 MATLAB數(shù)值計(jì)算 例 625 分別用 Jacobi迭代和 GaussSerdel迭代法求解下列線性方程組，看是否收斂。 b=[9,7,6]39。1,10,2。 設(shè)迭代初值為 0，迭代精度為 106。 n=n+1。 %迭代次數(shù) while norm(yx0)=eps x0=y。 y=G*x0+f。 %求 A的上三角陣 G=(DL)\U。 %求 A的對(duì)角矩陣 L=tril(A,1)。 第六章 MATLAB數(shù)值計(jì)算 GaussSerdel迭代法的 MATLAB函數(shù)文件 ： function [y,n]=gauseidel(A,b,x0,eps) if nargin==3 eps=。,) ??????????????610272109102132121xxxxxxx第六章 MATLAB數(shù)值計(jì)算 2． GaussSerdel迭代法 在 Jacobi迭代過程中，計(jì)算 xi(k+1)時(shí)， x1(k+1), … ,xi1(k+1) ,已經(jīng)得到 ,不必再用 x1(k), … ,xi1(k) ，即原來的迭代公式 Dx(k+1)=(L+U)x(k)+b可以 改進(jìn)為 Dx(k+1)=Lx(k+1)+Ux(k)+b，于是得到： x(k+1)=(DL)1Ux(k)+(DL)1b 該式即為 GaussSerdel迭代公式。 b=[9,7,6]39。1,10,2。 設(shè)迭代初值為 0，迭代精度為 106。 n=n+1。 %迭代次數(shù) while norm(yx0)=eps x0=y。 y=B*x0+f。 %求 A的上三角陣 B=D\(L+U)。 %求 A的對(duì)角矩陣 L=tril(A,1)。 第六章 MATLAB數(shù)值計(jì)算 Jacobi迭代法的 MATLAB函數(shù)文件 ： function [y,n]=jacobi(A,b,x0,eps) if nargin==3 eps=。 第六章 MATLAB數(shù)值計(jì)算 1． Jacobi迭代法 對(duì)于線性方程組 Ax=b，如果 A為非奇異方陣，即 aii≠0(i=1,2,… ,n)，則可將 A分解為 A=DLU，其中 D為對(duì)角陣，其元素為 A的對(duì)角元素， L與 U為 A的下三角陣和上三角陣，于是 Ax=b化為： x=D1(L+U)x+D1b 與之對(duì)應(yīng)的迭代公式為： x(k+1) =D1(L+U)x(k