2013年最新matlab语言及其在工程领域中的应用ppt课件_1

第二章 MATLAB在数值分析中的应用,t = dat(:,1); CA1 = dat(:,2); CA2 = dat(:,3); plot(t,CA1,o,t,CA2,) axis(0 t(end)+0.5 min(min(dat(:,2:3)-0.1 max(max(dat(:,2:3)+0.05) xlabel(Time (min) ylabel(C_A_1, C_A_2 (kmol/m3) legend(C_A_1,C_A_2) title(Concentration profiles),问题的提出,在平面上给定n个点(xk,yk),可以唯一确定一个最多n-1次的多项式通过这些点，这个多项式叫插值多项式,插值和拟合,P(xk ) = yk , k = 1,2,n,Lagrange插值多项式,插值多项式例子,x0 =0:3; y 0= -5 6 1 16; disp(x0; y0),x0= 0 1 2 3 y0=-5 -6 -1 16,Polyinterp (Lagrange 插值形式),function y = polyinterp(x0,y0,x),n = length(x0); m=length(x);,z=x(i); s=0.0;,p=1.0;,if j=k,p=p*(z-x0 (j)/(x0(k)-x0(j);,end,s= s + p*y0(k);,y(i)=s;,Lagrange插值多项式基函数,最后一列全为1，倒数第二列 为一个指定的向量，其他各列 是其后列与倒数第二列的点乘 积，这种矩阵称为范得蒙(Vandermonde) 矩阵。在MATLAB中，可用函 数vander(V)生成以向量V为基 础向量生成。,x0 =0:3; y0 = -5 6 1 16; A = vander(x0); X = AB;,对A*X=B X=AB,x0 = 0 1 2 3 B =y0=-5 -6 -1 16 ,x= -.25:.01:3.25; y1=polyval(X,x) plot(x0,y0,o,x,y1),x0 =0:3; y0 = -5 -6 -1 16; x = -.25:.001:3.25; %- tic A = vander(x0); X= Ay0; y1=polyval(X,x); t1=toc %- tic y2 = polyinterp(x0,y0,x); t2=toc %- plot(x0,y0,o,x,y1,-,x,y2,b) legend(data,vand method,polyinterp),x00=data(:,1); y00=data(:,2); x0 =da00(:,1); y0 =da00(:,2); x =0:.1:9; %- tic A = vander(x0); X= Ay0; y1=polyval(X,x); t1=toc %- tic y2 = polyinterp(x0,y0,x); t2=toc %- plot(x00,y00,*,x0,y0,o,x,y1,-,x,y2,b) legend(da00,data,vand method,polyinterp),Polyinterp(symbol),symx =sym(x) P = polyinterp(x0,y0,symx) pretty(P) P = simplify(P) P = x3-2*x-5,simplify(S) simplifies each element of the symbolic matrix S,8/3*z*(x-1)*(x-2)+1/2*x*(x-1)*(x-3)-3*x*(x-2)*(x-3)-5/6*(-x+1)*(x-2)*(x-3),symx= sym(x) creates the symbolic variable with name x,pretty(S) prints the symbolic expression S in a format that resembles type-set mathematics.,Polyinterp(另外的例子),x0 = 1:6; y0 = 16 18 21 17 15 12; x = .75:.05:6.25; y = polyinterp(x0,y0,x); plot(x0,y0,o,x,y,-),x0=-5:1:5; y0=1./(1+x0.2); x=-5:0.1:5; y=polyinterp(x0,y0,x); y1=1./(1+x.2); plot(x0,y0,r-) hold on plot(x,y1,b-),Lagrange插值多项式在全区间内不收敛的情况,Hermite插值,设x1,x2,xn对应的函数值y1,y2, ,yn，一阶导数值y1,y2, ,yn,Hermite插值公式,function y=hermite(x0,y0,y1,x),n=length(x0); m=length(x);,yy=0.0;,h=h*(x(k)-x0(j)/(x0(i)-x0(j)2; a=1/(x0(i)-x0(j)+a;,y(k)=yy;,分段插值,一维插值：可以分为最近插值、线性插值、三次样条插值，分段三次Hermite插值。 y=interp1(x0,y0,x) y=interp1(x0,y0,x,method) method= nearst: 最近插值 linear: 线性插值（默认值） spline: 分段三次样条插值 pchip:分段三次Hermite插值,分段线性插值,x = 1:6; y = 16 18 21 17 15 12;,k(x = x0(j) = j;,n = length(x0); k = ones(size(x);,s = x x0(k);,y = y0(k) + s.*d(k);,% Evaluate interpolant,% x0小于x的最近位值,% Find subinterval indices, x0(k) = x x0(k+1),d = diff(y0)./diff(x0);,function y= piecelin(x0,y0,x),% First divided difference,x0=1:1:3; x=1:0.5:3;,x0=-5:1:5; y0=1./(1+x0.2); x=-5:0.1:5; y=interp1(x0,y0,x,linear); y1=1./(1+x.2); y2= piecelin(x0,y0,x) plot(x0,y0,) hold on plot(x,y,o,x,y1, r,x,y2,-),x0= 1:6; y0 = 16 18 21 17 15 12; x = .75:.05:6.25; y= interp1(x0,y0,x,pchip); plot(x0,y0,o,x,y,-);,三次样条,三次样条也是分段三次插值函数。可以给出光滑的插值曲线（面）。样条函数必须满足一阶、二阶导数连续。,xk,在三次样条中，要寻找三次多项式，以逼近每对数据点间的曲线。在样条术语中，这些数据点称之为断点。因为，两点只能决定一条直线，而在两点间的曲线可用无限多的三次多项式近似。因此，为使结果具有唯一性。在三次样条中，增加了三次多项式的约束条件。通过限定每个三次多项式的一阶和二阶导数，使其在断点处相等，就可以较好地确定所有内部三次多项式。此外，近似多项式通过这些断点的斜率和曲率是连续的。然而，第一个和最后一个三次多项式在第一个和最后一个断点以外，没有伴随多项式。因此必须通过其它方法确定其余的约束。最常用的方法，也是函数spline所采用的方法，就是采用非扭结(not-a-knot)条件。这个条件强迫第一个和第二个三次多项式的三阶导数相等。对最后一个和倒数第二个三次多项式也做同样地处理。 基于上述描述，人们可能猜想到，寻找三次样条多项式需要求解大量的线性方程。实际上，给定N个断点，就要寻找N-1个三次多项式，每个多项式有4个未知系数。这样，所求解的方程组包含有4*(N-1)个未知数。把每个三次多项式列成特殊形式，并且运用各种约束，通过求解N个具有N个未知系数的方程组，就能确定三次多项式。这样，如果有50个断点，就有50个具有50个未知系数的方程组。幸好，用稀疏矩阵，这些方程式能够简明地列出并求解，这就是函数spline所使用的计算未知系数的方法。,x = 1:6; y = 16 18 21 17 15 12; u = .75:.05:6.25; v = interp1(x,y,u,spline); plot(x,y,o,u,v,-);,x00=data(:,1); y00=data(:,3); x0 =da00(:,1); y0 =da00(:,3); x =0:.1:9; tic A = vander(x0); X= Ay0; y1=polyval(X,x); t(1)=toc; tic y2 = polyinterp(x0,y0,x); t(2)=toc; tic y3 =interp1(x0,y0,x,linear); t(3)=toc; tic,y4 =interp1(x0,y0,x,pchip); t(4)=toc; tic y5 =interp1(x0,y0,x,spline); t(5)=toc; %- t subplot(1,2,1); plot(x00,y00,*,x0,y0,o,x,y1,-,x,y2,b,x,y5,k-) legend(18da,interda,vand method,polyinterp,spline) subplot(1,2,2); plot(x00,y00,*,x0,y0,o,x,y3,-,x,y4,b,x,y5,k-) legend(18da,interda,linear,hermite,spline),t = 0.0013 0.0021 0.0030 0.0058 0.0092,例二: 反应器出口示踪剂浓度随时间的变化,t=linspace(0,5,100); y=1-cos(3*t).*exp(-t); plot(t,y,b); grid; hold on； plot(t,0.95*ones(size(t),r); hold off it=min(find(y0.95); %确定原始数据中第一次使y0.95的元素“下标” T=(it-3):(it+3); %在第一穿越点两侧，各取3点，作为插值“基准” t_nearst=interp1(y(T),t(T),0.95,nearst); t_linear=interp1(y(T),t(T),0.95, linear); t_cubic=interp1(y(T),t(T),0.95, pchip); t_spline=interp1(y(T),t(T),0.95,spline); disp( t_nearst t_linear t_pchip t_spline) disp(t_nearst t_linear t_pchip t_spline),t_nearst t_linear t_cubic t_spline 0.5051 0.4965 0.4962 0.4962,例三:,插值,二维插值：可以分为最近插值、线性插值、三次样条插值，立方插值。 z=interp2(x0,y0,z0,x,y) z=interp2(x0,y0,z0,x,y,method) method= nearst: 最近插值 linear: 线性插值（默认值） spline: 三次样条插值 cubic:立方插值,多项式拟合,p=polyfit(x,y,n) p,s=polyfit(x,y,n) y=polyval(p,x),p,s = polyfit(x,y,n) 返回多项式系数 P 和一个用于估计误差的结构数组 a.,y = polyval(p,x) 当 p 是一个长度为N+1的向量时,向量元素为一个多项式的系数,y为多项式在x处的值. y = p(1)*xN + p(2)*x(N-1) + . + p(N)*x + p(N+1),clear all;clc; x=0.5 1.0 1.5 2.0 2.5 3.0; y=1.75 2.45 3.81 4.80 8.00 8.60; a=polyfit(x,y,2); b=polyfit(x,y,4); c=polyfit(x,y,6); x1=0.5:0.05:3; y1=polyval(a,x1); y2=polyval(b,x1); y3=polyval(c,x1); plot(x,y,o,Markersize,10) hold on plot(x1,y1,r-,x1,y2,b,x1,y3,k) legend(data,n=2,n=4,n=6),x0=0:0.1:1; y0=-.447,1.978,3.11,5.25,5.02, 4.66,4.01,4.58,3.45,5.35,9.22; n=3; P=polyfit(x0,y0,n) xx=0:0.01:1; yy=polyval(P,xx); plot(xx,yy,-b,x0,y0,.r,MarkerSize,20) legend(拟合曲线,原始数据,Location,SouthEast) xlabel(x),clear all; clc xi = 1.00 1.25 1.50 1.75 2.00; yi = 5.10 5.79 6.53 7.45 8.46; Lnyi = log(yi); p = polyfit(xi,Lnyi,1) % 线性拟合：lny=bx+lna b = p(1) Lna = p(2); a = exp(Lna) x = linspace(xi(1),xi(end),100); Lny = polyval(p,x) plot(xi,Lnyi,o,x,Lny,-) xlabel(x) ylabel(lny),对非线性模型 进行最小二乘拟合 两边取对数：lnylna+bx,最小二乘法拟合生成样条曲线,平滑生成三次样条函数 sp=csaps(x,y,p) % sp:拟合得到的样条函数 ys=csaps(x,y,p,xx) x,y 要处理的离散数据，P平滑参数，取值为0,1, 取0时，相当于最小二乘直线拟合；取1时，相当于三次样条插值,csapi(),spline()。xx:指定在给定点xx上计算其三次样条函数值（ys)。,拟合,平滑生成B样条函数 sp=spap2(knots, k, x, y) % sp:拟合得到的样条函数 sp=spap2(knots, k, x, y) knots: 节点序数，k，样条函数的阶次，一般取4。 sp=spaps(x,y,tol) % sp:拟合得到的平滑样条函数 sp,ys=spaps(x,y,tol,m) tol 光滑时允许的精度，m 默认值为2，即平滑生成三次B样条曲线。 ys=fnval(sp,xx) 计算给定函数sp在定点xx处的函数值,% FermData.m % 青霉素发酵的实验数据 % t (h) 青霉素浓度（单位/ml） FermData = . 0 0 20 106 40 1600 60 3000 80 5810 100 8600 120 9430 140 10950 160 10280 180 9620 200 9400 ;,% FermDataFit.m % 发酵实验数据的最小二乘样条拟合和最小二乘多项式拟合 clear all；clc % 调入青霉素发酵的实验数据: % - FermData t = FermData(:,1); % t (h) C = FermData(:,2); % 青霉素浓度（单位/ml）,% 三次样条拟合-不经过实验点的样条函数 knots = 6; k = 4; m = 6; cs = spap2(knots,k,t,C); cs = spap2(newknt(cs),k,t,C); % for a better approximation (with the same number of pieces) by sp = spap2(newknt(sp),k,X,Y); % 求插值点ti处的青霉素浓度（单位/ml） ti = 10,30,50,70,90,110,130,150,170,190; Ci1 = fnval(cs,ti); % 计算时间间隔很小的Ct离散数据, 以便绘制光滑的拟合曲线 t4plot = t(1):1:t(end); C4plot1 = fnval(cs,t4plot);,% 最小二乘多项式拟合: p = polyfit(t,C,m); % 求插值点ti处的青霉素浓度（单位/ml） Ci2 = polyval(p,ti); C4plot2 = polyval(p,t4plot); % 计算结果及图形输出: fprintf(n在%s%sn,ti = 10,30,50,70,90,110,130,150,170,190,处的青霉素浓度为:) fprintf(%sn,(a) 最小二乘样条拟合方法:) fprintf( %.0f,Ci1) fprintf(n%sn,(b) 最小二乘多项式拟合方法:) fprintf( %.0f,Ci2) plot(t,C,ko,t4plot,C4plot1,r-,t4plot,C4plot2,b-) legend(实验点,最小二乘样条拟合,最小二乘多项式拟合) xlabel(时间 (h) ylabel(青霉素浓度 (单位/ml),Knots=2, k=4, m=4,Knots=6, k=4, m=6,数值积分与数值微分,基本思想：将被积函数在积分区间上离散节点(足够多)处的函数值与其所在小区间长度的乘积之和作为定积分的近似值或者构造一个函数P(x)逼近原函数，使得 数值积分方法主要有Newton-Cotes数值积分法，Gauss积分法和Romberg积分法。,数值积分与数值微分,Newton-Cotes（牛顿柯蒂斯）数值积分法,n=1，梯形公式 n=2，复化simpson公式 n=3，Cotes公式,n1 梯形公式(trapezoid rule),n2 抛物（simpson，辛普森）复化求积:,用抛物线弧段去逼近曲线f(x)，基本思想： 在积分区间xi,xi+1上，增加一个中点xi+1/2, 根据给定的插值条件xi，f(xi)和xi+1，f(xi+1)，构造一个二次插值求积多项式P2(x),Matlab中quad用递归自适应的外插Simpon公式 q=quad(fun,a,b) q=quad(fun,a,b,tol) q=quad(fun,a,b,tol,trace,p1,p2,) p1,p2：直接传递给函数fun的已知参数,fun=inline(exp(-x.*x),x); Isim=quad(fun,0,1) Isim = 0.74682418072642 1/2*erf(1)*pi(1/2) 0.74682413281243,inline(EXPR) ：constructs an inline function object from the MATLAB expression contained in the string EXPR.,数值积分与数值微分,自适应：在函数值变化较大 的部分减小步长,syms x Isym=vpa(int(exp(-x2),x,0,1) Isym = 0.74682413281242702539946743613185 format long d=0.001;x=0:d:1; Itrapz=d*trapz(exp(-x.*x) Itrapz = 0.74682407149919 fx=exp(-x.2); Ic=quad(fx,0,1,1e-8) Ic = 0.74682413285445,符号积分,梯形积分,quad8:：自适应的Cotes公式 q=quad8(fun,a,b) q=quad8(fun,a,b,tol),quadl:：自适应Lobatto求积函数 q=quadl(fun,a,b) q=quadl(fun,a,b,tol),Matlab7.1不 再包括该函数,建议优先 采用该函数,quad()： 0.84111302415845 quadl()： 0.84111708263095 syms x y y=1/(x(1/2)+x(1/3) s=int(y,0,1) S=5-6*log(2) 0.84111691664033,syms x y s=vpa(int(int(xy,x,0,1),y,1,2) s = 0.40546510810816438197801311546432 format long s_n=dblquad(x,y)x.y,0,1,1,2) s_n = 0.40546626724351,双重积分,Gauss求积公式,数值微分,利用离散数据求数值微分 有限差分法： dy=diff(y)./diff(x) 多项式拟合法： p=polyfit(x,y,n) 用离散数据求向量p表示的多项式拟合函数p pp=polyder(p) 对向量p表示的多项式函数进行求导，返回导函数为pp dy=polyval(pp,x) 求导函数pp在xi处的导数值 三次样条插值法： cs=csapi(x,y) 用离散数据求三次样条插值函数cs cp=fnder(cs)对三次样条函数cs进行求导，返回导函数为cp dy=fnval(pp,x) 求导函数cp在xi处的导数值 最小二乘法(样条拟合法)： sp=csaps(x,y,p) 用离散数据求样条拟合函数sp pp=fnder(sp)计算函数sp的dorder阶导数，默认为1阶 dy=fnval(pp,x)计算导函数pp在xi处的导数值,function xDifferential clear all；clc % 用函数y = cos(x)生成离散数据 x =0:pi/20:pi; y = cos(x); dy = -sin(x) % 准确值 % 有限差分法 dy1 = diff(y)./diff(x) x1 = x(1:length(x)-1); % 多项式拟合方法 p = polyfit(x,y,5); % 进行多项式拟合，得到多项式p pp = polyder(p); % 对拟合得到的多项式求导，得到导函数pp dy2 = polyval(pp,x) % 计算导函数pp在x的导数值 % 三次样条插值方法 cs = csapi(x,y); % 生成三次样条插值函数cs pp = fnder(cs); % 计算三次样条插值函数cs的导函数pp dy3 = fnval(pp,x) % 计算导函数pp在x的导数值,% 最小二乘样条拟合方法 sp = csaps(x,y,1); % 生成三次样条插值函数cs pp = fnder(sp); % 计算三次样条插值函数cs的导函数pp dy4 = fnval(pp,x) % 计算导函数pp在x的导数值 % 多项式拟合方法,三次样条插值方法的结果,最小二乘样条拟合方法与准确值之差 disp(多项式拟合方法和三次样条插值方法的结果与准确值之差dy-dy2;dy-dy3;dy-dy4：) disp(dy-dy2; dy-dy3;dy-dy4) % 绘制导函数曲线 plot(x,dy,ko,x1,dy1,k.-,x,dy2,b-,x,dy3,k-,x,dy4,r) xlabel(x) ylabel(y) title(函数cos(x)的导函数曲线) legend(准确值,有限差分法,多项式拟合方法,三次样条插值方法,样条拟合法),数值积分与数值微分,双组分简单精馏塔理论板数计算 填料吸收塔填料高度计算 理想间歇反应器、管式反应器、循环反应器（反应时 间、空时、管长） 纯气体的逸度系数 反应器停留时间分布,氯仿-苯双组分精馏系统，xf=0.4，xd=0.9，xw=0.15，精馏段、提馏段回流比为R=5，R=4，计算理论板数。 x = 0.178 0.275 0.372 0.456 0.650 0.844 y = 0.243 0.382 0.518 0.616 0.795 0.931,双组分简单精馏塔理论板数计算,精馏段理论板数,提馏段理论板数,function DistStagesCal clear all；clc xi = 0.178 0.275 0.372 0.456 0.650 0.844; yi = 0.243 0.382 0.518 0.616 0.795 0.931; xf = 0.4; xd = 0.9; xw = 0.15; R = 5; R1 = 4; sp = csaps(xi,yi,1); %根据离散数据拟合三次样条函数 x4plot = linspace(xi(1),xi(end),200); x在其区间等分200 y4plot = fnval(sp,x4plot); 计算等分区间节点处的样条函数值 plot(xi,yi,o,x4plot,y4plot,-) 绘制离散数据点和拟合函数曲线 N = quadl(func1,xf,xd,sp,xd,R); 数值积分，获得理论板数 M = quadl(func2,xw,xf,sp,xw,R1); N = round(N+0.5); 取最近的整数值 M = round(M+0.5); fprintf( 共需理论板数为：%d %sn,N+M,（块）) 定义计算精馏段理论板数的被积函数 function f = func1(x,sp,xd,R) y = fnval(sp,x); f = 1./(y-x-(xd-y)/R); 定义计算提馏段理论板数的被积函数 function f = func2(x,sp,xw,R1) y = fnval(sp,x); f = 1./(y-x-(y-xw)/R1);,气相反应A3R在等温管式反应器内进行，反应温度T= 215C，p=506.625kPa（5atm），系统保持恒压。反应速率为（-rA）=kCA1/2 。 k=0.01，反应物中加入惰性气体，摩尔比为1：1。出口转化率为0.8，计算反应器的空时。,应用实例一,function PFRSpaceTime clear all；clc k = 0.01; T = 488.15; % K yA0 = 0.5; % 进料气中A的摩尔分率 P = 5; % atm R = 0.082; % 理想气体常数 CA0 = yA0*P/(R*T); % mol/liter epsilon = 1; xAf = 0.8; % the interval: a, b, the step: d a = 0; b = xAf; d = (b-1)/100;,% 自适应Lobatto求积法(Adaptive Lobatto quadrature, high order) I = quadl(func,a,b,CA0,k,epsilon); tau = CA0*I; fprintf(tSpace time: %.1f (s),tau) function f = func(xA,CA0,k,epsilon) CA = CA0*(1-xA)./(1+epsilon*xA); rate = k*sqrt(CA); f = 1./rate;,反应物A在一等温间歇反应器内进行，A产物，测量得到的反应器中不同时间下反应物A的浓度CA（mol/L）如表，试确定反应速率方程。,动力学模型 -rA=kCAn ln(-rA)=nlnCA+lnk,用微分法进行动力学数据分析,多项式拟合法： p=polyfit(t,CA,n) 用离散数据求向量p表示的多项式拟合函数p pp=polyder(p) 对向量p表示的多项式函数进行求导，返回导函数为pp dCAdt=polyval(pp,t) 求导函数pp在xi处的导数值,function KineticDataFit clear all clc % 动力学数据 t = 0 20 40 60 120 180 300; CA = 10 8 6 5 3 2 1; % 用最小二乘样条拟合法计算微分dCA/dt-使用不经过实验点的B样条插值函数 knots = 3; K = 3; % 三次B样条 sp = spap2(knots,K,t,CA) sp = spap2(newknt(sp),K,t,CA); pp = fnder(sp) % 计算B样条函数的导函数 dCAdt = fnval(pp,t) % 计算t处的导函数值,% 绘制图形 ti = linspace(t(1),t(end),200); CAi = fnval(sp,ti) plot(t,CA,ro,ti,CAi,b-) xlabel(t) ylabel(C_A) figure fnplt(pp) % dCAdti = fnval(pp,ti) % plot(ti,dCAdti,-) xlabel(t) ylabel(dC/dt) % 线性拟合 rA = dCAdt; y = log(-rA); x = log(CA); p = polyfit(x,y,1); k = exp(p(2) n = p(1),在t=0时刻,在一容器入口处突然向流进容器的流体脉冲注入一定量的示踪剂,同时在容器出口处测量物料中示踪剂浓度随时间的变化,实验数据如下表:,(1) 平均停留时间,(2) 方差,(3) 扩散特征准数,function RTD1 % t1,t2 - Time, min % C1,C2 - Concentration, (kmol/m3)*1e3 % dt - Mean residence time, min % ss - Variance (sigma square), min2 % DL2uL1,DL2uL2 - Dispersion number % % Author: 000 % Copyright 2007 UNILAB Research Center, % % $Revision: 1.0 $ $Date: 2007/10/04 $ clear all clc t = 0:20:200; C = 0, 0, 0, 0, 0.4, 5.5, 16.2, 11.1, 1.7, 0.1, 0;,% 观察样条插值/拟合效果 % By spline(): plot(t,C,o) sp = spline(t,C) hold on fnplt(sp,b-) xlabel(Time (s) ylabel(C (kmol/m3)103) hold off % By spaps(): figure plot(t,C,o) sp = spaps(t,C,1) hold on,fnplt(sp,b-) xlabel(Time (s) ylabel(C (kmol/m3)103) hold off % Method 1: 直接用式（1）、（2）和（3）计算 t0 = 0; tf = t(end); IC = quadl(Func,t0,tf,sp); ICt = quadl(Func1,t0,tf,sp); ICt2 = quadl(Func2,t0,tf,sp); tm1 = ICt/IC ss1 = ICt2/IC - tm12 DL2uL1 = ss1/tm12/2,% Method 2: 用式（4）、（5）和（3）进行计算 tm2 = sum(C.*t)/sum(C) ss2 = sum(C.*t.2)/sum(C) - tm22 DL2uL2 = ss2/tm22/2 % - function y = Func(x,sp) % f= C y = fnval(sp,x); % 定义被积函数 % - function y = Func1(x,sp) % f= Ct y = fnval(sp,x).*x; % - function y = Func2(x,sp) % f= Ct2 y = fnval(sp,x).*x.2;,代数方程（组）数值解,线性方程（组）数值解 直接解法：矩阵除法 迭代解法：Jacobi迭代，Gauss-Seidel迭代，松弛迭代法,非线性方程（组）数值解 二分法： 迭代解法：Newton迭代，割线迭代法，Muller 迭代法,fzero求解单变量方程 x=fzero(fun,x0) x=fzero(fun,x0,options,p1,p2,),fsolve求解多变量非线性方程组 x=fsolve (fun,x0) x=fsolve (fun,x0,options,p1,p2,),线性方程（组）数值解直接解法：矩阵除法,A =1 2 4 2 4 1 4 6 7 b=6 15 24 x=Ab x =3.8571 1.9286 -0.4286,A * x = b,x = A b,线性方程（组）数值解迭代解法：Jacobi迭代,需要一个初始向量x0 要构造一个迭代式，可由x(k-1)算出x(k) 判断由迭代式计算的向量序列x(k)是否收敛 收敛终止的条件,迭代基本思想： Ax=b 构造一个等价方程组 x=Gx+d 迭代式为： x(k+1)=Gx(k)+d,令,方程可写为 xd + Gx,Jacobi迭代式：,x(k+1)=Gx(k)+d GD(L+U) d=Db D=diag(diag(A),A=10 -1 0 -1 10 -2 0 -2 10; b=9;7;6; U=triu(A,1) U = 0 -1 0 0 0 -2 0 0 0 L=tril(A,-1) L = 0 0 0 -1 0 0 0 -2 0,D=diag(diag(A) D = 10 0 0 0 10 0 0 0 10 G=D(L+U) G = 0 -0.1000 0 -0.1000 0 -0.2000 0 -0.2000 0 d=Db d= 0.9000 0.7000 0.6000,xd + Gx,function x=jacobi(A,b,x0) D=diag(diag(A); U=triu(A,1); %上三角变换 L=tril(A,-1); %下三角变换 G=D(L+U); d=Db; x=G*x0+d; n=1; while norm(y-x0)=1.0e-6 x0=x; x=G*x0+d; n=n+1; end x n,A=10 -1 0;-1 10 -2;0 -2 10; b=9;7;6; format long jacobi(A,b.0;0;0) x = 0.84842106875000 0.51578953125000 0.49684213750000 n = 11,在计算 时， 已经计算出来。而 常常比 更准确些，因此可以用 来代替 作进一步计算：,Gauss-Seidel(高斯赛德尔)迭代式,线性方程（组）数值解迭代解法： Gauss-Seidel迭代,x(k+1)=Gx(k)+d G=(D-L)-1U d=(D-L)-1b,A=10 -1 0 -1 10 -2 0 -2 10; b=9;7;6; U=triu(A,1) U = 0 -1 0 0 0 -2 0 0 0 L=tril(A,-1) L = 0 0 0 -1 0 0 0 -2 0,xd + Gx,D=diag(diag(A) D = 10 0 0 0 10 0 0 0 10 G=(D-L)U G = 0 -0.1000 0 0 0.0100 -0.2000 0 -0.0020 0.0400 d=(D-L)b d = 0.9000 0.6100 0.4780,function x=seidel(A,b,x0) D=diag(diag(A); U=triu(A,1); L=tril(A,-1); G=(D-L)U; d=(D-L)b; x=G*x0+d; n=1; while norm(x-x0)=1.0e-6 x0=x; x=G*x0+d; n=n+1; end x n,A=10 -1 0 -1 10 -2 0 -2 10; b=9;7;6; format long seidel(A,b,0;0;0) x = 0.84842104968750 0.51578947515625 0.49684210496875 n = 7,A=10 -1 0 -1 10 -2 0 -2 10; b=9;7;6; U=triu(A,1) U = 0 -1 0 0 0 -2 0 0 0 L=tril(A,-1) L = 0 0 0,xd + Gx,-1 0 0 0 -2 0 D=diag(diag(A) D = 10 0 0 0 10 0 0 0 10 G=(D-w*L)(1-w)*D+w*U) G = -0.1030 -0.1103 0 0.0114 -0.0908 -0.2206 -0.0025 0.0200 -0.0543 d=(D-w*L)b*w d = 0.9927 0.6626 0.5156,线性方程（组）数值解迭代解法：松弛迭代法,x(k+1)=Gx(k)+d G=(D-w*L)(1-w)*D+w*U) d=(D-w*L)b*w,function x=sor(A,b,w,x0) D=diag(diag(A); U=triu(a,1); L=tril(a,-1); G=(D-w*L)(1-w)*D+w*U); d=(D-w*L)b*w; x=G*x0+d; n=1; while norm(x-x0)=1.0e-6 x0=x; y=G*x0+d; n=n+1; end x n,A=10 -1 0 -1 10 -2 0 -2 10; b=9;7;6; format long sor(A,b,1.103,0;0;0) x = 0.84842104192112 0.51578946663881 0.49684209907452 n = 8,x 0.84842106875000 0.51578953125000 0.49684213750000 n = 11,x = 0.84842104968750 0.51578947515625 0.49684210496875 n = 7,Jacobi,Guass-Sediel,x= 0.84842104192112 0.51578946663881 0.49684209907452 n = 8,Sor,function y=erfen(fun,a,b,esp) if narginesp if feval(fun,a)*feval(fun,c)0 b=c; c=(a+b)/2; elseif feval(fun,c)*feval(fun,b)0 a=c; c=(a+b)/2; else,非线性方程（组）数值解：二分法,y=c; esp=10000; end n=n+1; end y=c; elseif feval(fun,a)=0 y=a; elseif feval(fun,b)=0 y=b; else disp(these, may not be a root) end n, fun=inline(x2-x-1,x) erfen(fun,0,10,0.05) n = 56 ans = 1.6180,function y=newton(x0) x1=- n=1; while (abs(x1