MATLAB实现计算全息
用MATLAB 软件和液晶光阀实现傅立叶变换计算全息制作及其再现
姚雪灿
指导教师 阎晓娜
（上海大学理学院物理系，上海 200444）
摘要：利用MATLAB 语言制作了一个迂回相位编码的傅立叶变换全息图，使用电寻址的液晶光阀作为全息图的实时记录介质对得到的傅立叶计算全息图进行光学再现，并对编码过程中加随机相位和不加随机相位后的再现图进行了比较讨论。
关键词：计算全息 傅立叶变换全息 MATLAB 液晶光阀 迂回相位编码
全息制作包括二种方式，光学全息和计算全息。光学全息用光学干涉原理制作，计算全息是用计算机对物波场的数学描述进行抽样、计算、编码而制作。计算全息可以制作已存在物体的全息图，也可以制作不存在物体的全息图，只要物光波场可以用数学描述出来。制作的计算全息图要以适合光学再现的尺寸和方式来输出。由于计算全息图上每个抽样单元的尺寸在微米量级，需要专门的光学缩微照相系统或微光刻系统。在要求较低情况下也可用照相机将显示在计算机屏幕或打印输出的计算全息原图缩拍到高分辨感光胶片上，通过显影、定影等处理得到可用于光学再现的全息图。由于记录介质是照相胶片，这就限制了它在实时处理中的应用。 近年来，随着高分辨电寻址空间光调制器的发展，像元尺寸在微米量级，像素数超过100万的振幅型或相位型空间光调制器已经完全实用化。其中最具代表性的是液晶光阀，电寻址的液晶光阀是由驱动电路驱动的LCD ，根据寻址电信号改变每一液晶像素的透过率，从而把电信号转换成空间的光强分布。液晶光阀可以作为实时的信号处理和显示器件，代替全息干板可进行实现计算全息图的实时输出和再现。
本文提出一种利用电寻址液晶光阀作为实时记录介质的计算机制全息图的产生方法，实验结果证明了这种方法的可行性。
1 用Matlab 软件实现傅立叶变换计算全息图
傅立叶变换全息记录的复数波面是物光波的傅立叶变换。计算傅立叶变换全息图的制作包括：对物光波抽样、离散傅立叶变换、编码、画图、图像的输出。在制作全息图的过程中,编码是最关键的一步，通过编码把二维光场的复振幅分布变换为全息图的二维透过率分布。本文以迂回相位编码来介绍编码过程。
设抽样后物光波的复振幅经过离散傅里叶变换后的频谱分布为复数F(m,n), 记为
F(m,n) = R(m,n)+iI(m,n), F(m,n) = A(m,n)·exp[i φ(m,n)] (1)
其中, A(m,n)和φ(m,n)分别代表全息图上各点的幅值和相位,
A(m,n) =),(),(22n m I n m R +, φ(m,n) =arctg[I(m,n)/R(m,n)] (2)
由于光学模板的最大透过率为1,所以在编码前还应对A(m,n)的值进行归一化,使其最大值为
1。假定将物面分为N×N 个抽样单元, 抽样间距为δx 和δy, 其间距要遵循Nyquist 判据。采用罗曼Ⅲ型编码方法,通过改变每个抽样单元内通光孔径的面积来编码振幅,通过改变通光孔径中心与抽样单元中心的位置来编码相位。最后每个像素用一个矩形孔表示，矩形孔的宽度为Wδx, 其中W 为一常数。矩形孔径的高度为Lmnδy,与归一化振幅成正比, Pmnδx 是孔径中心与单元中心的距离,并与抽样点的位相成正比。孔径参数与复值函数的关系如下，
mn L =mn A , mn P =mn φ/2πK （3）
经过计算,取W =1/2, K =1。
根据以上二元傅里叶变换全息图的实现原理，采用以下的算法思想在MATLAB 中进行二元傅里叶变换计算全息图的制作，采用罗曼Ⅲ型编码方式且以字母K 为例。其编码如下： b=zeros(128,128); %采样点阵为128X128。
b(8:120,24:40)=1;
b(8+M,(96-M):(112-M))=1;
end
for N=1:56
b(64+N,(40+N):(56+N))=1;
end
%K图形用矩阵表示如上。也可选取简单图形F表示。
A=b;a=rand(128,128);
aa=exp(i*2*pi.*a);
％给矩阵图K一个随机相位，随机相位因子的作用是平滑傅立叶频谱。
AA=double(A).*aa;
Afft2=fft2(AA);
％用matlab工具箱对带有随机相位因子的K图进行快速傅立叶变换。
A1=abs(Afft2);
B2=angle(Afft2)/(2*pi);%对频谱的相位进行归一化
A1max=max(max(A1));
A1=A1/A1max; %对频谱的幅值进行归一化
s=1; ％设定正方形单元的边长
figure;
axis([0 128 0 128]);
hold on;
for J=1:128
y0=s/2+(J-1)*s;
for I=1:128
x0=s/2+(I-1)*s; ％计算抽样单元的中心位置
H=A1(J,I)*s; ％矩形高度直接等于归一化的频谱幅值
F1=B2(J,I)*s; ％偏离单元中心的量
W=s/2; ％矩形宽度
if abs (F1)<=s/4
x2=x0-W/2; x3=x0+W/2; y2=y0-H/2; y3=y0+H/2;
fill([x2,x2,x3,x3],[y2,y3,y3,y2],'k');
else if F1>s/4
x22=x0+F1-W/2;x33=x0+s/2;
y22=y0-H/2;y33=y0+H/2;
fill([x22,x22,x33,x33],[y22,y33,y33,y22],'k');
x222=x0-s/2;x333=x0+F1+W/2-s;
fill([x222,x222,x333,x333],[y22,y33,y33,y22],'k');
%根据迂回相位编码原理，对φ<π／2，则直接用单元的高度表达频谱的幅值，偏移中心量表达频谱的相位。
else x22=x0-s/2;x33=x0+F1+W/2;
y2=y0-H/2;y3=y0+H/2;
fill([x22,x22,x33,x33],[y2,y3,y3,y2],'k');
x222=x0+F1-W/2+s;x333=x0+s/2;
fill([x222,x222,x333,x333],[y2,y3,y3,y2],'k');
%对开孔矩形的位置偏离中心超过，则在进行位相编码时，当φ>π／2时，为了防止与邻近的矩形
%孔发生重叠而造成全息图再现时产生失真，依据光栅衍射理论，在程序中采用了“模式溢出校%正法”,即将溢出部分移至本抽样单元的另一侧。
end
end
end
axis('equal');
axis off;
由上述程序运行得到的CGH图形如右：
图1 迂回相型计算机制全息图（部分）
2 用液晶光阀作为计算全息图记录介质实现光学再现
液晶光阀是利用液晶的光学特性而制作的空间光调制器。液晶是一种有机高分子化合物,既有晶体的取向性,又有液体的流动性。当液晶分子有序排列时表现出光学各向异性:光矢量沿分子长轴方向时具有较大的非常光折射率n e;而垂直分子长轴方向为寻常光折射率n o。若我们把两块玻璃合在一起,中间用一定厚度的间隔层控制玻璃间的距离,再在间隔中充满液晶,便形成一液晶盒。液晶盒玻璃内表面经一定方法处理后,可以使盒中的液晶分子长轴沿表面的某一方向排列。此时的液晶盒相当于一块晶轴与表面平行的晶体片,晶轴方向即为分子长轴方向。若我们在组成液晶盒的两玻璃间加一定电压,盒里的液晶分子在电场的作用下会逐渐沿电场方向排列,即光轴向与表面垂直的方向偏转,其偏转程度与电场强度有关。当一束直线偏振光通过此液晶盒时,出射光的偏振态将由加在液晶盒上的电压决定.由此实现了电场对光的偏振态控制。液晶光阀正是利用此特点制作的器件。
一般的计算全息图像的输出使用照相缩小后再现，为了实现实时化处理，采用液光阀代替照相胶片作为记录介质。实现再现的光路如图2。如图2所示，透射型电寻址液晶光阀与计算机视频输出联结，接受其调制信号。计算机输出全息图的电信号到液晶光阀上，由驱动的LCD根据寻址电信号改变每一个液晶像素的透过率，从而把电信号转换成空间的光强分布。其中液晶光阀采用电寻址的液晶光阀。
图2 使用液晶光阀的计算全息再现光路图
如图2所示，透射型电寻址液晶光阀与计算机视频输出联结，接受其调制信号。计算机输出全息图的电信号加到液晶光阀上驱动电路上，根据寻址电信号改变每一个液晶像素的透过率，从而把电信号转换成空间的光强分布。
图3(a,b) 是用Matlab编码制作的傅立叶变换计算全息图通过液晶光阀实时显示，用光学再现得到的再现图像。其中图3(a)和3(b)分别是编码过程中对K字母矩阵加随机相位和不加随机相位后再现图像。由图可见，当加随机相位进行编码后所得到的图形较清晰。图b则只显现出字母
K图形的高频成分，而低频成分被滤去。由此可推断，在对全息图进行编码时，如果对其图形未加入随机相位，那么作出的全息图又具有高通滤波器的功能，最后通过液晶光阀所成的重现像突出了边缘，这在图像处理中又称作边缘增强。关于边缘增强的原理有待进一步讨论。
（a）加随机相位（b）不加随机相位
图3 傅立叶变换计算机全息图字母K的重现像
3、结束语
本文通过运用MATLAB软件和液晶光阀对傅立叶计算全息进行编码和再现，并通过在编码时对矩阵加随相位和不加随机相位，得到两种截然不同的效果，看到了边缘增强的效果。实践证明，利用一些现代化工具会对物理的研究更加有效，明朗，快捷。
参考文献：
[1] 苏显渝等信息光学第一版北京科学出版社1999年第159—175页。
[2] 虞祖良，金国藩计算机制全息学第一版北京清华大学出版社1984年第31—41页。
[3] 宋菲君,Jutamulia S.近代光学信息处理[M] 北京大学出版社1998年第33—34页
[4] 陈炳初,雷宏香,蔡志岗, 用MATLAB制作迂回相位编码的傅立叶变换计算全息图,中山大学研究生学刊(自然科学、医学版),27(2):60-65,2006.
《应用计算方法教程》matlab作业二
6-1 试验目的计算特征值，实现算法 试验容：随机产生一个10阶整数矩阵，各数均在-5和5之间。 （1） 用MATLAB 函数“eig ”求矩阵全部特征值。 （2） 用幂法求A 的主特征值及对应的特征向量。 （3） 用基本QR 算法求全部特征值（可用MATLAB 函数“qr ”实现矩阵的QR 分解）。 原理 幂法：设矩阵A 的特征值为12n ||>||||λλλ≥???≥并设A 有完全的特征向量系12,,,n χχχ???(它们线性无关)，则对任意一个非零向量0n V R ∈所构造的向量序列1k k V AV -=有11()lim ()k j k k j V V λ→∞ -=， 其中()k j V 表示向量的第j 个分量。 为避免逐次迭代向量k V 不为零的分量变得很大（1||1λ>时）或很小（1||1λ<时），将每一步的k V 按其模最大的元素进行归一化。具体过程如下： 选择初始向量0V ，令1max(),,,1k k k k k k k V m V U V AU k m +===≥，当k 充分大时1111,max()max() k k U V χλχ+≈ ≈。 QR 法求全部特征值： 111 11222 111 ,1,2,3,k k k k k A A Q R R Q A Q R k R Q A Q R +++==????==??=???? ??????==?? 由于此题的矩阵是10阶的，上述算法计算时间过长，考虑采用改进算法——移位加速。迭 代格式如下： 1 k k k k k k k k A q I Q R A R Q q I +-=?? =+? 计算k A 右下角的二阶矩阵() () 1,1 1,() (),1 ,k k n n n n k k n n n n a a a a ----?? ? ??? 的特征值()()1,k k n n λλ-，当()()1,k k n n λλ-为实数时，选k q 为()()1,k k n n λλ-中最接近(),k n n a 的。 程序
计算方法_全主元消去法_matlab程序
%求四阶线性方程组的MA TLAB程序 clear Ab=[0.001 2 1 5 1; 3 - 4 0.1 -2 2; 2 -1 2 0.01 3; 1.1 6 2.3 9 4];%增广矩阵 num=[1 2 3 4];%未知量x的对应序号 for i=1:3 A=abs(Ab(i:4,i:4));%系数矩阵取绝对值 [r,c]=find(A==max(A(:))); r=r+i-1;%最大值对应行号 c=c+i-1;%最大值对应列号 q=Ab(r,:),Ab(r,:)=Ab(i,:),Ab(i,:)=q;%行变换 w=Ab(:,c),Ab(:,c)=Ab(:,i),Ab(:,i)=w;%列变换 n=num(i),num(i)=num(c),num(c)=n;%列变换引起未知量x次序变化for j=i:3 Ab(j+1,:)=-Ab(j+1,i)*Ab(i,:)/Ab(i,i)+Ab(j+1,:);%消去过程 end end %最后得到系数矩阵为上三角矩阵 %回代算法求解上三角形方程组 x(4)=Ab(4,5)/Ab(4,4); x(3)=(Ab(3,5)-Ab(3,4)*x(4))/Ab(3,3); x(2)=(Ab(2,5)-Ab(2,3)*x(3)-Ab(2,4)*x(4))/Ab(2,2); x(1)=(Ab(1,5)-Ab(1,2)*x(2)-Ab(1,3)*x(3)-Ab(1,4)*x(4))/Ab(1,1); for s=1:4 fprintf('未知量x%g =%g\n',num(s),x(s)) end %验证如下 %A=[0.001 2 1 5 1; 3 -4 0.1 -2 2;2 -1 2 0.01 3; 1.1 6 2.3 9 4]; %b=[1 2 3 4]'; %x=A\b; %x1= 1.0308 %x2= 0.3144 %x3= 0.6267 %x4= -0.0513
实验五_MATLAB计算的可视化
实验五 MATLAB 计算的可视化（一） 实验目的 1. 熟练掌握MATLAB 二维曲线的绘制 2.掌握图形的修饰 3.掌握三维图的绘制 4.了解各种特殊图形的绘制 内容与步骤 1.在同一幅图形窗口中分别绘制y1=sin(t)和y2=cos(t)二条函数曲线，t 的取值范围为[0,10]。y1用红色虚线表示，y2用蓝色实线表示，横坐标轴名称为“时间t ”，纵坐标轴名称为“正弦、余弦”，整个图形的标题为“正弦和余弦曲线”。在坐标(1.7*pi ，-0.3)处添加文字“sin(t)”， 在坐标(1.6*pi ，0.8)处添加文字“cos(t)”，并在右上角添加图例，其运行界面图如下图所示。之后并尝试修改坐标轴刻度。 2.用subplot 命令在同一个窗口的不同子窗口绘制曲线y=sin(t)，y1=sin(t+0.25) y2=sin(t+0.5)，其中t=[0 10]。 3.绘制三维曲线：?? ? ??=≤≤==)cos()sin()200() cos()sin(t t t z t t y t x π （注意：用plot3命令） 4.三维网线图：绘制z=sin(y)cos(x) 三维网线图。 5. 三维曲面图 绘制22y x z +=的三维曲面图，x 在[-5,5]范围，y 在[-5,5]范围。将曲面图颜色用shading 命令连续变化，并用颜色标尺显示色图（使用函数colorbar 生成）。生成的图形如下图所示。
6.请绘制一个饼形图，数据如下表所示 7. 用semilogx命令绘制传递函数为1//(s+1)(0.5s+1)的对数幅频特性曲线，横坐标为w，纵坐标为Lw，w的范围为10-2-103，按对数分布。
王能超 计算方法——算法设计及MATLAB实现课后代码
第一章插值方法 1.1Lagrange插值 1.2逐步插值 1.3分段三次Hermite插值 1.4分段三次样条插值 第二章数值积分 2.1 Simpson公式 2.2 变步长梯形法 2.3 Romberg加速算法 2.4 三点Gauss公式 第三章常微分方程德差分方法 3.1 改进的Euler方法 3.2 四阶Runge-Kutta方法 3.3 二阶Adams预报校正系统 3.4 改进的四阶Adams预报校正系统 第四章方程求根 4.1 二分法 4.2 开方法 4.3 Newton下山法 4.4 快速弦截法 第五章线性方程组的迭代法 5.1 Jacobi迭代 5.2 Gauss-Seidel迭代 5.3 超松弛迭代 5.4 对称超松弛迭代 第六章线性方程组的直接法 6.1 追赶法 6.2 Cholesky方法 6.3 矩阵分解方法 6.4 Gauss列主元消去法
第一章插值方法 1.1Lagrange插值 计算Lagrange插值多项式在x=x0处的值. MATLAB文件：（文件名：Lagrange_eval.m）function [y0,N]= Lagrange_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Lagrange插值多项式在x0处的值 %N是Lagrange插值函数的权系数 m=length(X); N=zeros(m,1); y0=0; for i=1:m N(i)=1; for j=1:m if j~=i; N(i)=N(i)*(x0-X(j))/(X(i)-X(j)); end end y0=y0+Y(i)*N(i); end 用法》X=[…];Y=[…]; 》x0= ; 》[y0,N]= Lagrange_eval(X,Y,x0) 1.2逐步插值 计算逐步插值多项式在x=x0处的值. MATLAB文件：（文件名：Neville_eval.m）function y0=Neville_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Neville逐步插值多项式在x0处的值 m=length(X); P=zeros(m,1); P1=zeros(m,1); P=Y; for i=1:m P1=P; k=1; for j=i+1:m k=k+1;
实验4、matlab的计算可视化和GUI设计
p345 subplot(2,2,1) t1=0:0.1:2; y1=sin(2*pi*t1); plot(t1,y1); title('y=sin(2\pit)') 练习： subplot(2,2,2) t2=0:0.1:2; y2=[exp(-t2);exp(-2*t2);exp(-3*t2)]; plot(t2,y2) axis([0 2 -0.2 1.2]); title('y=e-t,y=e-2t,y=e-3t') 练习： subplot(2,2,3); t3=[0 1 1 2 2 3 4]; y3=[0 0 2 2 0 0 0]; plot(t3,y3); axis([0 4 -0.5 3]); title('脉冲信号') 练习： subplot(2,2,4); t4=0:0.1:2*pi; plot(sin(t4),cos(t4));
axis([-1.2 1.2 -1.2 1.2]); axis equal; title('圆') 练习： P346 x=0:0.1:20; zeta=0 y1=1-1/sqrt(1-zeta^2)*exp(-zeta*x).*sin(sqrt( 1-zeta^2)*x+acos(zeta)); plot(x,y1) zeta=0.3; y2=1-1/sqrt(1-zeta^2)*exp(-zeta*x).*sin(sqrt( 1-zeta^2)*x+acos(zeta)); hold on plot(x,y2,'r:') zeta=0.5; y3=1-1/sqrt(1-zeta^2)*exp(-zeta*x).*sin(sqrt( 1-zeta^2)*x+acos(zeta)); plot(x,y3,'g*') zeta=0.707; y4=1-1/sqrt(1-zeta^2)*exp(-zeta*x).*sin(sqrt( 1-zeta^2)*x+acos(zeta)); plot(x,y4,'m-') title('二阶系统曲线') legend('\zeta=0','\zeta=0.3','\zeta=0.5','\zeta=0. 707') grid on gtext('\zeta=0') gtext('\zeta=0.3') gtext('\zeta=0.5') gtext('\zeta=0.707') ginput(3) zeta = ans = 2.6037 0.9035 13.1106 2.0029 4.2166 1.0380 P347 h_fig=gcf h_axis=gca h_line1=gco h_title=get(gca,'title') h_text2=findobj(h_fig,'string','\zeta=0.3') h_fig = 1 h_axis = 151.0018 h_line1 = 1 h_title = 152.0018 h_text2 = Empty matrix: 0-by-1 set(h_line1,'linewidth',5)
(整理)matlab16常用计算方法.
常用计算方法 1．超越方程的求解 一超越方程为 x (2ln x – 3) -100 = 0 求超越方程的解。 [算法]方法一：用迭代算法。将方程改为 01002ln()3 x x =- 其中x 0是一个初始值，由此计算终值x 。取最大误差为e = 10-4，当| x - x 0| > e 时，就用x 的值换成x 0的值，重新进行计算；否则| x - x 0| < e 为止。 [程序]P1_1abs.m 如下。 %超越方程的迭代算法 clear %清除变量 x0=30; %初始值 xx=[]; %空向量 while 1 %无限循环 x=100/(2*log(x0)-3); %迭代运算 xx=[xx,x]; %连接结果 if length(xx)>1000,break ,end %如果项数太多则退出循环（暗示发散） if abs(x0-x)<1e-4,break ,end %当精度足够高时退出循环 x0=x; %替换初值 end %结束循环 figure %创建图形窗口 plot(xx,'.-','LineWidth',2,'MarkerSize',12)%画迭代线'.-'表示每个点用.来表示，再用线连接 grid on %加网格 fs=16; %字体大小 title('超越方程的迭代折线','fontsize',fs)%标题 xlabel('\itn','fontsize',fs) %x 标签 ylabel('\itx','fontsize',fs) %y 标签 text(length(xx),xx(end),num2str(xx(end)),'fontsize',fs)%显示结果 [图示]用下标作为自变量画迭代的折线。如P0_20_1图所示，当最大误差为10-4时，需要迭代19次才能达到精度，超越方程的解为27.539。 [算法]方法二：用求零函数和求解函数。将方程改为函数 100()2ln()3f x x x =-- MATLAB 求零函数为fzero ，fzero 函数的格式之一是 x = fzero(f,x0) 其中，f 表示求解的函数文件，x0是估计值。fzero 函数的格式之二是 x = fzero(f,[x1,x2])
matlab计算结果的可视化
第五讲计算结果的可视化 本节介绍MATLAB 的两种基本绘图功能：二维平面图形和三维立体图形。 5.1 二维平面图形 5.1.1 基本图形函数 plot 是绘制二维图形的最基本函数，它是针对向量或矩阵的列来绘制曲线的。也就是说，使用plot 函数之前，必须首先定义好曲线上每一点的x 及y 坐标，常用格式为：（1）plot(x) 当x 为一向量时，以x 元素的值为纵坐标，x 的序号为横坐标值绘制 曲线。当x 为一实矩阵时，则以其序号为横坐标，按列绘制每列元素值相对于其序号的曲线， 当x 为m× n 矩阵时，就由n 条曲线。 （2）plot(x,y) 以x 元素为横坐标值，y 元素为纵坐标值绘制曲线。 （3）plot(x,y1,x,y2,…) 以公共的x 元素为横坐标值，以y1,y2,… 元素为纵坐标值绘 制多条曲线。 例5.1.1 画出一条正弦曲线和一条余弦曲线。 >> x=0:pi/10:2*pi; >> y1=sin(x); >> y2=cos(x); >> plot(x,y1,x,y2) 图5.1.1 函数plot 绘制的正弦曲线 在绘制曲线图形时，常常采用多种颜色或线型来区分不同的数据组，MATLAB 软件专 门提供了这方面的参数选项（见表5.1.1），我们只要在每个坐标后加上相关字符串，就可实 现它们的功能。 - 2 - 表5.1.1 绘图参数表 色彩字符颜色线型字符线型格式标记符号数据点形式标记符号数据点形式 y 黄- 实线. 点<小于号 m 紫：点线o 圆s 正方形 c 青-. 点划线x 叉号 d 菱形 r 红- - 虚线+ 加号h 六角星 g 绿* 星号p 五角星 b 蓝v 向下三角形 w 白^ 向上三角形 k 黑>大于号 例如，在上例中输入 >> plot(x,y1,'r+-',x,y2,'k*:') 则得图5.1.2 图5.1.2 使用不同标记的plot 函数绘制的正弦曲线 5.1.2 图形修饰 MATLAB 软件为用户提供了一些特殊的图形函数，用于修饰已经绘制好的图形。 表5.1.2 图形修饰函数表
matlab用于计算方法的源程序
1、Newdon迭代法求解非线性方程 function [x k t]=NewdonToEquation(f,df,x0,eps) %牛顿迭代法解线性方程 %[x k t]=NewdonToEquation(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:原函数，定义为内联函数 ?:函数的倒数，定义为内联函数 %x0:初始值 %eps:误差限 % %应用举例： %f=inline('x^3+4*x^2-10'); ?=inline('3*x^2+8*x'); %x=NewdonToEquation(f,df,1,0.5e-6) %[x k]=NewdonToEquation(f,df,1,0.5e-6) %[x k t]=NewdonToEquation(f,df,1,0.5e-6) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquation(f,df,1) if nargin==3 eps="0".5e-6; end tic; k=0; while 1 x="x0-f"(x0)./df(x0); k="k"+1; if abs(x-x0) < eps || k >30 break; end x0=x; end t=toc; if k >= 30 disp('迭代次数太多。'); x="0"; t="0"; end
2、Newdon迭代法求解非线性方程组 function y="NewdonF"(x) %牛顿迭代法解非线性方程组的测试函数 %定义是必须定义为列向量 y(1,1)=x(1).^2-10*x(1)+x(2).^2+8; y(2,1)=x(1).*x(2).^2+x(1)-10*x(2)+8; return; function y="NewdonDF"(x) %牛顿迭代法解非线性方程组的测试函数的导数 y(1,1)=2*x(1)-10; y(1,2)=2*x(2); y(2,1)=x(2).^+1; y(2,2)=2*x(1).*x(2)-10; return; 以上两个函数仅供下面程序的测试 function [x k t]=NewdonToEquations(f,df,x0,eps) %牛顿迭代法解非线性方程组 %[x k t]=NewdonToEquations(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:方程组（事先定义） ?:方程组的导数（事先定义） %x0:初始值 %eps:误差限 % %说明：由于虚参f和df的类型都是函数，使用前需要事先在当前目录下采用函数M文件定义% 另外在使用此函数求解非线性方程组时，需要在函数名前加符号“@”，如下所示 % %应用举例： %x0=[0,0];eps=0.5e-6; %x=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)
中国矿业大学 实验六 MATLAB数据可视化
实验六MATLAB数据可视化 一、实验目的 掌握MATLAB 二维、三维图形绘制，掌握图形属性的设置和图形修饰；掌握图像文件的读取和显示。 二、实验内容 (1) 二维图形绘制。 (2) 三维曲线和三维曲面绘制。 三、实验步骤 1．二维图形绘制 (1) 二维图形绘制主要使用函数plot。 >> clear all; >> x=linspace(0,2*pi,100); >> y1=sin(x); >> plot(x,y) >> hold on >> y2=cos(x) >> plot(x,y) >> hold off
注：hold on 用于保持图形窗口中原有的图形，hold off解除保持。 (2) 函数plot 的参数也可以是矩阵。 >> close all >> x=linspace(0,2*pi,100); >> y1=sin(x); >> y2=cos(x); >> A=[y1 ; y2]'; >> B=[x ; x]' >> plot(B,A)
(3) 选用绘图线形和颜色。>> close all >> plot(x,y1,'g+',x,y2, 'r:') >> grid on
(4) 添加文字标注。 >> title('正弦曲线和余弦曲线') >> ylabel('幅度') >> xlabel('时间') >> legend('sin(x)', 'cos(x)') >> gtext('\leftarrowsinx')
(5) 修改坐标轴范围。 >> axis equal >> axis normal >> axis([0 pi 0 1.5]) 程序如下： x=linspace(0,2*pi,100); y1=sin(x); y2=cos(x); A=[y1 ; y2]'; B=[x ; x]' plot(B,A) plot(x,y1,'g+',x,y2, 'r:') axis equal axis normal axis([0 pi 0 1.5])
0计算方法及MATLAB实现简明讲义课件PPS8-1欧拉龙格法
第8章 常微分方程初值问题数值解法 8.1 引言 8.2 欧拉方法 8.3 龙格-库塔方法 8.4 单步法的收敛性与稳定性 8.5 线性多步法
8.1 引 言 考虑一阶常微分方程的初值问题 00(,),[,],(). y f x y x a b y x y '=∈=（1.1） （1.2） 如果存在实数 ，使得 121212(,)(,).,R f x y f x y L y y y y -≤-?∈（1.3） 则称 关于 满足李普希茨(Lipschitz ）条件， 称为 的李普希茨常数（简称Lips.常数）. 0>L f y L f (参阅教材386页)
计算方法及MATLAB 实现 所谓数值解法，就是寻求解 在一系列离散节点 )(x y <<<<<+121n n x x x x 上的近似值 . ,,,,,121+n n y y y y 相邻两个节点的间距 称为步长. n n n x x h -=+1 如不特别说明，总是假定 为定数， ),2,1( ==i h h i 这时节点为 . ) ,2,1,0(0 =+=i nh x x n 初值问题(1.1),(1.2)的数值解法的基本特点是采取 “步进式”. 即求解过程顺着节点排列的次序一步一步地向前推进. 00(,),[,], (). y f x y x a b y x y '=∈=
描述这类算法，只要给出用已知信息 ,,,21--n n n y y y 计算 的递推公式. 1+n y 一类是计算 时只用到前一点的值 ，称为单步法. 1+n y n y 另一类是用到 前面 点的值 ， 1+n y k 11,,,+--k n n n y y y 称为 步法. k 其次，要研究公式的局部截断误差和阶，数值解 与 精确解 的误差估计及收敛性，还有递推公式的计算 稳定性等问题. n y )(n x y 首先对方程 离散化，建立求数值解的递推 公式. ),(y x f y ='
用MATLAB实现结构可靠度计算.
用MATLAB实现结构可靠度计算 口徐华…朝泽刚‘u刘勇‘21 。 (【l】中国地质大学(武汉工程学院湖北?武汉430074; 12】河海大学土木工程学院江苏?南京210098 摘要:Matlab提供了各种矩阵的运算和操作,其中包含结构可靠度计算中常用的各种数值计算方法工具箱,本文从基本原理和相关算例分析两方面,阐述利用Matlab,编制了计算结构可靠度Matlab程.序,使得Matlab-语言在可靠度计算中得到应用。 关键词:结构可靠度Matlab软件最优化法 中图分类号:TP39文献标识码:A文章编号:1007-3973(200902-095-Ol 1结构可靠度的计算方法 当川概率描述结构的可靠性时,计算结构可靠度就是计算结构在规定时问内、规定条件F结构能够完成预定功能的概率。 从简单到复杂或精确稃度的不同,先后提出的可靠度计算方法有一次二阶矩方法、二次二阶矩方法、蒙特卡洛方法以及其他方法。一次■阶矩方法又分为。I-心点法和验算点法,其中验算点法足H前可靠度分析最常川的方法。 2最优化方法计算可靠度指标数学模型 由结构111n个任意分布的独立随机变量一,x:…以表示的结构极限状态方程为:Z=g(■.托…t=0,采用R-F将非正念变量当罱正态化,得到等效正态分布的均值o:和标准差虹及可靠度指标B,由可靠度指标B的几何意义知。o;辟
开始时验算点未知,把6看成极限状态曲面上点P(■,爿:---37,的函数,通过优化求解,找到B最小值。求解可靠皮指标aJ以归结为以下约束优化模型: rain睁喜t华,2 s.,.Z=g(工i,x2’,…,工:=0 如极限状态方栉巾某个变最(X。可用其他变量表示,则上述模型jfIJ‘转化为无约束优化模型: 。。B!:手f生丛r+阻:坚:坠:盐尘}二剐 t∞oY?’【叫,J 3用MATLAB实现结构可靠度计算 3.1Matlab简介 Matlab是++种功能强、效率高、便.丁.进行科学和工程计算的交互式软件包,汇集了人量数学、统计、科学和工程所需的函数,MATI.AB具有编程简甲直观、用户界mf友善、开放性强等特点。将MATLAB用于蒙特卡罗法的一个显著优点是它拥有功能强大的随机数发生器指令。 3.2算例 3.2.I例:已知非线形极限状态方程z=g(t r'H=567f r-0.5H2=0’f、r服从正态分布。IIf=0.6,o r=0.0786;la|_ 2.18,o r_0.0654;H服从对数正态分布。u H= 3218,O。 =0.984。f、r、H相互独立,求可靠度指标B及验算点(,,r’,H‘。 解:先将H当量正念化:h=ln H服从正态分布,且 ,‘-““了:等专虿’=,。49?口二-、『五ir面_。。3
【原创】MATLAB实验报告-第二次-用MATLAB实现计算数据可视化-北京交通大学
MATLAB 上机实验报告( 2 ) 实验内容: 一、试用如下几种方法来建立向量，观察结果 ( 1) x=1:5, x=(1:5) ' 实验结果：x=1:5是行向量,x=(1:5)是列向量.且1为初始值，5为终止值，默认的步长为 1. >> x=1:5 1 2 3 4 5 >> x=(1:5)' x = 1 2
3 4 5 ( 2) x=0:pi/4:pi 实验结果：x=0:pi/4:pi 指的是x=(0,0.25*pi,0.50*pi,0.75*pi,pi). 其中pi为圆周率初始值为0,终止值为pi,步长为pi/4. >> x=0:pi/4:pi x = 0 0.7854 1.5708 2.3562 3.1416 (3)x=(0：0.2：3) ', y=e-x)p.(*sin(x) 实验结果：x的初始值为0,终止值为3,步长为0.2.而函数y表示将x向量中的每一个数代入函数y=e%x)*sin(x)得到的函数值组成的向量. >> x=(0:0.2:3)', y=exp(-x).*sin(x)
x = 0.2000 0.4000 0.6000 0.8000 1.0000 1.2000 1.4000 1.6000 1.8000 2.0000 2.2000 2.4000 2.6000 2.8000 3.0000
0.1627 0.2610 0.3099 0.3223 0.3096 0.2807 0.2430 0.2018 0.1610 0.1231 0.0896 0.0613 0.0383 0.0204 0.0070 (4) k=linspace(-pi,pi,5), k=logspace(-3,-1,5) 实验结果：k=linspace(-pi,pi,5)产生的是初始值为-pi，终止值为 pi,元素总数为5的行向量，即k的步长为pi/2. k=logspace(-3,-1,5产生的是初始值为10八(-3),终止值为10八(-1),元素总数为5的列向量.
计算方法上机实验报告-MATLAB
《计算方法》实验报告 指导教师： 学院： 班级： 团队成员：
一、题目 例2.7应用Newton 迭代法求方程210x x --=在1x =附近的数值解 k x ，并使其满足8110k k x x ---< 原理： 在方程()0f x =解的隔离区间[],a b 上选取合适的迭代初值0x ，过曲线()y f x =的点()() 00x f x ，引切线 ()()()1000:'l y f x f x x x =+- 其与x 轴相交于点：()() 0100 'f x x x f x =-，进一步,过曲线()y f x =的 点()()11x f x , 引切线 ()()()2111: 'l y f x f x x x =+- 其与x 轴相交于点：() () 1211 'f x x x f x =- 如此循环往复，可得一列逼近方程()0f x =精确解*x 的点 01k x x x ,,,,,其一般表达式为： ()() 111 'k k k k f x x x f x ---=- 该公式所表述的求解方法称为Newton 迭代法或切线法。
程序： function y=f(x)%定义原函数 y=x^3-x-1; end function y1=f1(x0)%求导函数在x0点的值 syms x; t=diff(f(x),x); y1=subs(t,x,x0); end function newton_iteration(x0,tol)%输入初始迭代点x0及精度tol x1=x0-f(x0)/f1(x0);k=1;%调用f函数和f1函数 while abs(x1-x0)>=tol x0=x1;x1=x0-f(x0)/f1(x0);k=k+1; end fprintf('满足精度要求的数值为x(%d)=%1.16g\n',k,x1); fprintf('迭代次数为k=%d\n',k); end 结果：
计算方法及其MATLAB实现第二章作业
作者：夏云木子 1、 >> syms re(x) re(y) re(z) >> input('计算相对误差：'),re(x)=10/1991,re(y)=0.0001/1.991,re(y)=0.0000001/0.0001991 所以可知re(y)最小，即y精度最高 2、 >> format short,A=sqrt(2) >> format short e,B=sqrt(2) >> format short g,C=sqrt(2)
>> format long,D=sqrt(2) >> format long e,E=sqrt(2) >> format long g,F=sqrt(2) >> format bank,H=sqrt(2) >> format hex,I=sqrt(2) >> format +,J=sqrt(2) >> format,K=sqrt(2)
3、 >> syms A >> A=[sqrt(3) exp(7);sin(5) log(4)];vpa(pi*A,6) 4、1/6251-1/6252=1/6251*6252 5、（1）1/（1+3x）-(1-x)/(1+x)=x*(3*x-1)/[(1+3*x)*(1+x)] (2) sqrt(x+1/x)-sqrt(x-1/x)=2/x/[sqrt(x-1/x)+sqrt(x+1/x)] (3) log10(x1)-log(x2)=log10(x1/x2) (4) [1-cos(2*x)]/x =x^2/factorial(2)-x^4/factorial(4)+x^6/factorial(6)-…
实验二MATLAB计算的可视化
课程实验报告 学年学期2011-2012学年第1学期 课程名称MATLAB与科学计算 实验名称实验二MATLAB计算的可视化实验室测量测绘实验中心计算机室专业年级热动113 学生姓名白治朋 学生学号2011012106 提交时间2013年10月23日 成绩 任课教师许景辉、牛亚斌 水利与建筑工程学院
实验二 MATLAB 计算的可视化 1、目的和要求 （1）熟练掌握MATLAB 二维曲线、三维图形的绘制。 （2）熟练掌握各种特殊图形的绘制。 （3）熟练掌握三维图形绘制命令。 （4）了解GUI 设计的一般过程和方法。 2、内容和步骤 参见教材实验四。 3、实验报告提交要求 （1） x=[1 2 3]，y=[1 2;2 3;5 8]，z=[2 6 9;3 8 8;1 5 7]，绘制plot （x ，y ）、plot （x ，z ），说明其各 自绘制的内容。 （2） 绘制如下图形，建立figure （2），绘图同样曲线，但标题为“你的姓名（黑体，16号字）”， 在x 坐标和y 坐标上分别标识学号和班级名称，并将网格线打开。 数组X 的列个数与矩阵y 的行个数相同， plot （ x ，y ）绘制的是x 为横坐标y 的每列为纵坐标的图像。如图1。 图1 数组X 的列个数与方阵z 的行列个数相同，plot （x ，z ）绘制的是x 为横坐标z 的每列为纵坐标的图像。如图2。 图2
（3）演示P133页，例题4.17 。
（4）完成课本P336图S 4.1实验，并用.m文件显示其程序内容。 （5）完成P302第四章例题4.
（6）通过绘制二阶系统阶跃响应，综合演示图形标识，请注释每条命令的含义。 clf; %清除图形窗口 t=6*pi*(0:100)/100;y=1-exp(-0.3*t).*cos(0.7*t); % 数据准备 tt=t(find(abs(y-1)>0.05)); %找出符合条件(y-1)的绝对值>0.05的对应t，赋值给tt ts=max(tt); %ts为tt中最大值ts=9.6133 plot(t,y,'r-','LineWidth',3) %画曲线t-y,红色实线，线粗3磅 axis([-inf,6*pi,0.6,inf]) %设置坐标轴范围。x轴下限自动产生,上限为6*pi；y轴下限0.6,上限自动产生 set(gca,'Xtick',[2*pi,4*pi,6*pi],'Ytick',[0.95,1,1.05,max(y)]) %二维坐标刻度设置。x轴刻度线取2*pi,4*pi,6*pi,y轴取0.95,1,1.05,max(y) grid on %显示坐标刻度线 title('\it y = 1 - e^{ -\alphat}cos{\omegat}') %用斜体1书写图名 text(13.5,1.2,'\fontsize{12}{\alpha}=0.3') %图形标识，添加文字注释。在x=13.5,y=1.2处,字体大小12磅，标注α=0.3 text(13.5,1.1,'\fontsize{12}{\omega}=0.7') %图形标识，添加文字注释。在x=13.5,y=1.1处,字体大小12磅，标注ω=0.7 hold on; %保持原有图形 plot(ts,0.95,'bo','MarkerSize',10); %在x=ts，y=0.95处画蓝色的空心圆圈,大小为
数值计算方法matlab程序
function [x0,k]=bisect1(fun1,a,b,ep) if nargin<4 ep=1e-5; end fa=feval(fun1,a); fb=feval(fun1,b); if fa*fb>0 x0=[fa,fb]; k=0; return; end k=1; while abs(b-a)/2>ep x=(a+b)/2; fx=feval(fun1,x); if fx*fa<0 b=x; fb=fx; else a=x; fa=fx;
end end x0=(a+b)/2; >> fun1=inline('x^3-x-1'); >> [x0,k]=bisect1(fun1,1.3,1.4,1e-4) x0 = 1.3247 k = 7 >> 简单迭代法 function [x0,k]=iterate1(fun1,x0,ep,N) if nargin<4 N=500; end if nargin<3 ep=1e-5; end x=x0; x0=x+2*ep;
while abs(x-x0)>ep & k
计算方法及其MATLAB实现第一章作业
计算方法作业（作者：夏云木子） 1、help linspace type linspace 2、a1=[5 12 47;13 41 2;9 6 71];a2=[12 9;6 15;7 21];B=a1*a2, C=a1(:,1:2).*a2, D=a1.^2,
E=a1(:).^2 3、a1=[5 12 47;13 41 2;9 6 71];a2=[12 9;6 15;7 21];a1(4:5,1:3)=a2.';a1([4 5],:)=a1([5 4],:);b1=a1
c1=b1(4,1),c2=b1(5,3),D=b1(3:4,:)*a2 4、a1=[5 12 47;13 41 2;9 6 71]; E=eye(3,3); S = a1 + 5*a1' - E, S1=a1^3-rot90(a1)^2+6*E 5、a1=[5 12 47;13 41 2;9 6 71];s=5;A=s-a1,B=s*a1,C=s.*a1,D=s./a1,E=a1./s
6、c=[1 2 3 4;5 6 7 8;9 10 11 12;13 14 15 16];A=c^-4,B=(c^3)^-1,C=(3*c+5*c^-1)/5
7、a=[1 i 3;9i 2-i 8;7 4 8+i];A=a.' 8、abc=[-2.57 8.87;-0.57 3.2-5.5i];m1=sign(abc),m2=round(abc),m3=floor(abc) Sign为符号函数，round表示四舍五入取整，floor表示舍去小数部分取整
9、x=[1 4 3 2 0 8 10 5]';y=[8 0 0 4 2 1 9 11]';A=dot(x,y) 10、a=[3.82 5.71 9.62];b=[7.31 6.42 2.48];A=dot(a,b),B=cross(a,b) 11、P=[5 7 8 0 1];Pf=poly(P);Px=poly2str(Pf,'x') 12、P=[3 0 9 60 0 -90];K1=polyval(P,45),K2=polyval(P,-123),K3=polyval(P,579) 13、P1=[13 55 0 -17 9];P2=[63 0 26 -85 0 105];PP=conv(P1,P2);P1P2=poly2str(PP,'x'),[Q,r]=deconv(P2,P1)
层次分析法计算权重在matlab中的实现
信息系统分析与设计作业 层次分析法确定绩效评价权重在matlab中的实现 小组成员：孙高茹、王靖、李春梅、郭荣1 程序简要概述 编写程序一步实现评价指标特征值lam、特征向量w以及一致性比率CR的求解。 具体的操作步骤是：首先构造评价指标，用专家评定法对指标两两打分，构建比较矩阵，继而运用编写程序实现层次分析法在MATLAB中的应用。 通过编写MATLAB程序一步实现问题求解，可以简化权重计算方法与步骤，减少工作量，从而提高人力资源管理中绩效考核的科学化电算化。 2 程序在matlab中实现的具体步骤 function [w,lam,CR] = ccfx(A) %A为成对比较矩阵，返回值w为近似特征向量 % lam为近似最大特征值λmax，CR为一致性比率 n=length(A(:,1)); a=sum(A); B=A %用B代替A做计算 for j=1:n %将A的列向量归一化 B(:,j)=B(:,j)./a(j); end s=B(:,1); for j=2:n s=s+B(:,j); end c=sum(s);%计算近似最大特征值λmax w=s./c; d=A*w lam=1/n*sum((d./w)); CI=(lam-n)/(n-1);%一致性指标 RI=[0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45,1.49,1.51];%RI为随机一致
性指标 CR=CI/RI(n);%求一致性比率 if CR>0.1 disp('没有通过一致性检验'); else disp('通过一致性检验'); end end 3 案例应用 我们拟构建公司员工绩效评价分析权重，完整操作步骤如下： 3.1构建的评价指标体系 我们将影响员工绩效评定的指标因素分为：打卡、业绩、创新、态度与品德。 3.2专家打分，构建两两比较矩阵 A = 1.0000 0.5000 3.0000 4.0000 2.0000 1.0000 5.0000 3.0000 0.3333 0.2000 1.0000 2.0000 0.2500 0.3333 0.5000 1.0000 3.3在MATLAB中运用编写好的程序实现 直接在MATLAB命令窗口中输入 [w,lam,CR]=ccfx(A) 继而直接得出 d = 1.3035 2.0000 0.5145 0.3926 w = 0.3102 0.4691 0.1242 0.0966 lam =4.1687
第6章MATLAB计算结果可视化讲解
第六章MATLAB 计算结果可视化 6.1连续函数和离散函数的可视化 【例6-1】用图形表示离散函数1 ) 6(--=n y 。 n=0:12; %产生一组自变量数据 y=1./abs(n-6); %计算相应点的函数值 plot(n,y,'r*','MarkerSize',20) %用红花标出数据点 grid on %画坐标方格 【例6-2】用图形表示连续调制波形)9sin()sin(t t y =。 t1=(0:11)/11*pi; y1=sin(t1).*sin(9*t1); t2=(0:100)/100*pi; y2=sin(t2).*sin(9*t2); subplot(2,2,1),plot(t1,y1,'r.'),axis([0,pi,-1,1]),title('子图 (1)') subplot(2,2,2),plot(t2,y2,'r.'),axis([0,pi,-1,1]),title('子图 (2)') subplot(2,2,3),plot(t1,y1,t1,y1,'r.') axis([0,pi,-1,1]),title('子图 (3)') subplot(2,2,4),plot(t2,y2)
6.2二维曲线绘图的基本操作 6.2.1 plot 的基本调用格式 【例6-3】用图形表示连续调制波形)9sin()sin(t t y 及其包络线。 t=(0:pi/100:pi)'; %长度为101的时间采样列向量 y1=sin(t)*[1,-1]; %包络线函数值，是（101x2）的矩阵 y2=sin(t).*sin(9*t); %长度为101的调制波列向量 t3=pi*(0:9)/9; y3=sin(t3).*sin(9*t3);plot(t,y1,'r:',t,y2,'b',t3,y3,'bo') 6.2.2泛函绘图指令fplot 【例6-4】fplot 与一般绘图指令的绘图效果比较。 [x,y]=fplot('cos(tan(pi*x))',[-0.4,1.4],0.2e-3);n=length(x); subplot(1,2,1),plot(x,y) title('\fontsize{20}\fontname{隶书}泛函绘图指令效果') t=(-0.4:1.8/n:1.4)'; subplot(1,2,2),plot(t,cos(tan(pi*t))) 6.2.3曲线的色彩、线型和数据点形 【例6-5】用图形演示平面上一个方块四个顶点在仿射投影（Affine Projection ）下的位置、 形状变化。 %平面上的四个点和它们构成的方块