(完整版)MATLAB数值分析实例
MATLAB应用实例分析
Matlab 应用例题选讲仅举一些运用MATLAB 的例子，这些问题在数学建模中时常遇到，希望能帮助同学们在短时间内方便、快捷的使用MATLAB 解决数学建模中的问题，并善用这一工具。
常用控制命令：clc ：%清屏； clear ：%清变量； save ：%保存变量； load ：%导入变量一、利用公式直接进行赋值计算本金P 以每年n 次，每次i%的增值率（n 与i 的乘积为每年增值额的百分比）增加，当增加到r ×P 时所花费的时间T 为：（利用复利计息公式可得到下式）)01.01ln(ln )01.01(i n rT i P P r nT +=⇒+=⨯（12,5.0,2===n i r ）MATLAB 的表达形式及结果如下： >> r=2;i=0.5;n=12; %变量赋值 >> T=log(r)/(n*log(1+0.01*i)) 计算结果显示为：T = 11.5813即所花费的时间为T=11.5813 年。
分析：上面的问题是一个利用公式直接进行赋值计算问题，实际中若变量在某个范围变化取很多值时，使用MATLAB ，将倍感方便，轻松得到结果，其绘图功能还能将结果轻松的显示出来，变量之间的变化规律将一目了然。
若r 在[1,9]变化，i 在[0.5,3.5]变化；我们将MATLAB 的表达式作如下改动，结果如图1。
r=1:0.5:9; i=0.5:0.5:3.5; n=12;p=1./(n*log(1+0.01*i)); T=log(r')*p; plot(r,T)xlabel('r') %给x 轴加标题 ylabel('T') %给y 轴加标题q=ones(1,length(i));text(7*q-0.2,[T(14,1:5)+0.5,T(14,6)-0.1,T(14,7)-0.9],num2str(i'))rT图1从图1中既可以看到T 随r 的变化规律，而且还能看到i 的不同取值对T —r 曲线的影响（图中的六条曲线分别代表i 的不同取值）。
#数值分析Matlab作业
数值分析编程作业2012年12月第二章14.考虑梯形电阻电路的设计，电路如下：电路中的各个电流{i1,i2,…,i8}须满足下列线性方程组：121232343454565676787822/252025202520252025202520250i i V R i i i i i i i i i i i i i i i i i i i i -=-+-=-+-=-+-=-+-=-+-=-+-=-+=这是一个三对角方程组。
设V=220V ，R=27Ω，运用追赶法，求各段电路的电流量。
Matlab 程序如下：function chase () %追赶法求梯形电路中各段的电流量 a=input('请输入下主对角线向量a='); b=input('请输入主对角线向量b='); c=input('请输入上主对角线向量c='); d=input('请输入右端向量d='); n=input('请输入系数矩阵维数n='); u(1)=b(1); for i=2:nl(i)=a(i)/u(i-1); u(i)=b(i)-c(i-1)*l(i); endy(1)=d(1); for i=2:ny(i)=d(i)-l(i)*y(i-1); endx(n)=y(n)/u(n); i=n-1; while i>0x(i)=(y(i)-c(i)*x(i+1))/u(i); i=i-1;end x输入如下: >> chase请输入下主对角线向量a=[0,-2,-2,-2,-2,-2,-2,-2]; 请输入主对角线向量b=[2,5,5,5,5,5,5,5];请输入上主对角线向量c=[-2,-2,-2,-2,-2,-2,-2,0]; 请输入方程组右端向量d=[220/27,0,0,0,0,0,0,0]; 请输入系数矩阵阶数n=8 运行结果如下：x = 8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477第三章14.试分别用(1)Jacobi 迭代法；(2)Gauss-Seidel 迭代法解线性方程组1234510123412191232721735143231211743511512x x x x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=--⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥---⎣⎦⎣⎦⎣⎦ 迭代初始向量(0)(0,0,0,0,0)T x =。
(完整word版)含答案《MATLAB实用教程》
第二章 MATLAB 语言及应用实验项目实验一 MATLAB 数值计算三、实验内容与步骤1.创建矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=987654321a（1（2）用（3）用（42.矩阵的运算（1）利用矩阵除法解线性方程组。
⎪⎪⎩⎪⎪⎨⎧=+++=-+-=+++=+-12224732258232432143214321421x x x x x x x x x x x x x x x 将方程表示为AX=B ，计算X=A\B 。
（2）利用矩阵的基本运算求解矩阵方程。
已知矩阵A 和B 满足关系式A -1BA=6A+BA ，计算矩阵B 。
其中⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=7/10004/10003/1A ，Ps: format rata=[1/3 0 0;0 1/4 0;0 0 1/7];b=inv(a)*inv(inv(a)-eye(3))*6*a（3）计算矩阵的特征值和特征向量。
已知矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=1104152021X ，计算其特征值和特征向量。
（4）Page:322利用数学函数进行矩阵运算。
已知传递函数G(s)=1/(2s+1),计算幅频特性Lw=-20lg(1)2(2w )和相频特性Fw=-arctan(2w),w 的范围为[0.01，10]，按对数均匀分布。
3.多项式的运算（1）多项式的运算。
已知表达式G(x)=(x-4)(x+5)(x 2-6x+9),展开多项式形式，并计算当x 在[0，20]内变化时G(x)的值，计算出G(x)=0的根。
Page 324（2）多项式的拟合与插值。
将多项式G(x)=x 4-5x 3-17x 2+129x-180，当x 在[0，20]多项式的值上下加上随机数的偏差构成y1，对y1进行拟合。
对G(x)和y1分别进行插值，计算在5.5处的值。
Page 325 四、思考练习题1.使用logspace 函数创建0～4π的行向量，有20个元素，查看其元素分布情况。
Ps: logspace(log10(0),log10(4*pi),20) (2) sort(c,2) %顺序排列 3.1多项式1）f(x)=2x 2+3x+5x+8用向量表示该多项式，并计算f(10)值. 2）根据多项式的根[-0.5 -3+4i -3-4i]创建多项式。
数值分析matlab程序实例
1，秦九韶算法，求出P(x=3)=2+4x+5x^2+2x^3的值clear all;x=3;n=3;a(1)=2;a(2)=4;a(3)=5;a(4)=2v(1)=a(n+1);for k=2:(n+1);v(k)=x*v(k-1)+a(n-k+2);endp=v(n+1)p=,1132，一次线型插值程序：利用100.121.求115的开方。
clear all;x1=100;x2=121;y1=10;y2=11;x=115;l1=(x-x2)/(x1-x2);l2=(x-x1)/(x2-x1);p1=l1*y1+l2*y2p1=10.71433，分段插值程序，已知为S1（x)为（0，0），（1，1），（2，5）（3,8）上的分段一次插值，求S1（1.5）.clear allx=[0123];y=[0158];n=length(x);a=1.5;for i=2:nif(x(i-1)<=a<x(i));endendH1=y(i-1)+(y(i)-y(i-1))/(x(i)-x(i-1))*(a-x(i-1))H1=3.50004）曲线拟合：用一个5次多项式在区间[0，2π]内逼近函数sin(x)。
clear allX=linspace(0,2*pi,50);Y=sin(X);[P,S]=polyfit(X,Y,5)plot(X,Y,'k*',X,polyval(P,X),'k-')P=-0.00560.0874-0.39460.26850.87970.0102S=R:[6x6double]df:44normr:0.03375)求有理分式的导数clear allP=[3,5,0,-8,1,-5];Q=[10,5,0,0,6,0,0,7,-1,0,-100];[p,q]=polyder(P,Q)6）将以下数据按从小到大排序：4.3 5.7 5.2 1.89.4a=[4.35.75.21.89.4];b(1:100)=0;n=1;b(a*10)=1;for k=1:100a(n)=k/10;if b(k)>0a(n)=k/10;n=n+1;endendaa=1.8000 4.3000 5.2000 5.70009.400010.00007)用二分法求方程x 3-x-1=0在[1,2]内的近似根，要求误差不超过10-3。
第2讲 matlab的数值分析
第二讲MATLAB的数值分析2-1矩阵运算与数组运算矩阵运算和数组运算是MATLAB数值运算的两大类型，矩阵运算是按矩阵的运算规则进行的，而数组运算则是按数组元素逐一进行的。
因此，在进行某些运算（如乘、除）时，矩阵运算和数组运算有着较大的差别。
在MATLAB中，可以对矩阵进行数组运算，这时是把矩阵视为数组，运算按数组的运算规则。
也可以对数组进行矩阵运算，这时是把数组视为矩阵，运算按矩阵的运算规则进行。
1、矩阵加减与数组加减矩阵加减与数组加减运算效果一致，运算符也相同，可分为两种情况：（1）若参与运算的两矩阵（数组）的维数相同，则加减运算的结果是将两矩阵的对应元素进行加减，如A=[1 1 1;2 2 2;3 3 3];B=A;A+Bans=2 2 24 4 46 6 6（2）若参与运算的两矩阵之一为标量（1*1的矩阵），则加减运算的结果是将矩阵（数组）的每一元素与该标量逐一相加减，如A=[1 1 1;2 2 2;3 3 3];A+2ans=3 3 34 4 45 5 52、矩阵乘与数组乘（1）矩阵乘矩阵乘与数组乘有着较大差别，运算结果也完全不同。
矩阵乘的运算符为“*”，运算是按矩阵的乘法规则进行，即参与乘运算的两矩阵的内维必须相同。
设A、B为参与乘运算的=A m×k B k×n。
因此，参与运两矩阵，C为A和B的矩阵乘的结果，则它们必须满足关系C m×n算的两矩阵的顺序不能任意调换，因为A*B和B*A计算结果很可能是完全不一样的。
如：A=[1 1 1;2 2 2;3 3 3];B=A;A*Bans=6 6 612 12 1218 18 18F=ones(1,3);G=ones(3,1);F*Gans3G*Fans=1 1 11 1 11 1 1（2）数组乘数组乘的运算符为“.*”，运算符中的点号不能遗漏，也不能随意加空格符。
参加数组乘运算的两数组的大小必须相等（即同维数组）。
数值分析实验— MATLAB实现
数值分析实验——MATLAB实现姓名sumnat学号2013326600000班级13级应用数学2班指导老师2016年1月一、插值：拉格朗日插值 (1)1、代码： (1)2、示例： (1)二、函数逼近：最佳平方逼近 (2)1、代码： (2)2、示例： (2)三、数值积分：非反常积分的Romberg算法 (3)1、代码： (3)2、示例： (4)四、数值微分：5点法 (5)1、代码： (5)2、示例： (6)五、常微分方程：四阶龙格库塔及Adams加速法 (6)1、代码：四阶龙格库塔 (6)2、示例： (7)3、代码：Adams加速法 (7)4、示例： (8)六、方程求根：Aitken 迭代 (8)1、代码： (8)2、示例： (9)七、线性方程组直接法：三角分解 (9)1、代码： (9)2、示例： (10)八、线性方程组迭代法：Jacobi法及G-S法 (11)1、代码：Jacobi法 (11)2、示例： (12)3、代码：G-S法 (12)4、示例： (12)九、矩阵的特征值及特征向量：幂法 (13)1、代码： (13)2、示例： (13)一、插值：拉格朗日插值1、代码：function z=LGIP(x,y)%拉格朗日插值n=size(x);n=n(2);%计算点的个数syms a;u=0;%拉格朗日多项式f=1;%插值基函数for i=1:nfor j=1:nif j==if=f;elsef=f*(a-x(j))/(x(i)-x(j));endendu=u+y(i)*f;f=1;endz=expand(u);%展开2、示例：>> x=1:6;y1=x.^5+3*x.^2-6;y2=sin(x)+sqrt(x);>> f1=LGIP(x,y1)f1 =-6+3*a^2+a^5%可知多项式吻合得很好>> f2=vpa(LGIP(x,y2),3)f2 =.962e-1*a^4+1.38*a+.300*a^2+.504-.436*a^3-.616e-2*a^5二、函数逼近：最佳平方逼近1、代码：function z=BestF(u,a,b,n)%最佳平方逼近，用x^i逼近,n为逼近的次数n=n+1;syms xreal;old=findsym(u);u=subs(u,old,x); %将u中变量替换为xf=sym('');H=sym('');d=sym('');for i=1:n %生成函数系f(1,i)=x^(i-1);endfor i=1:n %生成内积Hfor j=1:nH(i,j)=int(f(1,i)*f(1,j),a,b);endendfor i=1:n %生成内积dd(i,1)=int(f(1,i)*u,a,b);enda=H\d;%解法方程Ha=dz=a'*f';2、示例：>> syms x>> f1=sqrt(x);>> f2=x^5+x^2;>> f3=exp(x);>> a=0 ;b=1;>> BestF(f1,a,b,5)ans =12/143+420/143*x-1120/143*x^2+2016/143*x^3-1800/143*x^4+56/13*x^5>> BestF(f2,a,b,5)ans =x^5+x^2>> BestF(f3,a,b,5)ans =-566826+208524*exp(1)+(16733010-6155730*exp(1))*x+(-115830120+42611520*exp(1))* x^2+(306348840-112699440*exp(1))*x^3+(-342469260+125987400*exp(1))*x^4+(136302012-5 0142708*exp(1))*x^5>> vpa(ans,3)ans =.1e4-.1e6*x-.1e7*x^3+.1e7*x^4三、数值积分：非反常积分的Romberg算法1、代码：function z=IntRom(f,a,b) %Romberg 算法e=1e-10;I{1}=linspace(a,b,2);%1等分I{2}=linspace(a,b,3);%2等分L=setdiff(I{2},I{1});%新得插值点h=b-a;T(1,1)=h/2*sum(subs(f,I{1}));T(2,1)=1/2*T(1,1)+h/2*sum(subs(f,L));T(2,2)=4/3*T(2,1)-1/3*T(1,1);k=2;while abs(T(k,k)-T(k-1,k-1))>e %精度要求k=k+1;I{k}=linspace(a,b,2^(k-1)+1);L=setdiff(I{k},I{k-1});%集合差运算，新得插值点h=h/2;T(k,1)=1/2*T(k-1,1)+h/2*sum(subs(f,L));%梯形for i=2:kT(k,i)=(4^(i-1)/(4^(i-1)-1))*T(k,i-1)-(1/(4^(i-1)-1))*T(k-1,i-1);%加速endEndz=T(k,k);2、示例：>> syms x>> f=x^4;>> a=-4;b=4;>> IntRom(f,a,b)T =1.0e+003 *2.04800000000000 0 0 01.02400000000000 0.68266666666667 0 00.57600000000000 0.42666666666667 0.40960000000000 00.45200000000000 0.41066666666667 0.40960000000000 0.40960000000000ans =4.096000000000000e+002>> vpa((int(f,a,b)-ans),3)ans =0.>> f=exp(x);>> a=0;b=1;>> IntRom(f,a,b)T =Columns 1 through 41.85914091422952 0 0 01.75393109246483 1.71886115187659 0 01.72722190455752 1.71831884192175 1.71828268792476 01.72051859216430 1.71828415469990 1.71828184221844 1.718281828794531.71884112857999 1.71828197405189 1.71828182867536 1.718281828460391.71842166031633 1.71828183756177 1.71828182846243 1.71828182845905Columns 5 through 60 00 00 00 01.71828182845908 01.71828182845905 1.71828182845905ans =1.71828182845905>> vpa((int(f,a,b)-ans),3)ans =0.四、数值微分：5点法1、代码：function z=VDiff(f,x0)%5点法求导数值e=1e-15;h=0.01;for i=0:4x(i+1)=x0+i*h;endy=subs(f,x);m(1)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));%5点导数公式h=h/2;for i=0:4x(i+1)=x0+i*h;endy=subs(f,x);m(2)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));h=h/2;for i=-0:4x(i+1)=x0+i*h;endy=subs(f,x);m(3)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));k=3;while abs(m(k)-m(k-1))<abs(m(k-1)-m(k-2)) & abs(m(k)-m(k-1))>e & (h/10)>0%控制收敛条件及精度要求及h非0h=h/2;k=k+1;for i=0:4x(i+1)=x0+i*h;endy=subs(f,x);m(k)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));ende=abs(m(k)-m(k-1));z=[m(k);e];2、示例：>> syms x>> f=exp(x);>> x0=2;>> VDiff(f,x0)ans =7.389056098949710.00000000002558五、常微分方程：四阶龙格库塔及Adams加速法1、代码：四阶龙格库塔function z=RGFour(f,y0,a,b)%4阶龙格库塔，f为函数句柄h=0.01;X=a:h:b;Y(1)=y0;n=size(X);n=n(2);for i=1:n-1K1=f([X(i) Y(i)]);K2=f([X(i)+h/2,Y(i)+h/2*K1]);K3=f([X(i)+h/2,Y(i)+h/2*K2]);K4=f([X(i)+h,Y(i)+h*K3]);Y(i+1)=Y(i)+h/6*(K1+2*K2+2*K3+K4);endz=Y;plot(X,Y);2、示例：>> f=@(x)sin(x(1));>> y0=0;a=0;b=2*pi;>> figure(1);>> RGFour(f,y0,a,b);3、代码：Adams加速法function z=myAdams(f,y0,a,b)h=0.01;p(4)=1;c(4)=1;X=a:h:b;Y(1)=y0;n=size(X);n=n(2);for i=1:3K1=f([X(i) Y(i)]);K2=f([X(i)+h/2,Y(i)+h/2*K1]);K3=f([X(i)+h/2,Y(i)+h/2*K2]);K4=f([X(i)+h,Y(i)+h*K3]);Y(i+1)=Y(i)+h/6*(K1+2*K2+2*K3+K4);endfor i=4:n-1p(i+1)=Y(i)+h/24*(55*f([X(i),Y(i)])-59*f([X(i-1),Y(i-1)])+37*f([X(i-2 ),Y(i-2)])-9*f([X(i-3),Y(i-3)]));m(i+1)=p(i+1)+251/720*(c(i)-p(i));m(i+1)=f([X(i+1),m(i+1)]);c(i+1)=Y(i)+h/24*(9*f([X(i+1),m(i+1)])+19*f([X(i),Y(i)])-5*f([X(i-1), Y(i-1)])+f([X(i-2),Y(i-2)]));Y(i+1)=c(i+1)-19/720*(c(i+1)-p(i+1));endz=Y;plot(X,Y);4、示例：>> f=@(x)exp(x(1));>> myAdams(f,0,0,2*pi);六、方程求根：Aitken 迭代1、代码：function z=myAitken(f,x0);%Aitken 迭代求方程的根e=1e-15;xu1=x0;xv1=subs(f,xu1);xv2=subs(f,xv1);if xv2-2*xv1+xu1==0%防止除0；xu2=xv2;elsexu2=xv2-(xv2-xv1)^2/(xv2-2*xv1+xu1);endwhile abs(xu2-xu1)>e%精度控制xu1=xu2;xv1=subs(f,xu1);xv2=subs(f,xv1);if xv2-2*xv1+xu1==0%防止除0；xu2=xv2;elsexu2=xv2-(xv2-xv1)^2/(xv2-2*xv1+xu1);%Aitken加速公式endendz=xu2;2、示例：>> syms x>> f=cos(x/2)+x;>> x0=3;>> myAitken(f,x0)ans =3.14159265358979>> f=x^2-2+x;>> x0=1;>> myAitken(f,x0)ans =1.41421356237309七、线性方程组直接法：三角分解1、代码：function z=myGuess(A,b);%线性方程组三角分解求根n=size(A);if n~=rank(A)z=['矩阵A线性相关，不符合要求'];return;endn=n(2);L=eye(n);for i=1:n-1for j=i+1:nL(j,i)=A(j,i)/A(i,i);A(j,:)=A(j,:)-L(j,i)*A(i,:);endendU=A;for i=1:n %解Ly=b，得ys=0;for j=1:i-1s=s+y(j)*L(i,j);endy(i)=(b(i)-s)/L(i,i);endfor i=n:-1:1 %解Ux=y，得xs=0;for j=i+1:ns=s+x(j)*U(i,j);endx(i)=(y(i)-s)/U(i,i);endLUz=x';2、示例：>> A=[1 2 3;2 1 5;11 17 21];>> b=[1 3 5]';>> myGuess(A,b)L =1.00000000000000 0 02.00000000000000 1.00000000000000 011.00000000000000 1.66666666666667 1.00000000000000U =1.000000000000002.000000000000003.000000000000000 -3.00000000000000 -1.000000000000000 0 -10.33333333333333ans =-0.06451612903226-0.580645161290320.74193548387097>> t=A\bt =-0.06451612903226-0.580645161290320.74193548387097八、线性方程组迭代法：Jacobi法及G-S法1、代码：Jacobi法function z=myJacobi(A,b)n=size(A);n=n(2);x{1}=zeros(n,1);%初始值e=1e-10;D=diag(diag(A));L=D-tril(A);U=D-triu(A);B=D\(L+U);f=D\b;Q=B'*B;[w,d]=eig(Q);p=max(abs(diag(d))) ; %谱半径if p>=1z=['迭代发散'];return;endx{2}=B*x{1}+f;k=2;while norm(x{k}-x{k-1})>ek=k+1;x{k}=B*x{k-1}+f;endz=x{k};2、示例：>> A=[8 -3 2;4 11 -1;6 3 12];>> b=[20 33 36]';>> myJacobi(A,b)ans =3.000000000013402.000000000012610.999999999992373、代码：G-S法function z=myGS(A,b)n=size(A);n=n(2);x{1}=zeros(n,1);e=1e-10;D=diag(diag(A));L=D-tril(A);U=D-triu(A);B=(D-L)\U;f=(D-L)\b;Q=B'*B;[w,d]=eig(Q);p=max(abs(diag(d))) ; %谱半径if p>=1z=['迭代发散'];return;endx{2}=B*x{1}+f;k=2;while norm(x{k}-x{k-1})>ek=k+1;x{k}=B*x{k-1}+f;endz=x{k};4、示例：>> A=[8 -3 2;4 11 -1;6 3 12];>> b=[20 33 36]';>> myGS(A,b)ans =3.000000000001351.999999999999160.99999999999954九、矩阵的特征值及特征向量：幂法1、代码：function z=myChar(A);%幂法求主特征值及对应的特征向量e=1e-10;n=size(A);n=n(2);v1=ones(n,1);u1=v1;v2=A*u1;a=min(v2);b=max(v2);if abs(a)>abs(b)c=a;elsec=b;endu2=v2/c;%规范化while norm(u2-u1)>eu1=u2;v2=A*u1;a=min(v2);b=max(v2);if abs(a)>abs(b)c=a;elsec=b;endu2=v2/c;%规范化endz{1}=c;z{2}=v2;2、示例：>> A=[8 -3 2;4 11 -1;6 3 12];>> u=myChar(A)u =[14.00000000046956] [3x1 double]>> u{1}ans =14.00000000046956 >> u{2}ans =4.20000000191478 0.93333333198946 14.00000000046956。
matlab编程实例100例(精编文档).doc
【最新整理，下载后即可编辑】1-32是：图形应用篇33-66是：界面设计篇67-84是：图形处理篇85-100是：数值分析篇实例1：三角函数曲线（1）function shili01h0=figure('toolbar','none',...'position',[198****0300],...'name','实例01');h1=axes('parent',h0,...'visible','off');x=-pi:0.05:pi;y=sin(x);plot(x,y);xlabel('自变量X');ylabel('函数值Y');title('SIN( )函数曲线');grid on实例2：三角函数曲线（2）function shili02h0=figure('toolbar','none',...'position',[200 150 450 350],...'name','实例02');x=-pi:0.05:pi;y=sin(x)+cos(x);plot(x,y,'-*r','linewidth',1);grid onxlabel('自变量X');ylabel('函数值Y');title('三角函数');实例3：图形的叠加function shili03h0=figure('toolbar','none',...'position',[200 150 450 350],...'name','实例03');x=-pi:0.05:pi;y1=sin(x);y2=cos(x);plot(x,y1,...'-*r',...x,y2,...'--og');grid onxlabel('自变量X');ylabel('函数值Y');title('三角函数');实例4：双y轴图形的绘制function shili04h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例04');x=0:900;a=1000;b=0.005;y1=2*x;y2=cos(b*x);[haxes,hline1,hline2]=plotyy(x,y1,x,y2,'semilogy','plot'); axes(haxes(1))ylabel('semilog plot');axes(haxes(2))ylabel('linear plot');实例5：单个轴窗口显示多个图形function shili05h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例05');t=0:pi/10:2*pi;[x,y]=meshgrid(t);subplot(2,2,1)plot(sin(t),cos(t))axis equalsubplot(2,2,2)z=sin(x)-cos(y);plot(t,z)axis([0 2*pi -2 2])subplot(2,2,3)h=sin(x)+cos(y);plot(t,h)axis([0 2*pi -2 2])subplot(2,2,4)g=(sin(x).^2)-(cos(y).^2);plot(t,g)axis([0 2*pi -1 1])实例6：图形标注function shili06h0=figure('toolbar','none',...'position',[200 150 450 400],...'name','实例06');t=0:pi/10:2*pi;h=plot(t,sin(t));xlabel('t=0到2\pi','fontsize',16);ylabel('sin(t)','fontsize',16);title('\it{从0to2\pi 的正弦曲线}','fontsize',16) x=get(h,'xdata');y=get(h,'ydata');imin=find(min(y)==y);imax=find(max(y)==y);text(x(imin),y(imin),...['\leftarrow最小值=',num2str(y(imin))],...'fontsize',16)text(x(imax),y(imax),...['\leftarrow最大值=',num2str(y(imax))],...'fontsize',16)实例7：条形图形function shili07h0=figure('toolbar','none',...'position',[200 150 450 350],...'name','实例07');tiao1=[562 548 224 545 41 445 745 512];tiao2=[47 48 57 58 54 52 65 48];t=0:7;bar(t,tiao1)xlabel('X轴');ylabel('TIAO1值');h1=gca;h2=axes('position',get(h1,'position'));plot(t,tiao2,'linewidth',3)set(h2,'yaxislocation','right','color','none','xticklabel',[]) 实例8：区域图形function shili08h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例08');x=91:95;profits1=[88 75 84 93 77];profits2=[51 64 54 56 68];profits3=[42 54 34 25 24];profits4=[26 38 18 15 4];area(x,profits1,'facecolor',[0.5 0.9 0.6],...'edgecolor','b',...'linewidth',3)hold onarea(x,profits2,'facecolor',[0.9 0.85 0.7],...'edgecolor','y',...'linewidth',3)hold onarea(x,profits3,'facecolor',[0.3 0.6 0.7],...'edgecolor','r',...'linewidth',3)hold onarea(x,profits4,'facecolor',[0.6 0.5 0.9],...'edgecolor','m',...'linewidth',3)hold offset(gca,'xtick',[91:95])set(gca,'layer','top')gtext('\leftarrow第一季度销量') gtext('\leftarrow第二季度销量') gtext('\leftarrow第三季度销量') gtext('\leftarrow第四季度销量') xlabel('年','fontsize',16);ylabel('销售量','fontsize',16);实例9：饼图的绘制function shili09h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例09');t=[54 21 35;68 54 35;45 25 12;48 68 45;68 54 69];x=sum(t);h=pie(x);textobjs=findobj(h,'type','text');str1=get(textobjs,{'string'});val1=get(textobjs,{'extent'});oldext=cat(1,val1{:});names={'商品一：';'商品二：';'商品三：'};str2=strcat(names,str1);set(textobjs,{'string'},str2)val2=get(textobjs,{'extent'});newext=cat(1,val2{:});offset=sign(oldext(:,1)).*(newext(:,3)-oldext(:,3))/2; pos=get(textobjs,{'position'});textpos=cat(1,pos{:});textpos(:,1)=textpos(:,1)+offset;set(textobjs,{'position'},num2cell(textpos,[3,2]))实例10：阶梯图function shili10h0=figure('toolbar','none',...'position',[200 150 450 400],...'name','实例10');a=0.01;b=0.5;t=0:10;f=exp(-a*t).*sin(b*t);stairs(t,f)hold onplot(t,f,':*')hold offglabel='函数e^{-(\alpha*t)}sin\beta*t的阶梯图'; gtext(glabel,'fontsize',16)xlabel('t=0:10','fontsize',16)axis([0 10 -1.2 1.2])实例11：枝干图function shili11h0=figure('toolbar','none',...'position',[200 150 450 350],...'name','实例11');x=0:pi/20:2*pi;y1=sin(x);y2=cos(x);h1=stem(x,y1+y2);hold onh2=plot(x,y1,'^r',x,y2,'*g');hold offh3=[h1(1);h2];legend(h3,'y1+y2','y1=sin(x)','y2=cos(x)') xlabel('自变量X');ylabel('函数值Y');title('正弦函数与余弦函数的线性组合'); 实例12：罗盘图function shili12h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例12');winddirection=[54 24 65 84256 12 235 62125 324 34 254];windpower=[2 5 5 36 8 12 76 14 10 8];rdirection=winddirection*pi/180;[x,y]=pol2cart(rdirection,windpower); compass(x,y);desc={'风向和风力','北京气象台','10月1日0：00到','10月1日12：00'};gtext(desc)实例13：轮廓图function shili13h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例13');[th,r]=meshgrid((0:10:360)*pi/180,0:0.05:1); [x,y]=pol2cart(th,r);z=x+i*y;f=(z.^4-1).^(0.25);contour(x,y,abs(f),20)axis equalxlabel('实部','fontsize',16);ylabel('虚部','fontsize',16);h=polar([0 2*pi],[0 1]);delete(h)hold oncontour(x,y,abs(f),20)实例14：交互式图形function shili14h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例14');axis([0 10 0 10]);hold onx=[];y=[];n=0;disp('单击鼠标左键点取需要的点'); disp('单击鼠标右键点取最后一个点'); but=1;while but==1[xi,yi,but]=ginput(1);plot(xi,yi,'bo')n=n+1;disp('单击鼠标左键点取下一个点');x(n,1)=xi;y(n,1)=yi;endt=1:n;ts=1:0.1:n;xs=spline(t,x,ts);ys=spline(t,y,ts);plot(xs,ys,'r-');hold off实例14：交互式图形function shili14h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例14');axis([0 10 0 10]);hold onx=[];y=[];n=0;disp('单击鼠标左键点取需要的点'); disp('单击鼠标右键点取最后一个点'); but=1;while but==1[xi,yi,but]=ginput(1);plot(xi,yi,'bo')n=n+1;disp('单击鼠标左键点取下一个点');x(n,1)=xi;y(n,1)=yi;endt=1:n;ts=1:0.1:n;xs=spline(t,x,ts);ys=spline(t,y,ts);plot(xs,ys,'r-');hold off实例15：变换的傅立叶函数曲线function shili15h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例15');axis equalm=moviein(20,gcf);set(gca,'nextplot','replacechildren')h=uicontrol('style','slider','position',...[100 10 500 20],'min',1,'max',20)for j=1:20plot(fft(eye(j+16)))set(h,'value',j)m(:,j)=getframe(gcf);endclf;axes('position',[0 0 1 1]);movie(m,30)实例16：劳伦兹非线形方程的无序活动function shili15h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例15');axis equalm=moviein(20,gcf);set(gca,'nextplot','replacechildren')h=uicontrol('style','slider','position',...[100 10 500 20],'min',1,'max',20)for j=1:20plot(fft(eye(j+16)))set(h,'value',j)m(:,j)=getframe(gcf);endclf;axes('position',[0 0 1 1]);movie(m,30)实例17：填充图function shili17h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例17');t=(1:2:15)*pi/8;x=sin(t);y=cos(t);fill(x,y,'r')axis square offtext(0,0,'STOP',...'color',[1 1 1],...'fontsize',50,...'horizontalalignment','center') 例18：条形图和阶梯形图function shili18h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例18');subplot(2,2,1)x=-3:0.2:3;y=exp(-x.*x);bar(x,y)title('2-D Bar Chart')subplot(2,2,2)x=-3:0.2:3;y=exp(-x.*x);bar3(x,y,'r')title('3-D Bar Chart')subplot(2,2,3)x=-3:0.2:3;y=exp(-x.*x);stairs(x,y)title('Stair Chart')subplot(2,2,4)x=-3:0.2:3;y=exp(-x.*x);barh(x,y)title('Horizontal Bar Chart')实例19：三维曲线图function shili19h0=figure('toolbar','none',...'position',[200 150 450 400],...'name','实例19');subplot(2,1,1)x=linspace(0,2*pi);y1=sin(x);y2=cos(x);y3=sin(x)+cos(x);z1=zeros(size(x));z2=0.5*z1;z3=z1;plot3(x,y1,z1,x,y2,z2,x,y3,z3)grid onxlabel('X轴');ylabel('Y轴');zlabel('Z轴');title('Figure1:3-D Plot')subplot(2,1,2)x=linspace(0,2*pi);y1=sin(x);y2=cos(x);y3=sin(x)+cos(x);z1=zeros(size(x));z2=0.5*z1;z3=z1;plot3(x,z1,y1,x,z2,y2,x,z3,y3)grid onxlabel('X轴');ylabel('Y轴');zlabel('Z轴');title('Figure2:3-D Plot')实例20：图形的隐藏属性function shili20h0=figure('toolbar','none',...'position',[200 150 450 300],...'name','实例20');subplot(1,2,1)[x,y,z]=sphere(10);mesh(x,y,z)axis offtitle('Figure1:Opaque')hidden onsubplot(1,2,2)[x,y,z]=sphere(10);mesh(x,y,z)axis offtitle('Figure2:Transparent') hidden off实例21PEAKS函数曲线function shili21h0=figure('toolbar','none',...'position',[200 100 450 450],...'name','实例21');[x,y,z]=peaks(30);subplot(2,1,1)x=x(1,:);y=y(:,1);i=find(y>0.8&y<1.2);j=find(x>-0.6&x<0.5);z(i,j)=nan*z(i,j);surfc(x,y,z)xlabel('X轴');ylabel('Y轴');zlabel('Z轴');title('Figure1:surfc函数形成的曲面') subplot(2,1,2)x=x(1,:);y=y(:,1);i=find(y>0.8&y<1.2);j=find(x>-0.6&x<0.5);z(i,j)=nan*z(i,j);surfl(x,y,z)xlabel('X轴');ylabel('Y轴');zlabel('Z轴');title('Figure2:surfl函数形成的曲面') 实例22：片状图function shili22h0=figure('toolbar','none',...'position',[200 150 550 350],...'name','实例22');subplot(1,2,1)x=rand(1,20);y=rand(1,20);z=peaks(x,y*pi);t=delaunay(x,y);trimesh(t,x,y,z)hidden offtitle('Figure1:Triangular Surface Plot'); subplot(1,2,2)x=rand(1,20);y=rand(1,20);z=peaks(x,y*pi);t=delaunay(x,y);trisurf(t,x,y,z)title('Figure1:Triangular Surface Plot'); 实例23：视角的调整function shili23h0=figure('toolbar','none',...'position',[200 150 450 350],...'name','实例23');x=-5:0.5:5;[x,y]=meshgrid(x);r=sqrt(x.^2+y.^2)+eps;z=sin(r)./r;subplot(2,2,1)surf(x,y,z)xlabel('X-axis')ylabel('Y-axis')zlabel('Z-axis')title('Figure1')view(-37.5,30)subplot(2,2,2)surf(x,y,z)xlabel('X-axis')ylabel('Y-axis')zlabel('Z-axis')title('Figure2')view(-37.5+90,30) subplot(2,2,3)surf(x,y,z)xlabel('X-axis')ylabel('Y-axis')zlabel('Z-axis')title('Figure3')view(-37.5,60)subplot(2,2,4)surf(x,y,z)xlabel('X-axis')ylabel('Y-axis')zlabel('Z-axis')title('Figure4')view(180,0)实例24：向量场的绘制function shili24h0=figure('toolbar','none',...'position',[200 150 450 350],...'name','实例24');subplot(2,2,1)z=peaks;ribbon(z)title('Figure1')subplot(2,2,2)[x,y,z]=peaks(15);[dx,dy]=gradient(z,0.5,0.5); contour(x,y,z,10)hold onquiver(x,y,dx,dy)hold offtitle('Figure2')subplot(2,2,3)[x,y,z]=peaks(15);[nx,ny,nz]=surfnorm(x,y,z);surf(x,y,z)hold onquiver3(x,y,z,nx,ny,nz)hold offtitle('Figure3')subplot(2,2,4)x=rand(3,5);y=rand(3,5);z=rand(3,5);c=rand(3,5);fill3(x,y,z,c)grid ontitle('Figure4')实例25：灯光定位function shili25h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例25');vert=[1 1 1;1 2 1;2 2 1;2 1 1;1 1 2;12 2;2 2 2;2 1 2];fac=[1 2 3 4;2 6 7 3;4 3 7 8;15 8 4;1 2 6 5;5 6 7 8];grid offsphere(36)h=findobj('type','surface');set(h,'facelighting','phong',...'facecolor',...'interp',...'edgecolor',[0.4 0.4 0.4],...'backfacelighting',...'lit')hold onpatch('faces',fac,'vertices',vert,...'facecolor','y');light('position',[1 3 2]);light('position',[-3 -1 3]);material shinyaxis vis3d offhold off实例26：柱状图function shili26h0=figure('toolbar','none',...'position',[200 50 450 450],...'name','实例26'); subplot(2,1,1)x=[5 2 18 7 39 8 65 5 54 3 2];bar(x)xlabel('X轴');ylabel('Y轴');title('第一子图');subplot(2,1,2)y=[5 2 18 7 39 8 65 5 54 3 2];barh(y)xlabel('X轴');ylabel('Y轴');title('第二子图');实例27：设置照明方式function shili27h0=figure('toolbar','none',...'position',[200 150 450 350],...'name','实例27');subplot(2,2,1)sphereshading flatcamlight leftcamlight rightlighting flatcolorbaraxis offtitle('Figure1')subplot(2,2,2)sphereshading flatcamlight leftcamlight rightlighting gouraudcolorbaraxis offtitle('Figure2')subplot(2,2,3)sphereshading interpcamlight rightcamlight leftlighting phongaxis offtitle('Figure3')subplot(2,2,4)sphereshading flatcamlight leftcamlight rightlighting nonecolorbaraxis offtitle('Figure4')实例28：羽状图function shili28h0=figure('toolbar','none',...'position',[200 150 450 350],...'name','实例28');subplot(2,1,1)alpha=90:-10:0;r=ones(size(alpha));m=alpha*pi/180;n=r*10;[u,v]=pol2cart(m,n);feather(u,v)title('羽状图')axis([0 20 0 10])subplot(2,1,2)t=0:0.5:10;y=exp(-x*t);feather(y)title('复数矩阵的羽状图')实例29：立体透视（1）function shili29h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例29');[x,y,z]=meshgrid(-2:0.1:2,...-2:0.1:2,...-2:0.1:2);v=x.*exp(-x.^2-y.^2-z.^2);grid onfor i=-2:0.5:2;h1=surf(linspace(-2,2,20),...linspace(-2,2,20),...zeros(20)+i);rotate(h1,[1 -1 1],30)dx=get(h1,'xdata');dy=get(h1,'ydata');dz=get(h1,'zdata');delete(h1)slice(x,y,z,v,[-2 2],2,-2)hold onslice(x,y,z,v,dx,dy,dz)hold offaxis tightview(-5,10)drawnowend实例30：立体透视（2）function shili30h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例30');[x,y,z]=meshgrid(-2:0.1:2,...-2:0.1:2,...-2:0.1:2);v=x.*exp(-x.^2-y.^2-z.^2); [dx,dy,dz]=cylinder;slice(x,y,z,v,[-2 2],2,-2)for i=-2:0.2:2h=surface(dx+i,dy,dz);rotate(h,[1 0 0],90)xp=get(h,'xdata');yp=get(h,'ydata');zp=get(h,'zdata');delete(h)hold onhs=slice(x,y,z,v,xp,yp,zp);axis tightxlim([-3 3])view(-10,35)drawnowdelete(hs)hold offend实例31：表面图形function shili31h0=figure('toolbar','none',...'position',[200 150 550 250],...'name','实例31');subplot(1,2,1)x=rand(100,1)*16-8;y=rand(100,1)*16-8;r=sqrt(x.^2+y.^2)+eps;z=sin(r)./r;xlin=linspace(min(x),max(x),33); ylin=linspace(min(y),max(y),33); [X,Y]=meshgrid(xlin,ylin);Z=griddata(x,y,z,X,Y,'cubic'); mesh(X,Y,Z)axis tighthold onplot3(x,y,z,'.','Markersize',20) subplot(1,2,2)k=5;n=2^k-1;theta=pi*(-n:2:n)/n;phi=(pi/2)*(-n:2:n)'/n;X=cos(phi)*cos(theta);Y=cos(phi)*sin(theta);Z=sin(phi)*ones(size(theta)); colormap([0 0 0;1 1 1])C=hadamard(2^k);surf(X,Y,Z,C)axis square实例32：沿曲线移动的小球h0=figure('toolbar','none',...'position',[198****8468],...'name','实例32');h1=axes('parent',h0,...'position',[0.15 0.45 0.7 0.5],...'visible','on');t=0:pi/24:4*pi;y=sin(t);plot(t,y,'b')n=length(t);h=line('color',[0 0.5 0.5],...'linestyle','.',...'markersize',25,...'erasemode','xor');k1=uicontrol('parent',h0,...'style','pushbutton',...'position',[80 100 50 30],...'string','开始',...'callback',[...'i=1;',...'k=1;,',...'m=0;,',...'while 1,',...'if k==0,',...'break,',...'end,',...'if k~=0,',...'set(h,''xdata'',t(i),''ydata'',y(i)),',...'drawnow;,',...'i=i+1;,',...'if i>n,',...'m=m+1;,',...'i=1;,',...'end,',...'end,',...'end']);k2=uicontrol('parent',h0,...'style','pushbutton',...'position',[180 100 50 30],...'string','停止',...'callback',[...'k=0;,',...'set(e1,''string'',m),',...'p=get(h,''xdata'');,',...'q=get(h,''ydata'');,',...'set(e2,''string'',p);,',...'set(e3,''string'',q)']); k3=uicontrol('parent',h0,...'style','pushbutton',...'position',[280 100 50 30],...'string','关闭',...'callback','close');e1=uicontrol('parent',h0,...'style','edit',...'position',[60 30 60 20]);t1=uicontrol('parent',h0,...'style','text',...'string','循环次数',...'position',[60 50 60 20]);e2=uicontrol('parent',h0,...'style','edit',...'position',[180 30 50 20]);t2=uicontrol('parent',h0,...'style','text',...'string','终点的X坐标值',...'position',[155 50 100 20]);e3=uicontrol('parent',h0,...'style','edit',...'position',[300 30 50 20]);t3=uicontrol('parent',h0,...'style','text',...'string','终点的Y坐标值',...'position',[275 50 100 20]);实例33：曲线转换按钮h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例33');x=0:0.5:2*pi;y=sin(x);h=plot(x,y);grid onhuidiao=[...'if i==1,',...'i=0;,',...'y=cos(x);,',...'delete(h),',...'set(hm,''string'',''正弦函数''),',...'h=plot(x,y);,',...'grid on,',...'else if i==0,',...'i=1;,',...'y=sin(x);,',...'set(hm,''string'',''余弦函数''),',...'delete(h),',...'h=plot(x,y);,',...'grid on,',...'end,',...'end'];hm=uicontrol(gcf,'style','pushbutton',...'string','余弦函数',...'callback',huidiao);i=1;set(hm,'position',[250 20 60 20]);set(gca,'position',[0.2 0.2 0.6 0.6])title('按钮的使用')hold on实例34：栅格控制按钮h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例34');x=0:0.5:2*pi;y=sin(x);plot(x,y)huidiao1=[...'set(h_toggle2,''value'',0),',...'grid on,',...];huidiao2=[...'set(h_toggle1,''value'',0),',...'grid off,',...];h_toggle1=uicontrol(gcf,'style','togglebutton',...'string','grid on',...'value',0,...'position',[20 45 50 20],...'callback',huidiao1);h_toggle2=uicontrol(gcf,'style','togglebutton',...'string','grid off',...'value',0,...'position',[20 20 50 20],...'callback',huidiao2);set(gca,'position',[0.2 0.2 0.6 0.6])title('开关按钮的使用')实例35：编辑框的使用h0=figure('toolbar','none',...'position',[200 150 350 250],...'name','实例35');f='Please input the letter';huidiao1=[...'g=upper(f);,',...'set(h2_edit,''string'',g),',...];huidiao2=[...'g=lower(f);,',...'set(h2_edit,''string'',g),',...];h1_edit=uicontrol(gcf,'style','edit',...'position',[100 200 100 50],...'HorizontalAlignment','left',...'string','Please input the letter',...'callback','f=get(h1_edit,''string'');',...'background','w',...'max',5,...'min',1);h2_edit=uicontrol(gcf,'style','edit',...'HorizontalAlignment','left',...'position',[100 100 100 50],...'background','w',...'max',5,...'min',1);h1_button=uicontrol(gcf,'style','pushbutton',...'string','小写变大写',...'position',[100 45 100 20],...'callback',huidiao1);h2_button=uicontrol(gcf,'style','pushbutton',...'string','大写变小写',...'position',[100 20 100 20],...'callback',huidiao2);实例36：弹出式菜单h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例36');x=0:0.5:2*pi;y=sin(x);h=plot(x,y);grid onhm=uicontrol(gcf,'style','popupmenu',...'string',...'sin(x)|cos(x)|sin(x)+cos(x)|exp(-sin(x))',...'position',[250 20 50 20]);set(hm,'value',1)huidiao=[...'v=get(hm,''value'');,',...'switch v,',...'case 1,',...'delete(h),',...'y=sin(x);,',...'h=plot(x,y);,',...'grid on,',...'case 2,',...'delete(h),',...'y=cos(x);,',...'h=plot(x,y);,',...'grid on,',...'case 3,',...'delete(h),',...'y=sin(x)+cos(x);,',...'h=plot(x,y);,',...'grid on,',...'case 4,',...'delete(h),',...'y=exp(-sin(x));,',...'h=plot(x,y);,',...'grid on,',...'end'];set(hm,'callback',huidiao)set(gca,'position',[0.2 0.2 0.6 0.6]) title('弹出式菜单的使用')实例37：滑标的使用h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例37');[x,y]=meshgrid(-8:0.5:8);r=sqrt(x.^2+y.^2)+eps;z=sin(r)./r;h0=mesh(x,y,z);h1=axes('position',...[0.2 0.2 0.5 0.5],...'visible','off');htext=uicontrol(gcf,...'units','points',...'position',[20 30 45 15],...'string','brightness',...'style','text');hslider=uicontrol(gcf,...'units','points',...'position',[10 10 300 15],...'min',-1,...'max',1,...'style','slider',...'callback',...'brighten(get(hslider,''value''))'); 实例38：多选菜单h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例38');[x,y]=meshgrid(-8:0.5:8);r=sqrt(x.^2+y.^2)+eps;z=sin(r)./r;h0=mesh(x,y,z);hlist=uicontrol(gcf,'style','listbox',...'string','default|spring|summer|autumn|winter',...'max',5,...'min',1,...'position',[20 20 80 100],...'callback',[...'k=get(hlist,''value'');,',...'switch k,',...'case 1,',...'colormap default,',...'case 2,',...'colormap spring,',...'case 3,',...'colormap summer,',...'case 4,',...'colormap autumn,',...'case 5,',...'colormap winter,',...'end']);实例39：菜单控制的使用h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例39');x=0:0.5:2*pi;y=cos(x);h=plot(x,y);grid onset(gcf,'toolbar','none')hm=uimenu('label','example');huidiao1=[...'set(hm_gridon,''checked'',''on''),',...'set(hm_gridoff,''checked'',''off''),',...'grid on'];huidiao2=[...'set(hm_gridoff,''checked'',''on''),',...'set(hm_gridon,''checked'',''off''),',...'grid off'];hm_gridon=uimenu(hm,'label','grid on',...'checked','on',...'callback',huidiao1);hm_gridoff=uimenu(hm,'label','grid off',...'checked','off',...'callback',huidiao2);实例40：UIMENU菜单的应用h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例40');h1=uimenu(gcf,'label','函数');h11=uimenu(h1,'label','轮廓图',...'callback',[...'set(h31,''checked'',''on''),',...'set(h32,''checked'',''off''),',...'[x,y,z]=peaks;,',...'contour3(x,y,z,30)']);h12=uimenu(h1,'label','高斯分布',...'callback',[...'set(h31,''checked'',''on''),',...'set(h32,''checked'',''off''),',...'mesh(peaks);,',...'axis tight']);。
MATLAB_智能算法30个案例分析
MATLAB_智能算法30个案例分析MATLAB是一种强大的数值计算和编程工具，教育和科研领域中广泛应用于数据分析、机器学习和智能算法的研究。
在本文中，我们将介绍30个MATLAB智能算法的案例分析，并探讨其用途和优势。
分析的案例包括分类、回归、聚类、神经网络和遗传算法等不同类型的智能算法。
1. K均值聚类：利用MATLAB中的kmeans函数对一组数据进行聚类分析，得到不同的簇。
2. 随机森林：利用MATLAB中的TreeBagger函数构建一个随机森林分类器，并通过测试数据进行分类预测。
3. 人工神经网络：使用MATLAB中的feedforwardnet函数构建一个人工神经网络，并通过训练集进行预测。
4. 遗传算法：利用MATLAB中的ga函数对一个优化问题进行求解，找到最优解。
5. 支持向量机：使用MATLAB中的svmtrain和svmclassify函数构建一个支持向量机分类器，并进行分类预测。
6. 极限学习机：使用MATLAB中的elmtrain和elmpredict函数构建一个极限学习机分类器，并进行分类预测。
7. 逻辑回归：使用MATLAB中的mnrfit和mnrval函数构建一个逻辑回归模型，并进行预测。
8. 隐马尔可夫模型：使用MATLAB中的hmmtrain和hmmdecode函数构建一个隐马尔可夫模型，对一系列观测数据进行预测。
9. 神经进化算法：利用MATLAB中的ne_train函数构建一个基于神经进化算法的神经网络分类器，并进行分类预测。
10. 朴素贝叶斯分类器：使用MATLAB中的NaiveBayes对象构建一个朴素贝叶斯分类器，并进行分类预测。
11. 高斯过程回归：使用MATLAB中的fitrgp函数构建一个高斯过程回归模型，并进行回归预测。
12. 最小二乘支持向量机：使用MATLAB中的fitcsvm函数构建一个最小二乘支持向量机分类器，并进行分类预测。
13. 遗传网络：利用MATLAB中的ngenetic函数构建一个基于遗传算法和人工神经网络的分类器，并进行分类预测。
数值分析matlab
plot(90,389,'*r')
hold on
grid
A=[30,50,70,90];
B=[133,221,305,Βιβλιοθήκη 89];norm(A,2)
norm(B,2)
程序运行如图：
x =
40 60 80
y =
177 263 352
p2 =
0.0037 3.9250 14.0000
二.解题方法描述、程序代码、结果展示
x=[40,60,80]
y=[177,263,352]
plot(x,y,'*r')
hold on
grid
p2=polyfit(x,y,2)
py2=polyval(p2,x)
plot(x,py2,'b')
a=polyfit(x,y,length(x)-1);
poly2sym(a)
一.问题
已知E和C的关系如下表，E为自变量，C为函数
n
0
1
2
3
4
5
6
En
30
40
50
60
70
80
90
Cn
133
177
221
263
305
352
389
求：以n=1、3、5，三点为测试点进行2阶多项式插值和线性拟合，在一幅图中分别画出7个离散的点和两条曲线（插值和拟合各一条），并以n=0、2、4、6为检验点进行二范数误差分析。
plot(x,y);
图像：
求得多项式为0.0037x^2+3.9250x+14.0000.
用matlab解析实际案例
用MATLAB解析实际案例简介MATLAB是一种面向科学计算、可视化和数学建模的高级技术计算软件。
它被广泛应用于研究、工程和行业中，如信号处理、图像处理、机器学习、控制系统设计等。
本文将使用MATLAB来解析一个实际案例，并说明如何使用MATLAB进行数据分析和可视化。
案例背景某国外航空公司的一架飞机在飞行途中突然失去了控制，最终在一片森林地带坠毁，导致机上所有人员全部遇难。
该事件引起了广泛关注和调查。
在调查过程中，分析了飞机事故中的数据记录器，也被称为“黑匣子”，以便更好地了解整个事件。
数据记录器记录了飞机在飞行途中的多种参数，包括高度、速度、航向、油门设置等。
本文将使用MATLAB对这些数据进行分析和可视化，以更好地理解飞机失控的原因。
数据处理首先，我们需要将数据导入MATLAB。
这个过程可以通过多种方式实现，例如读取CSV文件、导入Excel文件或使用数据库连接器。
在本例中，原始数据以CSV 文件的形式提供。
为了导入数据，我们在MATLAB中使用readtable函数：data = readtable('data.csv');这会将数据读入一个MATLAB表中，其中每一列代表一个不同的参数。
我们可以使用summary函数查看数据摘要：summary(data)此命令将显示数据的摘要信息、变量名称和其他统计数据。
接下来，我们可以使用MATLAB的plot函数来绘制数据。
例如，要显示速度随时间变化的情况：plot(data.Time, data.Speed);这将在MATLAB中绘制速度与时间的关系图。
在绘制图形之前，我们需要对数据进行一些预处理。
在本例中，我们需要使用曲线拟合技术来找到可能导致飞机事故的异常点。
MATLAB提供了多种曲线拟合技术，例如多项式拟合、趋势线拟合等。
在本例中，我们将使用Loess曲线拟合技术来找到异常点。
首先，我们需要使用MATLAB的loess函数来计算数据的平滑曲线：smooth_data = loess(data.Time, data.Speed);然后，我们可以使用MATLAB的detrend函数来计算数据的趋势：detrended_data = detrend(data.Speed - smooth_data);这将返回数据的趋势，并将其从原始数据中减去，以获得数据的异常部分。
基于MATLAB的数值分析(4)
（1）syms x y u fai t;
f=x/(1+u^2);g=cos(y+fai);fg1=compose(f,g)
fg1 = cos(y+fai)/(1+u^2) （2）fg2=compose(f,g,u,fai,t) ） fg2 = x/(cos(y+t)^2+1)
三、置换操作
1、 子表达式置换操作 、 [RS,ssub]=subexpr(S,ssub) 运用符号变量ssub置换子表达式，重写 为RS. 置换子表达式， 运用符号变量 置换子表达式 重写S为 把复杂表达式中所含的多个相同子表达式用一个符号代替， 【例】把复杂表达式中所含的多个相同子表达式用一个符号代替， 使表达简洁。 使表达简洁。 clear all syms a b c d W [V,D]=eig([a b;c d])； ； [RVD,W]=subexpr([V;D],W) RVD = [ -(1/2*d-1/2*a-1/2*W)/c, -(1/2*d-1/2*a+1/2*W)/c] [ 1, 1] [ 1/2*d+1/2*a+1/2*W, 0] [ 0, 1/2*d+1/2*a-1/2*W] W = (d^2-2*a*d+a^2+4*b*c)^(1/2)
4.3 符号微积分
一、 符号序列的求和 s=symsum(f,v,a,b) : a,b缺省时，自变量区间 缺省时， 缺省时 为[0,v-1]。 。
【例】求
∑ [t
t −1 t =0
k3
]
1 ，∑ (2k − 1) 2 k =1
∞
(−1) k k
syms k t;f1=[t k^3];f2=[1/(2*k-1)^2,(-1)^k/k]; s1=simple(symsum(f1)) s2=simple(symsum(f2,1,inf)) s1 =[ 1/2*t*(t-1), k^3*t] s2 =[ 1/8*pi^2, -log(2)]
2-MATLAB在数值分析中的应用
分段插值
一维插值：可以分为最近插值、线性插值、三次样 条插值，分段三次Hermite插值。 y=interp1(x0,y0,x) y=interp1(x0,y0,x,’method’) method= nearst: 最近插值 linear: 线性插值（默认值） spline: 分段三次样条插值 pchip:分段三次Hermite插值
-9-
Polyinterp(symbol)
symx =sym(‘x’) P = polyinterp(x0,y0,symx) pretty(P) P = simplify(P) P= x^3-2*x-5 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)
0.5 0 1 2 3 4 5 6 7 8 9
-2-
插值和拟合
在平面上给定n个点(xk,yk),可以唯一确定一个最多n-1次 的多项式通过这些点，这个多项式叫插值多项式
P(xk ) = yk , k = 1,2,…,n Lagrange插值多项式
x xj P( x ) k 1 j k xk x j
2 3
sk ,0 yk sk ,2 P xk 2
sk ,1 dk sk ,3
hk 2P xk P xk 1
6 P xk 1 P xk 6hk
-21-
x = 1:6; y = [16 18 21 17 15 12];
Lagrange插值多项式基函数
P( x) c1 xn1 c2 xn2 ... cn1 x cn
n 1 x1 n 1 x2 x n 1 n n2 x1 n2 x2
数值分析matlab实验报告
数值分析matlab实验报告数值分析 Matlab 实验报告一、实验目的数值分析是研究各种数学问题数值解法的学科，Matlab 则是一款功能强大的科学计算软件。
本次实验旨在通过使用 Matlab 解决一系列数值分析问题，加深对数值分析方法的理解和应用能力，掌握数值计算中的误差分析、数值逼近、数值积分与数值微分等基本概念和方法，并培养运用计算机解决实际数学问题的能力。
二、实验内容（一）误差分析在数值计算中，误差是不可避免的。
通过对给定函数进行计算，分析截断误差和舍入误差的影响。
例如，计算函数＄f(x) ＝＼sin(x)＄在＄x ＝ 05$ 附近的值，比较不同精度下的结果差异。
（二）数值逼近1、多项式插值使用拉格朗日插值法和牛顿插值法对给定的数据点进行插值，得到拟合多项式，并分析其误差。
2、曲线拟合采用最小二乘法对给定的数据进行线性和非线性曲线拟合，如多项式曲线拟合和指数曲线拟合。
（三）数值积分1、牛顿柯特斯公式实现梯形公式、辛普森公式和柯特斯公式，计算给定函数在特定区间上的积分值，并分析误差。
2、高斯求积公式使用高斯勒让德求积公式计算积分，比较其精度与牛顿柯特斯公式的差异。
（四）数值微分利用差商公式计算函数的数值导数，分析步长对结果的影响，探讨如何选择合适的步长以提高精度。
三、实验步骤（一）误差分析1、定义函数｀compute_sin_error` 来计算不同精度下的正弦函数值和误差。
｀｀｀matlabfunction value, error ＝ compute_sin_error(x, precision)true_value ＝ sin(x)；computed_value ＝ vpa(sin(x)， precision)；error ＝ abs(true_value computed_value)；end｀｀｀2、在主程序中调用该函数，分别设置不同的精度进行计算和分析。
（二）数值逼近1、拉格朗日插值法｀｀｀matlabfunction L ＝ lagrange_interpolation(x, y, xi)n ＝ length(x)；L ＝ 0;for i ＝ 1:nli ＝ 1;for j ＝ 1:nif j ～＝ ili ＝ li （xi x(j)）／（x(i) x(j)）；endendL ＝ L ＋ y(i) li;endend｀｀｀2、牛顿插值法｀｀｀matlabfunction N ＝ newton_interpolation(x, y, xi)n ＝ length(x)；％计算差商表D ＝ zeros(n, n)；D(：， 1) ＝ y'；for j ＝ 2:nfor i ＝ j:nD(i, j) ＝（D(i, j 1) D(i 1, j 1)）／（x(i) x(i j ＋ 1)）；endend％计算插值结果N ＝ D(1, 1)；term ＝ 1;for i ＝ 2:nterm ＝ term （xi x(i 1)）；N ＝ N ＋ D(i, i) term;endend｀｀｀3、曲线拟合｀｀｀matlab％线性最小二乘拟合p ＝ polyfit(x, y, 1)；y_fit_linear ＝ polyval(p, x)；％多项式曲线拟合p ＝ polyfit(x, y, n)；％ n 为多项式的次数y_fit_poly ＝ polyval(p, x)；％指数曲线拟合p ＝ fit(x, y, ＇exp1'）；y_fit_exp ＝ p(x)；｀｀｀（三）数值积分1、梯形公式｀｀｀matlabfunction T ＝ trapezoidal_rule(f, a, b, n)h ＝（b a) ／ n;x ＝ a:h:b;y ＝ f(x)；T ＝ h （（y(1) ＋ y(end)）／ 2 ＋ sum(y(2:end 1)））；end｀｀｀2、辛普森公式｀｀｀matlabfunction S ＝ simpson_rule(f, a, b, n)if mod(n, 2) ～＝ 0error(＇n 必须为偶数'）；endh ＝（b a) ／ n;x ＝ a:h:b;y ＝ f(x)；S ＝ h ／ 3 （y(1) ＋ 4 sum(y(2:2:end 1)）＋ 2 sum(y(3:2:end 2)）＋ y(end)）；end｀｀｀3、柯特斯公式｀｀｀matlabfunction C ＝ cotes_rule(f, a, b, n)h ＝（b a) ／ n;x ＝ a:h:b;y ＝ f(x)；w ＝ 7, 32, 12, 32, 7 ／ 90;C ＝ h sum(w y)；end｀｀｀4、高斯勒让德求积公式｀｀｀matlabfunction G ＝ gauss_legendre_integration(f, a, b)x, w ＝ gauss_legendre(5)；％选择适当的节点数t ＝（b a) ／ 2 x ＋（a ＋ b) ／ 2;G ＝（b a) ／ 2 sum(w f(t)）；end｀｀｀（四）数值微分｀｀｀matlabfunction dydx ＝ numerical_derivative(f, x, h)dydx ＝（f(x ＋ h) f(x h)）／（2 h)；end｀｀｀四、实验结果与分析（一）误差分析通过不同精度的计算，发现随着精度的提高，误差逐渐减小，但计算时间也相应增加。
数值分析在生活中的应用举例及Matlab实现
一、最小二乘法，用MATLAB实现1. 数值实例下面给定的是乌鲁木齐最近1个月早晨7：00左右（新疆时间）的天气预报所得到的温度，按照数据找出任意次曲线拟合方程和它的图像。
下面用MATLAB编程对上述数据进行最小二乘拟合。
下面用MATLAB编程对上述数据进行最小二乘拟合2、程序代码x=[1:1:30];y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1];a1=polyfit(x,y,3) %三次多项式拟合%a2= polyfit(x,y,9) %九次多项式拟合%a3= polyfit(x,y,15) %十五次多项式拟合%b1=polyval(a1,x)b2=polyval(a2,x)b3=polyval(a3,x)r1= sum((y-b1).^2) %三次多项式误差平方和%r2= sum((y-b2).^2) %九次次多项式误差平方和%r3= sum((y-b3).^2) %十五次多项式误差平方和%plot(x,y,'*') %用*画出x,y图像%hold onplot(x,b1, 'r') %用红色线画出x,b1图像%hold onplot(x,b2, 'g') %用绿色线画出x,b2图像%hold onplot(x,b3, 'b:o') %用蓝色o线画出x,b3图像%3、数值结果不同次数多项式拟合误差平方和为：r1=67.6659r2=20.1060r3=3.7952r1、r2、r3分别表示三次、九次、十五次多项式误差平方和。
4、拟合曲线如下图二、 线性方程组的求解（ 高斯-塞德尔迭代算法 ）1、实例： 求解线性方程组（见书P233页）⎪⎪⎩⎪⎪⎨⎧=++=-+=+-3612363311420238321321321x x x x x x x x x 记A x=b, 其中⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=363320,,12361114238321b x A x x x任取初始值()()Tx0000=，进行迭代。
数值分析matlab数值试验
实验一：误差传播及算法稳定性实验1.21、试验程序：function charpt1_2% 误差传播及算法稳定性实验clc;clear all;promps={'请选择递推关系式，若选E1=1/e,En=1-nEn-1,请输入1，若选EN=0，En-1=(1-En)/n,请输入2：'};I=1;while Iresult=inputdlg(promps,'charpt1_2',1,{'1'});Nb=str2num(char(result));if ((Nb~=1)|(Nb~=2))I=0;endend%%%%%%%%%%%%%%%%%I=1;while Iresult=inputdlg('请输入递推步数 n>=1:','charpt1_2',1,{'10'});steps=str2num(char(result));if (steps>0)&(steps==fix(steps)) %% 如果steps大于0且为整数I=0;endend%%%%%%%%%%%%%%%%%result=inputdlg('请输入计算中所采用的有效数字位数n:','charpt1_2',1,{'5'});Sd=str2num(char(result));format long %% 设置显示精度result=zeros(1,steps); %% 存储计算结果err=result; %% 存储计算的绝对误差值func=result; %% 存储用quadl计算的近似值%%%%%%%%%%%%%%%%%%% 用quadl计算积分近似值for n=1:stepsfun=@(x) x.^n.*exp(x-1);func(n)= quadl(fun,0,1);end%%%%%%%%%%%%%%%%%%% 用自定义算法计算if(Nb==1)digits(Sd);result(1)=subs(vpa(1/exp(1)));for n=2:stepsresult(n)=subs(vpa(1-n*result(n-1)));enderr=abs(result-func);elseif(Nb==2)digits(Sd);result(steps)=0;for n=(steps-1):-1:1result(n)=subs(vpa((1-result(n+1))/(n+1)));enderr=abs(result-func);end%%%%%%%%%%%%%%%%%%% 输出结果数值及图像clf;disp('库函数计算值:');disp(sprintf('%e ',func));disp('递推值:');disp(sprintf('%e ',result));disp('误差值:');disp(sprintf('%e ',err));if(Nb==1)plot([1:steps],result,'-rs',[1:steps],func,':k*',[1:steps],err,'-.bo' );elseif(Nb==2)plot([steps:-1:1],result,'-rs',[steps:-1:1],func,':k*',[steps:-1:1],e rr,'-.bo');endxlabel('第n步');ylabel('计算值');legend('自定义算法结果','库函数计算结果','误差值');grid on2、试验结果：选择递推关系式1，递推步数为10，有效数字为5位，计算结果如下：库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-0011.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002递推值:3.678800e-001 2.642400e-001 2.072800e-001 1.708800e-001 1.456000e-0011.264000e-001 1.152000e-001 7.840000e-0022.944000e-001 -1.944000e+000误差值:5.588280e-007 1.117662e-006 3.352927e-006 1.341222e-0056.705713e-005 4.023702e-004 2.816427e-003 2.253226e-002 2.027877e-001 2.027877e+00012345678910第n 步计算值选择递推关系式2，递推步数为10，有效数字为5位，计算结果如下： 库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-001 1.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002 递推值:3.678800e-001 2.642400e-001 2.072800e-001 1.708900e-001 1.455300e-001 1.267900e-001 1.125000e-001 1.000000e-001 1.000000e-001 0.000000e+000 误差值:5.588280e-007 1.117662e-006 3.352927e-006 3.412224e-006 2.942873e-006 1.237016e-005 1.164270e-004 9.322618e-004 8.387707e-003 8.387707e-002第n 步计算值选择递推关系式1，递推步数为10，有效数字为6位，计算结果如下： 库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-001 1.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002 递推值:3.678790e-001 2.642420e-001 2.072740e-001 1.709040e-001 1.454800e-001 1.271200e-001 1.101600e-001 1.187200e-001 -6.848000e-002 1.684800e+000 误差值:4.411720e-007 8.823378e-007 2.647073e-006 1.058778e-0055.294287e-005 3.176298e-004 2.223573e-003 1.778774e-002 1.600923e-001 1.600923e+00012345678910第n 步计算值选择递推关系式2，递推步数为10，有效数字为6位，计算结果如下： 库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-001 1.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002 递推值:3.678800e-001 2.642410e-001 2.072770e-001 1.708930e-001 1.455360e-001 1.267860e-001 1.125000e-001 1.000000e-001 1.000000e-001 0.000000e+000 误差值:5.588280e-007 1.176622e-007 3.529274e-007 4.122239e-007 3.057127e-006 1.637016e-005 1.164270e-004 9.322618e-004 8.387707e-003 8.387707e-002第n 步计算值3、结果分析：很明显第二种递推式结果要比第一种好，式1在第七步后有明显误差，而式2在第三步后基本与近似解一致。
(完整word版)matlab数值分析例题
1、 在MATLAB 中用Jacobi 迭代法讨论线性方程组,1231231234748212515x x x x x x x x x -+=⎧⎪-+=-⎨⎪-++=⎩（1）给出Jacobi 迭代法的迭代方程，并判定Jacobi 迭代法求解此方程组是否收敛。
（2）若收敛，编程求解该线性方程组.解（1）：A=[4 -1 1;4 —8 1;-2 1 5］ %线性方程组系数矩阵A =4 -1 1 4 -8 1 —2 1 5>> D=diag(diag(A)）D =4 0 0 0 —8 0 0 0 5〉〉 L=—tril （A,-1) ％ A 的下三角矩阵L =0 0 0 —4 0 0 2 —1 0〉〉U=-triu(A，1）% A的上三角矩阵U =0 1 —10 0 —10 0 0B=inv（D）*（L+U）% B为雅可比迭代矩阵B =0 0.2500 —0。
25000.5000 0 0.12500。
4000 —0.2000 0〉〉r=eigs(B,1）%B的谱半径r =0。
3347 〈1Jacobi迭代法收敛。
(2）在matlab上编写程序如下:A=[4 —1 1;4 -8 1;—2 1 5]；〉〉b=［7 —21 15］'；>〉x0=［0 0 0]’;〉〉［x,k］=jacobi（A,b，x0，1e—7)x =2。
00004.00003。
0000k =17附jacobi迭代法的matlab程序如下：function [x,k］=jacobi（A,b，x0，eps)% 采用Jacobi迭代法求Ax=b的解％A为系数矩阵％b为常数向量％x0为迭代初始向量％eps为解的精度控制max1= 300; %默认最多迭代300，超过300次给出警告D=diag（diag（A）)；％求A的对角矩阵L=-tril(A，—1); %求A的下三角阵U=—triu（A,1); %求A的上三角阵B=D\(L+U);f=D\b;x=B＊x0+f；k=1；%迭代次数while norm（x-x0）>=epsx0=x;x=B＊x0+f；k=k+1;if（k〉=max1)disp(’迭代超过300次，方程组可能不收敛’）；return；endend2、设有某实验数据如下：（1)在MATLAB中作图观察离散点的结构,用多项式拟合的方法拟合一个合适的多项式函数；（2）在MATLAB中作出离散点和拟合曲线图。
MATLAB数值实验一(数据的插值运算及其应用完整版)
佛山科学技术学院实 验 报 告课程名称 数值分析 实验项目 插值法与数据拟合 专业班级 机械工程 姓 名 余红杰 学 号 2111505010 指导教师 陈剑 成 绩 日 期 月 日一、实验目的1、学会Lagrange 插值、牛顿插值和三次样条插值等基本插值方法；2、讨论插值的Runge 现象3、学会Matlab 提供的插值函数的使用方法，会用这些函数解决实际问题。
二、实验原理1、拉格朗日插值多项式2、牛顿插值多项式3、三次样条插值 三、实验步骤1、用MA TLAB 编写独立的拉格朗日插值多项式函数2、用MA TLAB 编写独立的牛顿插值多项式函数3、用MA TLAB 编写独立的三次样条函数（边界条件为第一、二种情形）4、已知函数在下列各点的值为：根据步骤1,2,3编好的程序，试分别用4次拉格朗日多项式4()L x 、牛顿插值多项式4()P x 以及三次样条函数()S x (自然边界条件)对数据进行插值，并用图给出 ｛(,),0.20.08,0,1,2,,10i i i x y x i i =+= ｝，4()L x 、4()P x 和()S x 。
5、在区间[-1,1]上分别取10,2n =用两组等距节点对龙格函数21(),(11)125f x x x =-≤≤+作多项式插值，对不同n 值，分别画出插值函数及()f x 的图形。
6、下列数据点的插值可以得到平方根函数的近似，在区间[0,64]上作图。
（1）用这9个点作8次多项式插值8()L x 。
（2）用三次样条（第一边界条件）程序求()S x 。
7、对于给函数21()125f x x =+在区间[-1,1]上取10.2(0,1,,10)i x i i =-+= ，试求3次曲线拟合，试画出拟合曲线并打印出方程，与第5题的结果比较。
四、实验过程与结果：1、Lagrange 插值多项式源代码：function ya=lag(x,y,xa) %x 所有已知插值点 %y 插值点对应函数值 %xa 所求点，自变量 %ya 所求点插值估计量 ya=0; mu=1; %初始化%循环方式求L 系数，并求和： for i = 1:length(y) for j = 1:length(x) if i ~= jmu = mu * (xa - x(j) ) / ( x(i) - x(j) ); else continue end endya = ya + y(i) * mu ; mu = 1; end2、Newton 源代码：function ya = newton(x,y,xa) %x 所有已知插值点 %y 插值点对应函数值 %xa 所求点，自变量 %ya 所求点插值估计量 %建立系数零矩阵D 及初始化：D = zeros(length(x)-1);ya = y(1);xi = 1;%求出矩阵D，该矩阵第一行为牛顿插值多项式系数：for i=1:(length(x)-1)D(i,1) = (y(i+1) -y(i))/(x(i+1) -x(i));endfor j=2:(length(x)-1)for i=1:(length(x)-j)D(i,j) = (D(i+1,j-1) - D(i,j-1)) / (x(i+j) - x(i)); endend%xi为单个多项式（x-x(1))(x-x(2))...的值for i=1:(length(x)-1)for j=1:ixi = xi*(xa - x(j));endya = ya + D(1,i)*xi;xi = 1;end3、三次样条插值多项式（1）（第一边界条件）源代码：function y=yt1(x0,y0,f_0,f_n,x) _____________(1)%第一类边界条件下三次样条插值；%xi 所求点；%yi 所求点函数值；%x 已知插值点；%y 已知插值点函数值；%f_0左端点一次导数值；%f_n右端点一次导数值；n = length(x0);z = length(y0);h = zeros(n-1,1);k=zeros(n-2,1);l=zeros(n-2,1);S=2*eye(n);for i=1:n-1h(i)= x0(i+1)-x0(i);endfor i=1:n-2k(i)= h(i+1)/(h(i+1)+h(i));l(i)= 1-k(i);end%对于第一种边界条件：k = [1;k]; _______________________(2)l = [l;1]; _______________________(3)%构建系数矩阵S：for i = 1:n-1S(i,i+1) = k(i);S(i+1,i) = l(i);end%建立均差表：F=zeros(n-1,2);for i = 1:n-1F(i,1) = (y0(i+1)-y0(i))/(x0(i+1)-x0(i));endD = zeros(n-2,1);for i = 1:n-2F(i,2) = (F(i+1,1)-F(i,1))/(x0(i+2)-x0(i));D(i,1) = 6 * F(i,2);end%构建函数D：d0 = 6*(F(1,2)-f_0)/h(1); ___________(4)dn = 6*(f_n-F(n-1,2))/h(n-1); ___________(5)D = [d0;D;dn]; ______________(6)m= S\D;%寻找x所在位置，并求出对应插值：for i = 1:length(x)for j = 1:n-1if (x(i)<=x0(j+1))&(x(i)>=x0(j))y(i) =( m(j)*(x0(j+1)-x(i))^3)/(6*h(j))+...(m(j+1)*(x(i)-x0(j))^3)/(6*h(j))+...(y0(j)-(m(j)*h(j)^2)/6)*(x0(j+1)-x(i))/h(j)+... (y0(j+1)-(m(j+1)*h(j)^2)/6)*(x(i)-x0(j))/h(j) ; break;else continue;endendend（2）（自然边界条件）源代码：仅仅需要对上面部分标注的位置做如下修改:__(1):function y=yt2(x0,y0,x)__(2):k=[0;k]__(3):l=[l;0]__(4)+(5):删除—(6):D=[0:D:0]4、——————————————PS:另建了一个f 方程文件，后面有一题也有用到。
数值分析matlab方法插值法
其中,
n
【注】
x [a, b] w(x) (x x j ) j 0
（1）误差估计
Rn (x)
M n1 (n 1)!
w( x)
M n1
max
x( a ,b )
f
(n1) (x)
（2）余项与 x、M n1 节点的位置、个数 n 有关
（3）当 f (x)是 n 的多项式时Ln (x) f (x) n
M2 2!
(x
x0 )(x
x1 )
其中，
M2
max
x( x0 , x1 )
f
(x)
x
[
6
,
4
]
，所以
R1
(
5
24
)
sin
2!
4 (5
24
)( 5
6 24
)
4
0.0061
2）
抛物插值误差估计.因为
R2 (x)
M3 3!
(x
x0 )(x
x1)(x
x2 )
其中，
M3
max
x( x0 ,x2 )
f (x)
yiynewtonbackwardxyxicos035yi09394数值分析插值法55埃尔米特插值2n12n2数值分析插值法551埃尔米特插值多项式的存在唯一性2n1数值分析插值法数值分析插值法552埃尔米特插值余项553三次埃尔米特插值多项式maxsinxsin1数值分析插值法56561高次插值的病态性质0908原函数150706050405分段线性插值0302543214321055数值分析插值法562分段低次插值方法563分段低次插值余项090807060504030201090807060504030201543214321数值分析插值法57三次样条插值571三次样条插值572三弯矩法数值分析插值法573三次样条插值的误差估计与收敛性58插值运算的matlab函数581一维插值函数interp1yiinterp1xyximethod?linear?yiinterp1xyxilinear?1200135019
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2016-2017 第一学期数值分析
上机实验报告
姓名： xxx
学号： 20162……．
学院：土木工程学院
导师：………..
联系电话：…………..
指导老师：………..
目录
第一题 (1)
1.1题目要求 (1)
1.2程序编写 (1)
1.3计算结果及分析 (2)
第二题 (4)
2.1题目要求 (4)
2.2程序编写 (4)
2.3计算结果及分析 (6)
第三题 (7)
3.1题目要求 (7)
3.2程序编写 (7)
3.3计算结果及分析 (8)
第四题 (9)
4.1题目要求 (9)
4.2程序编写 (9)
4.3计算结果及分析 (10)
第五题 (11)
5.1题目要求 (11)
5.2程序编写 (12)
5.3计算结果及分析 (13)
第六题 (17)
6.1题目要求 (17)
6.2程序编写 (17)
6.3计算结果及分析 (18)
6.4程序改进 (18)
第一题
选做的是第（1）小问。
1.1题目要求
编出不动点迭代法求根的程序；把??3+4??2-10=0写成至少四种??=g(??)的形式，取初值??0=1.5,进行不动点迭代求根，并比较收敛性及收敛速度。
1.2程序编写
1.3计算结果及分析
①第一种迭代公式：??=??3+4??2+??-10；matlab计算结果如下：（以下为命令行窗口的内容）
2；matlab计算结果如下：
②第二种迭代公式：??=√(10-??3)/4
（以下为命令窗口内容）
2；matlab计算结果如下：
③第三种迭代公式：??=√10(??+4)
?
（以下为命令窗口内容）
④第四种迭代公式：??=10(??2+4??)
?；matlab计算结果如下：
（以下为命令窗口内容）
上述4种迭代公式，1、4两种由于在x真实值附近|g`(x)|>1,不满足迭代局部收敛条件，所以迭代序列不收敛。
对于2、3两种式子，由于在x真实值附近|g`(x)|<=L<1,满足迭代局部收敛条件,所以迭代序列收敛。
对于2、3两迭代公式，由于L3<L2,所以第3个迭代公式比第2个迭代公式收敛更快。
第二题
选做的是第（2）小问
2.1题目要求
编写有效程序解线性方程组A??=b。
2.2程序编写
2.3计算结果及分析
（命令行窗口显示内容）数据太多截取一部分
第三题
3.1题目要求
在区间[-1,1]上取????=-1+0.2??（??=0,1,2?,10)，
对函数??(??)=1
1+25??2
（a）对函数进行多项式插值和三次样条插值，并画出插值函数及??(??)的函数；
（b）对函数求其三次拟合曲线并画出拟合曲线的图像，与（a）中结果进行比较。
3.2程序编写
3.3计算结果及分析
matlab画出的图像如下：
从上图可以清晰的看到，三次样条插值曲线与原始曲线最为接近，作为f(x)的插值函数比较理想。
而三次多项式拟合的话，与原始结果偏差很大。
第四题
选择的是第一小问
4.1题目要求
编写Gauss -Legendre 求积公式程序，并计算下列积分
（1）∫??-(????????)2????1-1
（2）∫??-(????????)2????20
（3）∫??10
|sin?(1???)|????（4）∫ln?(1+??)???1
0????4.2程序编写
4.3计算结果及分析
(a)∫??-(????????)2????
1-1（以下为命令行窗口的内容）
(b)∫??-(????????)2????
20（以下为命令行窗口的内容）
(c)∫??1
0|sin?(1???)|????（以下为命令行窗口的内容）
(d)∫log?(1+??)???1
0????（以下为命令行窗口的内容）
第五题
5.1题目要求
给定初值问题
??(??)={y`=-50y+50??2+2??,0≤x≤1
y(0)=13
????
用经典的四阶Lunge-Kutta方法求解，步长分别取为?=0.1,0.025,0.01，计算并打印??=0.1??(??=1,2,?,10)各点的值，并与准确解y(??)=1
3
??-50??+??2作比较。
5.2程序编写
5.3计算结果及分析
（以下为命令行窗口的内容）请输入步长h：0.1
从上面的结果可以看到，当步长为0.1时，求解是很不准确的，可能是步长过大导致的，所以应适当减小步长，再尝试进行计算。
下面的运算取步长为0.025，再进行试算。
（以下为命令行窗口的内容）
请输入步长h：0.025
从上面的结果可以看到，由于步长h取值减小,误差变得很小，可以看出其误差限是呈线性增加的。
（以下为命令行窗口的内容）
请输入步长h：0.01
从上面的结果可以看到，当步长取0.01时，误差减小了10多倍。
此外，我也尝试了采用更小的步长进行计算，计算所得的相对误差随步长减小也越来越小，
限于篇幅，也就不再一一罗列计算结果。
所以减小步长，可以很好地提高该算法的精确度，若我们把计算结果看成关于步长h的序列，该序列的收敛速度还比较快。
而在这一点上，我们也可以通过理论证明，当步长h趋于零时，计算结果就趋于真实结果，即收敛性。
第六题
6.1题目要求
编程计算∑1????
∞??=1,其中c =4.4942×10307，给出并观察计算结果，若有问题，分析之。
要求逐项相加编程求解，不允许使用
matlab 中的现有函数。
6.2程序编写
6.3计算结果及分析
根据我们所学的数学知识，这个关于n的求和序列实际是不收敛的，但是由于计算机计算时的大数吃掉了小数后，数值无法累加。
所以显示为上面的收敛结果。
6.4程序改进
上面的算法发生了大数吃小数的现象，为了避免大数吃小数，可以第一步设
n个数相加，然后把这些数分为两个一组，每组两个数相加，然后再分为两个一组，循环下去。
这样应该能避免大数吃小数。