下载文档原格式
MATLAB目录MATLAB的简介发展历程应用Matlab的优势和特点Matlab常用工具箱常用函数[编辑本段]MATLAB的简介MATLAB 是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLA B和Simulink两大部分。
MATLAB是矩阵实验室（Matrix Laboratory）的简称，和Mathematica、Mapl e并称为三大数学软件。
它在数学类科技应用软件中在数值计算方面首屈一指。
MAT LAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。
MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完相同的事情简捷得多，并且mathwork也吸收了像Maple等软件的优点,使MATLAB成为一个强大的数学软件。
在新的版本中也加入了对C，FORTRAN，C++ ，JAVA的支持。
可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。
[编辑本段]发展历程20世纪70年代，美国新墨西哥大学计算机科学系主任Cleve Moler为了减轻学生编程的负担，用FORTRAN编写了最早的MATLAB。
1984年由Little、Moler、St eve Bangert合作成立了的MathWorks公司正式把MATLAB推向市场。
到20世纪9 0年代，MATLAB已成为国际控制界的标准计算软件。
版本更新[编辑本段]应用MATLAB 产品族可以用来进行以下各种工作：● 数值分析● 数值和符号计算● 工程与科学绘图● 控制系统的设计与仿真● 数字图像处理● 数字信号处理● 通讯系统设计与仿真● 财务与金融工程MATLAB 的应用范围非常广，包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。
MATLAB与数值分析课程总结姓名：董建伟学号：2015020904027一：MATLAB部分1.处理矩阵－容易矩阵的创建（1）直接创建注意 a中括号里可以用空格或者逗号将矩阵元素分开 b矩阵元素可以是任何MATLAB表达式，如实数复数等c可以调用赋值过的任何变量，变量名不要重复，否则会被覆盖（2）用MATLAB函数创建矩阵如：a空阵［］ b rand/randn——随机矩阵 c eye——单位矩阵 d zeros ——0矩阵 e ones——1矩阵 f magic——产生n阶幻方矩阵等向量的生成（1）用冒号生成向量（2）使用linspace和logspace分别生成线性等分向量和对数等分向量矩阵的标识和引用（1）向量标识（2）“0 1”逻辑向量或矩阵标识（3）全下标，单下标，逻辑矩阵方式引用字符串数组（1）字符串按行向量进行储存（2）所有字符串用单引号括起来（3）直接进行创建矩阵运算（1）注意与数组点乘，除与直接乘除的区别，数组为乘方对应元素的幂（2）左右除时斜杠底部靠近谁谁是分母（3）其他运算如，inv矩阵求逆，det行列式的值, eig特征值，diag 对角矩阵2.绘图－轻松plot－绘制二维曲线（1）plot(x)绘制以x为纵坐标的二维曲线plot(x,y) 绘制以x为横坐标，y为纵坐标的二维曲线x,y为向量或矩阵（2）plot（x1,y1,x2,y2，。
）绘制多条曲线，不同字母代替不同颜色:b蓝色，y黄色，r红色，g绿色（3）hold on后面的pl ot图像叠加在一起hold off解除hold on命令，plot将先冲去窗口已有图形（4）在hold后面加上figure，可以绘制多幅图形（5）subplot在同一窗口画多个子图三维图形的绘制(1)plot3(x,y,z,’s’) s是指定线型，色彩，数据点形的字符串(2)［X,Y］＝meshgrid(x,y)生成平面网格点(3)mesh(x,y,z，c)生成三维网格点，c为颜色矩阵(4)三维表面处理mesh命令对网格着色，surf对网格片着色(5)contour绘制二维等高线(6)axis(［x1,xu,y1,yu］)定义x,y的显示范围3．编程－简洁（1）变量命名时可以由字母，数字，下划线，但是不得包含空格和标点（2）最常用的数据类型只有双精度型和字符型，其他数据类型只在特殊条件下使用（3）为得到高效代码，尽量提高代码的向量化程度，避免使用循环结构（4）为得到最快的运行速度，在循环指令前尽量对数组进行预定义（5）流程控制语句如a break命令可以使包含break的最内层的for或while语句强制终止，跳出循环结构，执行end后面的语句。
数值分析编程作业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 =。
数值分析中求解非线性方程的MATLAB求解程序1. fzero函数：fzero函数是MATLAB中最常用的求解非线性方程的函数之一、它使用了割线法、二分法和反复均值法等多种迭代算法来求解方程。
使用fzero函数可以很方便地求解单变量非线性方程和非线性方程组。
例如，要求解方程f(x) = 0，可以使用以下语法：``````2. fsolve函数：fsolve函数是MATLAB中求解多维非线性方程组的函数。
它是基于牛顿法的迭代算法来求解方程组。
使用fsolve函数可以非常方便地求解非线性方程组。
例如，要求解方程组F(x) = 0，可以使用以下语法：``````3. root函数：root函数是MATLAB中求解非线性方程组的函数之一、它采用牛顿法或拟牛顿法来求解方程组。
使用root函数可以非常方便地求解非线性方程组。
例如，要求解方程组F(x) = 0，可以使用以下语法：``````4. vpasolve函数：vpasolve函数是MATLAB中求解符号方程的函数。
它使用符号计算的方法来求解方程，可以得到精确的解。
vpasolve函数可以求解多变量非线性方程组和含有符号参数的非线性方程。
例如，要求解方程组F(x) = 0，可以使用以下语法：```x = vpasolve(F(x) == 0, x)```vpasolve函数会返回方程组的一个精确解x。
5. fsolve和lsqnonlin结合：在MATLAB中，可以将求解非线性方程转化为求解最小二乘问题的形式。
可以使用fsolve函数或lsqnonlin函数来求解最小二乘问题。
例如，要求解方程f(x) = 0，可以将其转化为最小二乘问题g(x) = min，然后使用fsolve或lsqnonlin函数来求解。
具体使用方法可以参考MATLAB官方文档。
6. Newton-Raphson法手动实现：除了使用MATLAB中的函数来求解非线性方程，还可以手动实现Newton-Raphson法来求解。
牛顿插值法是一种常用的数值分析方法，用于构造一个多项式函数，以便在给定的数据点上进行插值。
这个主题在数学和工程领域中有着广泛的应用，特别是在数据拟合和函数逼近方面。
牛顿插值法的核心思想是通过不断地添加新的数据点来构造一个多项式，并利用已知数据点来确定多项式的系数，从而实现对未知数据点的插值预测。
在Matlab中，实现牛顿插值法并不困难，我们可以利用已有的函数和工具来简化计算过程。
下面，我们将通过一个具体的例题来讲解如何使用Matlab编写牛顿插值法的程序，并分析其结果。
我们需要明确牛顿插值法的数学原理。
给定n个互不相同的节点\(x_0, x_1, ... , x_n\)，以及在这些节点上的函数值\(f(x_0), f(x_1), ... , f(x_n)\)，我们希望构造一个n次插值多项式p(x)，满足p(x_i) = f(x_i)，i=0,1,...,n。
牛顿插值多项式的一般形式为：\[p(x) = a_0 + a_1(x - x_0) + a_2(x - x_0)(x - x_1) + ... + a_n(x -x_0)(x - x_1)...(x - x_{n-1})\]其中，\[a_i\]表示插值多项式的系数。
通过牛顿插值法的迭代过程，可以逐步求解出这些系数，进而得到插值多项式的表达式。
接下来，我们将以一个具体的例题来演示如何在Matlab中实现牛顿插值法。
假设我们有如下的数据点和函数值：\(x = [1, 2, 3, 4]\)\(f(x) = [1, 4, 9, 16]\)我们希望利用这些数据点来构造一个插值多项式，并在给定的区间上进行插值计算。
在Matlab中，可以通过interp1函数来进行插值计算，该函数支持多种插值方法，包括牛顿插值法。
下面是一个简单的Matlab程序示例：```matlabx = [1, 2, 3, 4];y = [1, 4, 9, 16];xi = 2.5;yi = interp1(x, y, xi, 'spline');disp(['在x=',num2str(xi),'处的插值结果为：',num2str(yi)]);```在这段代码中，我们首先定义了给定的数据点x和对应的函数值y，然后利用interp1函数对x=2.5处的插值结果进行计算。
数值分析试验幂法与反幂法m a t l a b一、问题的描述及算法设计（一）问题的描述我所要做的课题是：对称矩阵的条件数的求解设计1、求矩阵A 的二条件数问题 A=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----210121012 2、设计内容：1）采用幂法求出A 的. 2）采用反幂法求出A 的.3）计算A 的条件数 ⅡA Ⅱ2* ⅡA -1Ⅱ2=cond2（A ）=/.（精度要求为10-6）3、设计要求1）求出ⅡA Ⅱ2。
2）并进行一定的理论分析。
（二）算法设计1、幂法算法（1）取初始向量u )0(（例如取u )0(=(1,1,…1)T ）,置精度要求ε，置k=1.（2）计算v )(k =Au )1(-k ,m k =max(v )(k ), u )(k = v )(k / m k（3）若| m k = m 1-k |<ε，则停止计算（m k 作为绝对值最大特征值1λ，u )(k 作为相应的特征向量）否则置k=k+1，转（2）2、反幂法算法（1）取初始向量u )0(（例如取u )0(=(1,1,…1)T ）,置精度要求ε，置k=1.（2）对A 作LU 分解，即A=LU（3）解线性方程组 Ly )(k =u )1(-k ,Uv )(k =y )(k（4）计算m k =max(v )(k ), u )(k = v )(k / m k（5）若|mk =m1-k|<ε，则停止计算（1/mk作为绝对值最小特征值nλ，u)(k作为相应的特征向量）;否则置k=k+1，转（3）.二、算法的流程图（一）幂法算法的流程图（二）反幂法算法的流程图三、算法的理论依据及其推导（一）幂法算法的理论依据及推导幂法是用来确定矩阵的主特征值的一种迭代方法，也即，绝对值最大的特征值。
稍微修改该方法，也可以用来确定其他特征值。
幂法的一个很有用的特性是它不仅可以生成特征值，而且可以生成相应的特征向量。
实际上，幂法经常用来求通过其他方法确定的特征值的特征向量。
程序代码：(代码标记[code]...[/code] ) 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');h3=[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)xlabel('实部','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实例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轴');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实例21：PEAKS函数曲线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轴');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 phongcolorbaraxis 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;x=0.05+i;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 on'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,',...];'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','高斯分布',...。
1.% 最小二乘法拟合数据点方法1：% 左除右除：xA=B ==> x=B/A | Ax=B ==> x=B\A% A 表示拟合函数的组合，如：多项式插值，A=[1,x,x.^2,...,x.^n],表示拟合函数为% 多项式：s(x)=a0+a1*x+...+an*x^n;又如：A=[log(x),cos(x),exp(x)]%则表示拟合函数为s(x)=a0*ln(x)+a1*cos(x)+a2*exp(x)% 法方程为：A'*A*z=A'*y ==> A*z=y ==> z=A\y z=(a0,...,an)'% Date: 2012-1-1clear;clc;x=[0 0.25 0.5 0.75 1]';y=[1 1.284 1.6487 2.1170 2.7183]';%索要拟合的数据点x1=ones(size(x),1);A=[x1 x x.^2];%拟合函数Z=A\y %A中每列函数的参数% 最小二乘法拟合数据方法2：采用polyfit函数% Date:2012-1-1clear;clc;x=[0 0.25 0.5 0.75 1]';y=[1 1.284 1.6487 2.1170 2.7183]';%索要拟合的数据点p=polyfit(x,y,2) %polyfit(x,y,n),(x,y)为数据点坐标，n为拟合多项式阶数，p 为% 所求拟合多项式的幂次从高到低排列的系数。
2.% 复合梯形公式计算积分值% 输入：fun--积分函数；a，b--积分区间；n--区间等分数% 输出：I--数值积分结果% 调用格式(ex)：re=ftrapz(@fun1,0,1,10)% 2012-1-1function I=ftrapz(fun,a,b,n)h=(b-a)/n;%区间等分x=linspace(a,b,n+1);%将a到b的区间等分成(n+1)-1个区间，数据点有(n+!)个y=feval(fun,x);I=h*(0.5*y(1)+sum(y(2:n))+0.5*y(n+1));%积分原函数% Date：2012-1-1function y=fun1(x)y=exp(-x);3.% 复合simpson公式求积分% 输入：fun--积分函数；a，b--积分区间；n--区间等分数% 输出：I--数值积分结果% 调用格式(ex)：re=fsimpson(@fun1,0,1,10)% 2012-1-1function I=fsimpson(fun,a,b,n)h=(b-a)/n;x=linspace(a,b,2*n+1);y=feval(fun,x);I=(h/6)*(y(1)+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n))+y(2*n+1));4.% 两点GS-Legendre公式求积分% 输入：fun--积分函数；a，b--积分区间；% 输出：I--数值积分结果% 调用格式(ex)：re=GSLege(@fun2,0,1)% 2012-1-1function I=GSLege(fun,a,b)%将区间[a.b]通过变量替换x=(a+b)/2-(b-a)/2*t变到[-1,1]%其中t取GS-Lege的高斯点m1=feval(fun,(a+b)/2+(b-a)/2*(-1/sqrt(3)));m2=feval(fun,(a+b)/2+(b-a)/2*(1/sqrt(3)));I=(b-a)/2*(m1+m2);% GS-Legendre公式的积分函数% Date:2012-1-1function y=fun2(x)y=sin(x)/x;5.% 追赶法求解线性方程组Ax=b，其中A是三对角方阵%function x=tridiagsolver(A,b)clear;clc;A=[2 -1 0 0;-1 3 -2 0;0 -2 4 -3;0 0 -3 5];%三对角矩阵，线性方程组系数矩阵b=[6,1,-2,1]';%[n,n]=size(A);for i=1:nif(i<2)l(i)=A(i,i);y(i)=b(i)/l(i);u(i)=A(i,i+1)/l(i);elseif i<nl(i)=A(i,i)-A(i,i-1)*u(i-1);y(i)=(b(i)-y(i-1)*A(i,i-1))/l(i);u(i)=A(i,i+1)/l(i);elsel(i)=A(i,i)-A(i,i-1)*u(i-1);y(i)=(b(i)-y(i-1)*A(i,i-1))/l(i);endendx(n)=y(n);for j=n-1:-1:1x(j)=y(j)-u(j)*x(j+1);endx6.% SOR迭代求解非线性方程组Ax=b% 输入：A--系数矩阵；b--；omega--松弛因子（0~2);tol--精度% 输出：x--方程的解向量;iter--迭代次数;% 调用格式(ex)：[x,iter=sor(A,b,1.1,1e-4)% 2012-1-1%function [x,iter]=sor(A,b,omega,tol)%{%}clear;clc;A=[2 -1 0;-1 3 -1;0 -1 2];b=[1 8 -5]';omega=1.1;tol=1e-4;D=diag(diag(A));%diag(A)返回的是A的对角元组成的列向量；diag(b)返回的是以列向量b为对角元的方阵；L=D-tril(A);%tril(A)返回的是矩阵A的下三角矩阵；U=D-triu(A);%triu(A)返回的是矩阵A的上三角矩阵；x=zeros(size(b));%给定迭代初始值零向量iter=1;while iter<500x=(D-omega*L)\(omega*b+(1-omega)*D*x+omega*U*x);%迭代格式error=norm(b-A*x)/norm(b);%收敛条件if error<tolbreak;enditer=iter+1;endif iter>= 500fprintf('root not found!');end7.% 牛顿法求解非线性方程的根% 输入：fun--非线性函数；dfun--非线性函数导数；x0--初始值；tol--精度；% 输出：x--非线性方程数值根% 调用格式(ex)：x=newton(@fun3,@dfun3,6,1e-3)% 2012-1-1function x=newton(fun,dfun,x0,tol)iter=1;if abs(feval(fun,x0))<tolx=x0;return;endwhile iter<500if iter==1x=x0-feval(fun,x0)/feval(dfun,x0);if abs(feval(fun,x))<tolbreak;endelsex=x-feval(fun,x)/feval(dfun,x);if abs(feval(fun,x))<tolbreak;endenditer=iter+1;endif iter==500fprintf('not successful!');x=NaN;end% newton的函数文件% Date：2012-1-1function y=fun3(x)y=3.*x.^3-8.*x.^2-8.*x-11;% newton的导函数文件% Date：2012-1-1function y=dfun3(x)y=9.*x.^2-16.*x-8;8.% 两点割线法求解非线性方程的根% 输入：fun--非线性函数；a,b--两个初始值；tol--精度；% 输出：x--非线性方程数值根% 调用格式(ex)：x=gexian(@fun3,3,6,1e-3)% 2012-1-1function x=gexian(fun,a,b,tol)iter=1;xk1=a;xk2=b;while iter<500xk3=xk2-feval(fun,xk2)*(xk2-xk1)/(feval(fun,xk2)-feval(fun,xk1));if abs(feval(fun,xk3))<tolbreak;elsexk1=xk2;xk2=xk3;enditer=iter+1;endif iter==500fprintf('not successful!');x=NaN;elsex=xk3;end9.% 乘幂法求矩阵的按模最大特征值及其特征向量% 输入：A--要求的矩阵；v0--初始非零向量；tol--精度；% 输出：x--特征向量；lam--按模最大特征值% 调用格式(ex)：[lam,x]=power(A,v0,1e-2)% 2012-1-1function [lam,x]=matrixpower(A,v0,tol)v1=A*v0;v2=A*v1;sum=0;p=0;for i=1:size(v1)if v1(i)~=0sum=sum+v2(i)/v1(i);p=p+1;endendlam0=sum/p;iter=2;vk1=v2;while iter<500vk2=A*vk1;sum=0;p=0;for i=1:size(vk1)if vk1(i)~=0sum=sum+vk2(i)/vk1(i);p=p+1;endendlam=sum/p;if abs(lam-lam0)<tolbreak;elselam0=lam;vk1=vk2;endendif iter==500fprintf('not successful!');lam=NaN;elsex=vk2;end9_2% 改进乘幂法求矩阵的按模最大特征值及其特征向量% 输入：A--要求的矩阵；v0--初始非零向量；tol--精度；% 输出：x--特征向量；lam--按模最大特征值% 调用格式(ex)：[lam,x]=pMatrixPower(A,v0,1e-2)% Date:2012-1-2function [lam,x]=pMatrixPower(A,v0,tol)[ty,ti]=max(abs(v0));%返回v0中元素绝对值最大的元素值与下标,ti为下标lam0=v0(ti);u0=v0/lam0;iter=1;while iter<500v1=A*u0;[tv,ti]=max(abs(v1));lam1=v1(ti);u0=v1/lam1;if abs(lam0-lam1)<tolbreak;elselam0=lam1;iter=iter+1;endendif iter>=500fprintf('not successful!');lam=NaN;elselam=lam1;x=u0;end10.% 反幂法求矩阵的按模最小特征值及其特征向量% 输入：A--要求的矩阵；v0--初始非零向量；tol--精度；% 输出：x--特征向量；lam--按模最大特征值% 调用格式(ex)：[lam,x]=invMatrixPower(A,v0,1e-2)% Date:2012-1-2function [lam,x]=invMatrixPower(A,v0,tol)[ty,ti]=max(abs(v0));lam0=v0(ti);u0=v0/lam0;iter=1;while iter<500v1=A\u0;[tv,ti]=max(abs(v1));lam1=v1(ti);u0=v1/lam1;if abs(1/lam0-1/lam1)<tolbreak;elselam0=lam1;iter=iter+1;endendif iter>=500fprintf('not successful!');lam=NaN;elselam=1/lam1;x=u0;end11.% 欧拉方法求解一阶常微分方程初值问题% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;% 输出：y--常微分方程在区间[a,b]上各点的数值解；% 调用格式(ex)：y=odeEuler(@fun4,0,1,1,10)% Date:2012-1-2function y=odeEuler(fun,a,b,y0,n)h=(b-a)/n;y(1)=y0;x=a:h:b;for i=1:ny(i+1)=y(i)+h*feval(fun,x(i),y(i));end% 常微分方程fun4=0% Date：2012-1-2function re=fun4(x,y)re=y-2*x/y;12.% 改进欧拉公式（预估--校正）方法求解一阶常微分方程初值问题% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;% 输出：y--常微分方程在区间[a,b]上各点的数值解；% 调用格式(ex)：y=pOdeEuler(@fun4,0,1,1,10)% Date:2012-1-2function y=pOdeEuler(fun,a,b,y0,n)h=(b-a)/n;y(1)=y0;x=a:h:b; for i=1:nyp=y(i)+h*feval(fun,x(i),y(i));yc=y(i)+h*feval(fun,x(i+1),yp);y(i+1)=0.5*(yp+yc);end13.% 梯形方法求解一阶常微分方程初值问题% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;% 输出：y--常微分方程在区间[a,b]上各点的数值解；% 调用格式(ex)：y=trapezium(@fun4,0,1,1,10)% Date:2012-1-2function y=trapezium(fun,a,b,y0,n)h=(b-a)/n;y(1)=y0;x=a:h:b;tol=1e-6;for i=1:n%用不动点迭代的方法求解非线性方程：%y(i+1)=y(i)+h*feval(fun,x(i),y(i))/2+h*feval(fun,x(i),y(i+1))/2;iter=1;yy0=1+(i-1)*h;%迭代初始值while iter<500yy1=y(i)+h*feval(fun,x(i),y(i))/2+h*feval(fun,x(i)+h,yy0)/2;if abs(yy1-yy0)<tolbreak;elseyy0=yy1;iter=iter+1;endendif iter>=500printf('not successful!');y=NaN;return;elsey(i+1)=yy1;endend14.% 标准四阶四段龙格-库塔方法求解一阶常微分方程初值问题% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;% 输出：y--常微分方程在区间[a,b]上各点的数值解；% 调用格式(ex)：y=longekuta(@fun4,0,1,1,10)% Date:2012-1-2function y=longekuta(fun,a,b,y0,n)h=(b-a)/n;y(1)=y0;for k=2:n+1x=a+(k-2)*h;k1=h*feval(fun,x,y(k-1));k2=h*feval(fun,x+h/2,y(k-1)+k1/2);k3=h*feval(fun,x+h/2,y(k-1)+k2/2);k4=h*feval(fun,x+h,y(k-1)+k3);y(k)=y(k-1)+(k1+2*k2+2*k3+k4)/6;end15.插值% 拉格朗日插值% 输入：x,y--插值数据点（x，y均为行向量）；xh--要插值的点; % 输出：yh--插值结果；% 调用格式(ex)：yh=lagrange(x,y,xh)% Date:2012-1-2function yh=lagrange(x,y,xh)n=length(x);m=length(xh);yh=zeros(1,m);c1=ones(n-1,1);c2=ones(1,m);for i=1:nxp=x([1:i-1 i+1:n]);yh=yh+y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2));%prod对输入的一个向量返回其所有分量的乘积end%拉格朗日调用clear;clc;x=[11 12];y=[2.3979 2.4849];xh=11.75;yh=lagrange(x,y,xh)15_2% 牛顿插值% 输入：x,y--插值数据点（x，y均为行向量）；xh--要插值的点; % 输出：yh--插值结果；% 调用格式(ex)：yh=newtonPol(x,y,xh)% Date:2012-1-2function yh=newtonPol(x,y,xh)n=length(x);p(:,1)=x;p(:,2)=y;for j=3:n+1p(1:n+2-j,j)=diff(p(1:n+3-j,j-1))./(x(j-1:n)-x(1:n+2-j))';%求差商表endq=p(1,2:n+1)';%求牛顿法的系数--取第一行yh=0;m=1;yh=q(1);for i=2:nm=q(i);for j=2:im=m*(xh-x(j-1));%求牛顿法中各多项式值(xh-x0)…(xh-x n-1) endyh=yh+m;%求和end%牛顿插值调用clear;clc;x=[11 12 13];y=[2.3979 2.4849 2.5649];xh=11.75;yh=newtonPol(x,y,xh)。
