伊辛模型MATLAB源代码

合集下载

Matlab源程序代码

Matlab源程序代码

Matlab源程序代码正弦波的源程序:(一),用到的函数1,f2t函数function x=f2t(X)global dt df t f T N%x=f2t(X)%x为时域的取样值矢量%X为x的傅氏变换%X与x长度相同并为2的整幂%本函数需要一个全局变量dt(时域取样间隔) X=[X(N/2+1:N),X(1:N/2)];x=ifft(X)/dt;end2,t2f函数。

function X=t2f(x)global dt df N t f T%X=t2f(x)%x为时域的取样值矢量%X为x的傅氏变换%X与x长度相同,并为2的整幂。

%本函数需要一个全局变量dt(时域取样间隔) H=fft(x);X=[H(N/2+1:N),H(1:N/2)]*dt;end(二),主程序。

1,%(1)绘出正弦信号波形及频谱global dt df t f Nclose allk=input('取样点数=2^k, k取10左右');if isempty(k), k=10; endf0=input('f0=取1(kz)左右');if isempty(f0), f0=1; endN=2^k;dt=0.01; %msdf=1/(N*dt); %KHzT=N*dt; %截短时间Bs=N*df/2; %系统带宽f=[-Bs+df/2:df:Bs]; %频域横坐标t=[-T/2+dt/2:dt:T/2]; %时域横坐标s=sin(2*pi*f0*t); %输入的正弦信号S=t2f(s); %S是s的傅氏变换a=f2t(S); %a是S的傅氏反变换a=real(a);as=abs(S);subplot(2,1,1) %输出的频谱plot(f,as,'b');gridaxis([-2*f0,+2*f0,min(as),max(as)]) xlabel('f (KHz)') ylabel('|S(f)| (V/KHz)') %figure(2) subplot(2,1,2) plot(t,a,'black') %输出信号波形画图gridaxis([-2/f0,+2/f0,-1.5,1.5])xlabel('t(ms)')ylabel('a(t)(V)')gtext('频谱图')最佳基带系统的源程序:(一),用到的函数f2t函数和t2f函数。

因子分析MATLAB程序源代码

因子分析MATLAB程序源代码

因子分析MATLAB程序源代码因子分析是一种统计分析方法,用于确定一组观测变量的潜在因子结构。

在MATLAB中,可以使用`factoran`函数进行因子分析。

下面是一个示例的MATLAB代码,用于执行因子分析,并使用经典的因子旋转方法-变简单结构法(varimax)进行因子旋转:```matlab%数据准备data = [1 2 3 4 5; 2 3 4 5 6; 3 4 5 6 7; 4 5 6 7 8; 5 6 7 8 9];%因子分析num_factors = 2; % 欲提取的因子个数method = 'varimax'; % 因子旋转方法[loadings, spec_var, T, stats] = factoran(data, num_factors, 'rotate', method);%输出结果disp('因子载荷矩阵:');disp(loadings);disp('特殊因子方差:');disp(spec_var);disp('变换矩阵:');disp(T);disp('统计信息:');disp(stats);```在这个示例中,我们首先准备了一个包含5个观测变量的数据矩阵`data`,然后指定欲提取的因子个数`num_factors`为2,并选择`varimax`作为因子旋转方法。

调用`factoran`函数进行因子分析后,会返回以下结果:- `loadings`:因子载荷矩阵,描述所有观测变量与因子之间的关系;- `spec_var`:特殊因子方差,每个观测变量独自解释的方差;-`T`:变换矩阵,将数据旋转到简单结构后的载荷矩阵;- `stats`:统计信息,包括共同度、特殊方差和解释总方差等。

使用`disp`函数将结果打印到控制台上。

以上就是一个简单的因子分析的MATLAB代码示例,你可以根据自己的需求修改和扩展代码。

Matlab第一章代码

Matlab第一章代码

阅读使人快乐,成长需要时间
江阴室内设计培训学校排名
一、江阴问鼎教育(澄江中路5号东都大厦8楼)
培训方式:小班制授课,10人以内,高端室内设计培训机构;
影响力:江阴高性价比的室内设计培训品牌;
特点:环境特优雅,教学质量上乘;
上课模式:小班化教学,不限学时,包学包会;
不足:广告做的不多,主要靠口碑推广
面向对象:成人、学生。

二、江阴逸仙教育
培训方式:面授为主,辅助在线预习等学习;
影响力:一家留学生开的机构,广告打的挺多的;
特点:课堂比较漂亮,环境比较好。

学习时间自由,要自制力比较好的学员;
不足:价格比较贵,要1万以上;课时数比较少;
上课模式:课堂面授一对4+大班课+在线学习;
面对对象:成人。

三、远大培训机构
影响力:成人培训的影响力不如上面几家培训做的那么闻名遐迩,不过品牌效用还是很好;特点:老师口才比较好,上课气氛比较好,价钱也还好;
不足:老生常谈,创新度没其他机构高,上课人数比较多;
师资:兼职老师;
上课模式:课堂面授、大班为主;
面对对象:大学生、中学生。

四、英华培训机构
培训方式:课堂面授;
影响力:打入江苏的培训机构,知名度还可以;
特点:全部学员使用3阶段学习法,自学可以借鉴这个方法。

学习时间自由,要人监督才行;不足:课程设置比较制式化,主题比较少,然后也比较贵,1万以上了;
大概描述:提供了三阶段的学习方法;
师资:中教、外教;
上课模式:课堂面授、一对4;
面对对象:企业和个人。

pearson建模实例matlab代码

pearson建模实例matlab代码

在数据分析领域中,Pearson相关系数是一种用来衡量两个变量之间线性相关程度的统计量。

它的取值范围在-1到1之间,0表示没有线性相关,-1表示完全负相关,1表示完全正相关。

Pearson相关系数被广泛应用于数据分析、机器学习和统计学中,对于研究变量之间的相关性、趋势和关联性都有着重要的作用。

在本文中,我们将以Pearson建模实例为主题,结合Matlab代码进行深入探讨。

通过示例代码的分析和讨论,旨在为读者提供对Pearson建模及其在Matlab中的应用有全面、深刻的理解和应用能力。

1. Pearson相关系数让我们简要回顾一下Pearson相关系数的计算公式:\[ r = \frac{n(\sum{xy}) - (\sum{x})(\sum{y})}{\sqrt{(n\sum{x^2} - (\sum{x})^2)(n\sum{y^2} - (\sum{y})^2)}} \]在这个公式中, \( \sum \) 代表总和, \( \sum{xy} \) 表示 x 和 y 变量对应数据的乘积之和, \( \sum{x} \) 和 \( \sum{y} \) 分别表示 x 和y 变量的数据之和, \( \sum{x^2} \) 和 \( \sum{y^2} \) 分别表示 x 和 y 变量数据的平方和,n 代表样本数量。

2. Pearson建模实例现在,我们将通过一个具体的实例来说明如何使用Matlab进行Pearson建模的实践。

假设我们有两个变量 x 和 y,我们想要计算它们之间的Pearson相关系数,并用Matlab代码实现。

在Matlab中,我们可以使用 corrcoef 函数来计算两个变量之间的Pearson相关系数。

以下是一个简单的示例代码:```Matlab% 假设我们有两个变量 x 和 yx = [1, 2, 3, 4, 5];y = [2, 4, 6, 8, 10];% 使用 corrcoef 函数计算Pearson相关系数r = corrcoef(x, y);% 显示计算结果disp('Pearson相关系数:');disp(r(1, 2));```在这段示例代码中,我们首先定义了两个变量 x 和 y,然后使用corrcoef 函数计算它们之间的Pearson相关系数,并最后输出计算结果。

人工神经网络matlab源程序代码_光环大数据人工智能培训

人工神经网络matlab源程序代码_光环大数据人工智能培训

人工神经网络matlab源程序代码_光环大数据人工智能培训人工神经网络matlab 源程序代码%产生指定类别的样本点,并在图中绘出X = [0 1; 0 1]; % 限制类中心的范围clusters = 5; % 指定类别数目points = 10; % 指定每一类的点的数目std_dev = 0.05; % 每一类的标准差P = nngenc(X,clusters,points,std_dev);plot(P(1,:),P(2,:),'+r');title('输入样本向量');xlabel('p(1)');ylabel('p(2)');%建立网络net=newc([0 1;0 1],5,0.1); %设置神经元数目为5%得到网络权值,并在图上绘出figure;plot(P(1,:),P(2,:),'+r');w=net.iw{1}hold on;plot(w(:,1),w(:,2),'ob');hold off;title('输入样本向量及初始权值');xlabel('p(1)');ylabel('p(2)');figure;plot(P(1,:),P(2,:),'+r');hold on;%训练网络net.trainParam.epochs=7;net=init(net);net=train(net,P);%得到训练后的网络权值,并在图上绘出w=net.iw{1}plot(w(:,1),w(:,2),'ob');hold off;title('输入样本向量及更新后的权值');xlabel('p(1)');ylabel('p(2)');a=0;p = [0.6 ;0.8];a=sim(net,p)-------------------example8_2%随机生成1000个二维向量,作为样本,并绘出其分布P = rands(2,1000);plot(P(1,:),P(2,:),'+r')title('初始随机样本点分布');xlabel('P(1)');ylabel('P(2)');%建立网络,得到初始权值net=newsom([0 1; 0 1],[5 6]);w1_init=net.iw{1,1}%绘出初始权值分布图figure;plotsom(w1_init,yers{1}.distances)%分别对不同的步长,训练网络,绘出相应的权值分布图for i=10:30:100net.trainParam.epochs=i;net=train(net,P);figure;plotsom(net.iw{1,1},yers{1}.distances)end%对于训练好的网络,选择特定的输入向量,得到网络的输出结果p=[0.5;0.3];a=0;a = sim(net,p)--------------------------example8_3%指定输入二维向量及其类别P = [-3 -2 -2 0 0 0 0 +2 +2 +3;0 +1 -1 +2 +1 -1 -2 +1 -1 0];C = [1 1 1 2 2 2 2 1 1 1];%将这些类别转换成学习向量量化网络使用的目标向量T = ind2vec(C)%用不同的颜色,绘出这些输入向量plotvec(P,C),title('输入二维向量');xlabel('P(1)');ylabel('P(2)');%建立网络net = newlvq(minmax(P),4,[.6 .4],0.1);%在同一幅图上绘出输入向量及初始权重向量figure;plotvec(P,C)hold onW1=net.iw{1};plot(W1(1,1),W1(1,2),'ow')title('输入以及权重向量');xlabel('P(1), W(1)');ylabel('P(2), W(2)');hold off;%训练网络,并再次绘出权重向量figure;plotvec(P,C);hold on;net.trainParam.epochs=150;net.trainParam.show=Inf;net=train(net,P,T);plotvec(net.iw{1}',vec2ind(net.lw{2}),'o');%对于一个特定的点,得到网络的输出p = [0.8; 0.3];a = vec2ind(sim(net,p))为什么大家选择光环大数据!大数据培训、人工智能培训、Python培训、大数据培训机构、大数据培训班、数据分析培训、大数据可视化培训,就选光环大数据!光环大数据,聘请专业的大数据领域知名讲师,确保教学的整体质量与教学水准。

matlab ann源代码

matlab ann源代码

matlab ann源代码如何编写一个简单的人工神经网络(ANN)的MATLAB源代码。

人工神经网络(Artificial Neural Networks,简称ANN)是一种模拟生物神经网络行为和结构的计算模型,它是一种广泛应用于机器学习和模式识别领域的人工智能技术。

通过对大量数据进行训练和学习,ANN 可以通过多层次的神经元之间的连接,进行复杂的非线性函数拟合和任务解决。

在本文中,我们将以MATLAB为编程语言,一步一步地介绍如何编写简单的人工神经网络的源代码。

下面是源代码的主要步骤:第一步:导入数据集在使用ANN进行训练和测试之前,我们需要准备一个数据集。

数据集应该包含输入特征以及对应的目标输出。

在MATLAB中,我们可以使用两个矩阵来表示输入和输出数据,其中每一行表示一个样本,每一列表示一个特征。

导入数据集的代码如下:matlab导入数据集load iris_dataset.mat第二步:设置网络结构在设计ANN之前,我们需要设置网络的结构,包括输入层、隐藏层和输出层的神经元数量。

一般来说,输入层的神经元数量应该与数据集的特征数量相同,而输出层的神经元数量应该与问题的类别数量相同。

隐藏层的数量和神经元数量可以根据具体问题的复杂性进行调整。

设置网络结构的代码如下:matlab设置网络结构input_size = size(inputs,2);output_size = size(targets,2);hidden_size = 10; 设置隐藏层神经元数量第三步:初始化权重和偏置在ANN中,神经元之间的连接强度由权重和偏置决定。

我们需要为所有神经元之间的连接,以及每个神经元的偏置,随机初始化一个初始值。

初始化权重和偏置的代码如下:matlab初始化权重和偏置W1 = rand(input_size,hidden_size); 初始化输入层到隐藏层的权重b1 = rand(1,hidden_size); 初始化隐藏层的偏置W2 = rand(hidden_size,output_size); 初始化隐藏层到输出层的权重b2 = rand(1,output_size); 初始化输出层的偏置第四步:定义激活函数激活函数是ANN中非线性转换的关键。

计算材料学_Ising模型实验报告

计算材料学_Ising模型实验报告

Monte Carlo实验报告一、项目名称:Ising 模型二、项目内容概要1、编译和运行进入实验的文件夹:cd□~/sourcecode/2D_Ising文件夹里有源代码mc2d.f和输入文件in.2d阅读理解并编辑输入文件:gedit□in.2d之后编译mc2d.ff95 mc2d.f -o mc2d.exe运行可执行文件./mc2d.exe查看刚刚生成的四个输出文件,四个文件的内容如下:file1.out:温度;时间;单位原子能量;单位原子磁化强度file2.out:温度;单位原子能量;能量变化;单位原子磁化强度;磁化强度变化;单位原子热容file3.out:温度;自旋构型file4.out:温度;能量升高而被接受的数目;能量下降而被接受的数目;被拒绝的数目2、gnuplot 作图作温度与能量图:p “file2.out”u 1:2 w p ps 3 pt 5 作出file2.out 中第1 列与第2 列数据;作温度与磁化强度图:p “file2.out”u 1:4 w p ps 3 pt 5 作出file2.out 中第1 列与第4 列数据作温度与热容图:p “file2.out”u 1:6 w p ps 3 pt 5 作出file2.out 中第1 列与第6 列数据三、项目实施方法/原理1925 年,伊辛提出描写铁磁体的简化模型:设有N 个自旋组成的d 维晶格(d=1,2,3),第i 格点自旋为Si=±1(i=1,2,…N; ±代表上下)。

只考虑最近邻作用,相互作用能为±J(J>0 为铁磁性, J<0 为反铁磁性),平行为-J,反平行为J。

伊辛模型的蒙特卡洛模拟基本步骤如下:四、项目实施结果:1.各种情况下能量温度曲线能量温度 能量能量温度/K能量铁磁正方形点阵温度和能量曲线 铁磁三角形点阵能量与温度曲线能量温度/K能量温度/K反铁磁性正方形点阵能量温度曲线 反铁磁性正方形点阵外场为1时能量温度曲线能量能量/K势能温度/K反铁磁性正方形点阵外场为0.5时能量温度曲线2.各种情况下磁化强度和温度的关系曲线磁化强度温度/K磁化强度温度/K铁磁正方形点阵磁化强度能量曲线 铁磁三角形点阵磁化强度温度曲线磁化强度温度/K磁化强度温度/K磁化强度-温度反铁磁性正方形点阵磁化强度温度曲线 反铁磁性正方形点阵磁化强度温度曲线(外场为0.5)磁化强度温度/K磁化强度温度/K反铁磁性正方形点阵磁化强度温度曲线(外场为1) 铁磁性正方形点阵磁化强度温度曲线(外场为0.5)磁化强度温度/K铁磁性正方形点阵磁化强度温度曲线(外场为0.5)4.各种情况下热容和温度的关系图热容温度/K热容温度/K铁磁正方形点阵热容能量曲线铁磁三角形点阵热容能量曲线能量温度热容温度/K反铁磁正方形点阵热容能量曲线反铁磁正方形点阵热容能量曲线(外场为1)热容温度/K热容温度/K反铁磁正方形点阵热容能量曲线(外场为0.5) 铁磁正方形点阵热容能量曲线(外场为0.5)热容温度/K铁磁正方形点阵热容能量曲线(外场为1)五、项目小结:1.在保持原参数不变的情况下,可以得出,温度越高,原子热运动越剧烈,因此单个原子的能量也就越高。

GM(1,1)模型的MATLAB源代码2

GM(1,1)模型的MATLAB源代码2
x1=zeros(size(x0));
for k=1:n
for i=1:k
x1(k)=x1(k)+x0(i); %计算累加值,并将值赋予矩阵be
end
end
z=zeros(size(x0));
for k=2:n
z(k)=0.5*(x1(k)+x1(k-1)); %计算数据矩阵B的第一列数据
%px0为预测数列,rel为平均相对误差,rel为平均相对误差(为百分比)
%默认的number参数为原数组大小
if nargin==1 %对输入矩阵进行判断,如不是一维列矩阵,进行转置变换
number=max(size(x0));
end
n=max(size(x0)); %取输入数据的样本量
z(k)=0.5*(x1(k)+x1(k-1));
end
x0;
z;
B=[-(z(2:n))' (z(2:n).^2)'];
B;
Y=(x0(2:n))';
Y;
ab=inv(B'*B)*B'*Y;
a=ab(1);b=ab(2);
for k=1:number
px0(k)=(a*x1(1))/(b*x1(1)+(a-b*x1(1)).*exp(a*(k-1));1);
end
if nargin==1
number=max(size(x0));
end
n=max(size(x0)); %数组大小..
[px0,ab,rel]=gm11(x0,number);
wucha=x0-px0(1:n);

(完整word版)matlab经典代码大全

(完整word版)matlab经典代码大全

(完整word版)matlab经典代码大全哈哈哈MATLAB显示正炫余炫图:plot(x,y1,'* r',x,y2,'o b')定义【0,2π】;t=0:pi/10:2*pi;定义函数文件:function [返回变量列表]=函数名(输入变量列表)顺序结构:选择结构1)if-else-end语句其格式为:if 逻辑表达式程序模块1;else程序模块2;End图片读取:%选择图片路径[filename, pathname] = ...uigetfile({'*.jpg';'*.bmp';'*.gif'},'选择图片');%合成路径+文件名str=[pathname,filename];%为什么pathname和filename要前面出现的位置相反才能运行呢%读取图片im=imread(str);%使用图片axes(handles.axes1);%显示图片imshow(im);边缘检测:global imstr=get(hObject,'string');axes (handles.axes1);switch strcase ' 原图'imshow(im);case 'sobel'BW = edge(rgb2gray(im),'sobel');imshow(BW);case 'prewitt'BW = edge(rgb2gray(im),'prewitt');imshow(BW);case 'canny'BW = edge(rgb2gray(im),'canny');imshow(BW);Canny算子边缘定位精确性和抗噪声能力效果较好,是一个折中方案end;开闭运算:se=[1,1,1;1,1,1;1,1,1;1,1,1]; %Structuring ElementI=rgb2gray(im);imshow(I,[]);title('Original Image');I=double(I);[im_height,im_width]=size(I);[se_height,se_width]=size(se);halfheight=floor(se_height/2);halfwidth=floor(se_width/2);[se_origin]=floor((size(se)+1)/2);image_dilation=padarray(I,se_origin,0,'both'); %Image to be used for dilationimage_erosion=padarray(I,se_origin,256,'both'); %Image to be used for erosion %%%%%%%%%%%%%%%%%% %%% Dilation %%%%%%%%%%%%%%%%%%%%%for k=se_origin(1)+1:im_height+se_origin(1)for kk=se_origin(2)+1:im_width+se_origin(2)dilated_image(k-se_origin(1),kk-se_origin(2))=max(max(se+image_dilation(k-se_origin(1):k+halfh eight-1,kk-se_origin(2):kk+halfwidth-1)));endendfigure;imshow(dilated_image,[]);title('Image after Dilation'); %%%%%%%%%%%%%%%%%%%% Erosion %%%%%%%%%%%%%%%%%%%%se=se';for k=se_origin(2)+1:im_height+se_origin(2)for kk=se_origin(1)+1:im_width+se_origin(1)eroded_image(k-se_origin(2),kk-se_origin(1))=min(min(image_erosion(k-se_origin(2):k+halfwidth -1,kk-se_origin(1):kk+halfheight-1)-se));endendfigure;imshow(eroded_image,[]);title('Image after Erosion'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% Opening(Erosion first, then Dilation) %%% %%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%se=se';image_dilation2=eroded_image; %Image to be used for dilationfor k=se_origin(1)+1:im_height-se_origin(1)for kk=se_origin(2)+1:im_width-se_origin(2)opening_image(k-se_origin(1),kk-se_origin(2))=max(max(se+image_dilation2(k-se_origin(1):k+hal fheight-1,kk-se_origin(2):kk+halfwidth-1)));endendfigure;imshow(opening_image,[]);title('OpeningImage'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%% Closing(Dilation first, then Erosion) %%% %%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%se=se';image_erosion2=dilated_image; %Image to be used for erosionfor k=se_origin(2)+1:im_height-se_origin(2)for kk=se_origin(1)+1:im_width-se_origin(1)closing_image(k-se_origin(2),kk-se_origin(1))=min(min(image_erosion2(k-se_origin(2):k+halfwidt h-1,kk-se_origin(1):kk+halfheight-1)-se));endendfigure;imshow(closing_image,[]);title('Closing Image');Warning: Image is too big to fit on screen; displaying at 31% scale.> In truesize>Resize1 at 308In truesize at 44In imshow at 161图像的直方图归一化:I=imread(‘red.bmp’);%读入图像figure;%打开新窗口[M,N]=size(I);%计算图像大小[counts,x]=imhist(I,32);%计算有32个小区间的灰度直方图counts=counts/M/N;%计算归一化灰度直方图各区间的值stem(x,counts);%绘制归一化直方图图像平移:I=imread('shuichi.jpg');se=translate(strel(1),[180 190]);B=imdilate(I,se);figure;subplot(1,2,1),subimage(I);title('原图像');subplot(1,2,2),subimage(B);title('平移后图像');图像的转置;A=imread('nir.bmp');tform=maketform('affine',[0 1 0;1 0 0;0 0 1]);B=imtransform(A,tform,'nearest');figure;imshow(A);figure;imshow(B);imwrite(B,'nir转置后图像.bmp');图像滤波:B = imfilter(A,H,option1,option2,...)或写作g = imfilter(f, w, filtering_mode, boundary_options, size_options)其中,f为输入图像,w为滤波掩模,g为滤波后图像。

matlab实验代码(总)

matlab实验代码(总)

matlab实验代码(总)% 使⽤两种⽅法,创建⼀稀疏矩阵% 使⽤函数sparse,可以⽤⼀组⾮零元素直接创建⼀个稀疏矩阵。

该函数调⽤格式为:% S=sparse(i,j,s,m,n)% 其中i和j都为⽮量,分别是指矩阵中⾮零元素的⾏号与列号,% s是⼀个全部为⾮零元素⽮量,元素在矩阵中排列的位置为(i,j)% m为输出的稀疏矩阵的⾏数,n为输出的稀疏矩阵的列数。

%⽅法1A9=[0 0 1;0 3 0;2 4 0]B9=sparse(A9)C9=full(B9)%⽅法2A10=sparse([1 3 2 4],[2 3 1 4],[1 2 3 4],4,4)C10=full(A10)A11=[1 2 3];B11=[4 5 6];C11=3.^A11D11=A11.^B11%使⽤函数,实现矩阵左旋90°或右旋90°的功能。

A=[ 1 2 3 ; 4 5 6 ; 7 8 9 ]B=rot90(A,1)C=rot90(A,-1)%求S=2^0+2^1+2^2+2^3+2^4+……+2^10的值(提⽰:利⽤求和函数与累乘积函数。

)A=2*ones(1,10)%10个2B=cumprod(A)%平⽅C=sum(B)+1%加上2^0%建⽴⼀个字符串向量,删除其中的⼤写字母(提⽰:利⽤find函数和空矩阵。

)str='AAAbCcd'b=find(str>='A' & str<='Z');str(b)=[];% 输⼊⼀个百分制成绩,要求输出成绩等级A、B、C、D、E。

其中90分~100分为A,80分~89分为B,70分~79为C,60分~69分为D,60分以下为E。

switch(score)case num2cell(90:0.5:100)disp(['成绩等级为:A']);case num2cell(80:0.5:89.5)disp(['成绩等级为:B']);case num2cell(70:0.5:79.5)disp(['成绩等级为:C']);case num2cell(60:0.5:69.5)disp(['成绩等级为:D']);case num2cell(0:0.5:59.5)disp(['成绩等级为:E']);otherwisedisp(['输⼊成绩不合理!']);end%设计程序,完成两位数的加、减、乘、除四则运算,%即产⽣两个两位随机整数,再输⼊⼀个运算符号,做相应的运算,显⽰相应的结果,并要求结果显⽰类似于“a=x+y=34”。

matlab设计dmc模型代码

matlab设计dmc模型代码

标题:MATLAB设计DMC模型代码摘要:本文将详细介绍如何使用MATLAB软件编写DMC(Dynamic Matrix Control)模型代码。

DMC是一种基于模型的控制方法,可以有效地应用于工业过程的控制和优化中。

通过本文的学习,读者将能够掌握DMC模型的基本原理和编写代码的方法,为工业控制领域的应用提供参考。

正文:一、DMC模型概述DMC是一种基于模型的预测控制方法,它通过建立过程的数学模型,预测未来的系统响应,并根据预测结果进行控制。

DMC模型具有良好的鲁棒性和适应性,适用于各种工业过程的控制。

二、DMC模型代码编写步骤1. 数据采集和建模需对待控制的工业过程进行数据采集,并利用MATLAB中的系统辨识工具箱进行建模。

建立过程的动态数学模型是DMC模型设计的基础。

2. DMC模型参数确定根据建立的数学模型,确定DMC模型的参数,包括控制时域长度、预测时域长度、控制权重矩阵、预测权重矩阵等。

3. 编写MATLAB代码使用MATLAB软件,在命令窗口或脚本文件中编写DMC模型的代码。

主要包括以下步骤:(1)导入过程的数学模型(2)初始化模型参数(3)设定控制目标(4)进行预测控制(5)实时更新控制器输出三、DMC模型代码示例以下是一个简单的DMC模型代码示例,用于控制一个二阶惯性过程:```MATLAB导入过程模型G = tf([1],[1,1,0]);初始化模型参数N = 10;Nu = 3;lambda = 1;D = 100;设定控制目标r = ones(D,1);预测控制for k=1:Dif k<=10deltau(k)=1;elsedeltau(k)=0;endendfor k=1:Dy(k) = G * deltau(k);if k<=100u(k) = u(k-1)+deltau(k); elseu(k) = u(k-1);endend```四、DMC模型代码应用实例以一个温度控制系统为例,通过MATLAB编写DMC模型代码,对温度进行控制。

因子分析MATLAB程序源代码

因子分析MATLAB程序源代码

因⼦分析MATLAB程序源代码clear all;DATA=load('D:0.m');DATA=double(DATA);DATA=DATA';TESTDATA=load('D:14f.m');TESTDATA=double(TESTDATA);% DATA=load('D:正常.txt');% DATA=double(DATA);% DATA=DATA(:,3:12);% TESTDATA=load('D:异常.txt');% TESTDATA=double(TESTDATA);% TESTDATA=TESTDATA(:,3:12);[Kp,T2]=tztq(DATA,TESTDATA);function [contribution,T2,SPE,t2cl,s_cl] = PCA_model(Xtrain,Xtest)X_mean = mean(Xtrain);X_std = std(Xtrain);[X_row ,X_col]= size(Xtrain);for i = 1:X_colXtrain(:,i) = (Xtrain(:,i)-X_mean(i))./X_std(i);Xtest(:,i) = (Xtest(:,i)-X_mean(i))./X_std(i);end[U,S,V]=svd(Xtrain./sqrt(size(Xtrain,1)-1),0);D= S^2;lamda=diag(D);num_pc=1;while sum(lamda(1:num_pc))/sum(lamda)<0.9num_pc=num_pc+1;endD=diag(lamda);P=V(:,1:num_pc);[a,b]=size(Xtest);[r,y]=size(P*P');I=eye(r,y);e=Xtest*(I-P*P');for i=1:aT2(i)=Xtest(i,:)*P*inv(D(1:num_pc,1:num_pc))*P'*Xtest(i,:)';endfor l=1:aSPE(l)=e(l,:)*e(l,:)';endfor j=1:bcontribution(j)=(norm(e(:,j)))^2;endt2cl=num_pc*(X_row-1)*(X_row+1)*icdf('f',0.99,num_pc,X_row-num_pc)/(X_row*(X_row-num_pc)); for i=1:3theta(i)=trace((D(num_pc+1:X_col,num_pc+1:X_col))^i);end% 另⼀种SPE控制线算法% h=(theta(1)^2)/theta(2);% g=theta(2)/theta(1);% conf=0.95;% df=round(h);% delta2a1=g*pinv(df,2);h0=1-2*theta(1)*theta(3)/(3*theta(2)^2);ca=icdf('norm',0.99,0,1);s_cl=theta(1)*(ca*sqrt(2*theta(2)*h0^2)/theta(1)+1+theta(2)*h0*(h0-1)/theta(1)^2)^(1/h0); clear all;X1=load('D:0.m');Xtrain=X1';Xtrain=double(Xtrain);X2=load('D:14f.m');Xtest=X2;Xtest=double(Xtest);% X1=load('D:正常br.txt');% Xtrain=X1(:,3:62);% Xtrain=double(Xtrain);% X2=load('D:异常br.txt');% Xtest=X2(:,3:62);% Xtest=double(Xtest);[contribution,T2,SPE,t2cl,s_cl]=PCA_model(Xtrain,Xtest);figurex=size(Xtest,1);plot(1:x,SPE,'k');hold on;plot(1:x,s_cl,'r-');title('SPE');hold off;figureplot(1:x,T2,'K');title('T2');hold on;plot(1:x,t2cl,'r-');hold off;figurebar(contribution,'group')title('贡献图');function [Kp,T2]=tztq(ax,ay)[Nx]=size(ax);mean_X = mean(ax);axb=ax;std_X=std(ax);ax=ax-mean_X(ones(Nx,1),:);std_X(find(std_X==0))=1;%数据预处理ax=ax./std_X(ones(Nx,1),:);c=10000;% gama=0.05;% ni=1;% F1=ax(1,:);% F=F1';% for ib=2:Nx% for i=1:ni% for j=1:ni%% batar1(ib).block(i,j)=exp(-norm(ax(i,:)-ax(j,:))^2/c);% end% batar2(i,1)=exp(-norm(ax(i,:)-ax(ib,:))^2/c);% batar3(1,i)=exp(-norm(ax(ib,:)-ax(i,:))^2/c);% end% s1=exp(-norm(ax(ib,:)-ax(ib,:))^2/c);% batar= batar3(1,i)*inv(batar1(ib).block)* batar2(i,1);% s=(s1- batar)/s1;% if s>sqrt(gama)% ni=ni+1;% F=[F ax(ib,:)'];% end%%% end% AX=F'%训练集基的提取结束[N]=size(ax,1);for i=1:Nfor j=1:NK(i,j)=exp(-norm(ax(i,:)-ax(j,:))^2/c);%求核矩阵endendn1=ones(N,N);N1=1/N*n1;Kp=K-N1*K-K*N1+N1*K*N1;[u,s,v]=svd(Kp/N);lamda=s;P=v;lamda=diag(lamda);B=length(find(lamda>1e-10)); %求⾮零的特征值个数%求主元个数npc=1;while sum(lamda(1:npc))/sum(lamda(1:B))<0.9npc=npc+1;endnpc%求特征空间有效维数DimFS=1;while sum(lamda(1:DimFS))/sum(lamda(1:B))<=0.99 DimFS=DimFS+1;endlamda=diag(lamda);for i=1:B% P(:,i)=P(:,i)/norm(P(:,i)*s(i,i));P(:,i)=P(:,i)/(norm(P(:,i))*sqrt(N*lamda(i,i)));end[Ny]=size(ay,1);mean_X =mean(axb);std_X = std(axb);[num_sample] = Ny;ay = ay-mean_X(ones(num_sample,1),:);ay = ay./std_X(ones(num_sample,1),:);% mean_y = mean(ay);% std_y=std(ay);% ay = ay-mean_y(ones(Ny,1),:);% std_y(find(std_y==0))=1;%数据处理% ay = ay./std_y(ones(Ny,1),:);for i=1:Nyfor j=1:NKy(i,j)=exp(-norm(ay(i,:)-ax(j,:))^2/c);endendt1=ones(1,N);t11=1/N*t1;for i=1:Nykp1(i,:)= Ky(i,:)-t11*K- Ky(i,:)*N1+t11*K*N1;endfor i=1:Nyfor k=1:Bt(i,k)=P(:,k)'*kp1(i,:)';endend% 求T2,SPE% covtyb=inv(t'*t);for i=1:NyT2(i)=t(i,1:npc)*inv(lamda(1:npc,1:npc))*t(i,1:npc)'; %也可以% SPE(i)=t(i,1:npc)*t(i,1:npc)';% T2(1,i)=t(i,1:npc)*(covtyb(1:npc,1:npc))*t(i,1:npc)';SPE(i)=t(i,(npc+1):B)*t(i,(npc+1):B)';end%T2,SPE控制线t2cl=npc*(N-1)*(N+1)*icdf('f',0.99,npc,N-npc)/(N*(N-npc));for i=1:3theta(i)=trace((lamda(npc+1:DimFS,npc+1:DimFS))^i);endh0=1-2*theta(1)*theta(3)/(3*theta(2)^2);ca=icdf('norm',0.99,0,1);s_cl=theta(1)*(ca*sqrt(2*theta(2)*h0^2)/theta(1)+1+theta(2)*h0*(h0-1)/theta(1)^2)^(1/h0); figureplot(1:Ny,T2,'k');hold on;plot(1:Ny,t2cl,'r');title('T2');hold off;figureplot(1:Ny,SPE,'k')hold on;plot(1:Ny, s_cl,'r');title('SPE');hold off;。

GM(1-1)模型中的MATLAB程序

GM(1-1)模型中的MATLAB程序

GM (1,1)模型中的MATLAB 程序一、GM (1,1)模型的建立:(1)、一次累加生成序列的MA TLAB 命令:>> X0=[142 340 200 500 900 800 490 980 463 1100];>> X1(1)=X0(1)X1 =142>> for k=2:10X1(k)=X1(k-1)+X0(k)endX1 =142 482 682 1182 20822882 3372 4352 4815 5915(2)、由一次累加生成序列紧邻均值生成的)1(Z 的MA TLAB 命令:>> X0=[142 340 200 500 900 800 490 980 463 1100];>> X1(1)=X0(1);>> for k=2:10X1(k)=X1(k-1)+X0(k)Z(k)=(1/2)*(X1(k)+X1(k-1))EndX1 =142 482 682 1182 20822882 3372 4352 4815 5915Z =1.0e+003 *0 0.3120 0.5820 0.9320 1.6320 2.4820 3.12703.86204.58355.3650(3)、GM (1,1)的灰微分方程模型为:b k aZ k X =+)()()1()0(。

设∧α为待估计参数向量,⎥⎦⎤⎢⎣⎡=∧b a α。

利用最小二乘法得到Y B B B ')'(1-∧=α,MA TLAB 程序如下:>> B=[-Z(2:10)',ones(9,1)];>> Y=(X0(2:10))';>> alfa=inv(B'*B)*B'*Yalfa =-0.1062371.6018(4)、GM (1,1)的灰微分方程模型 b k aZ k X =+)()()1()0(的时间相应序列为:ab e a b X k X ak +⋅-=+-∧))1(()1()0()1( 由.6018.371,1062.0=-=b a 令.)1(,)0(u X v a b u -== 计算得到1.3499-=u , 1.3641=v 。

模型的matlab代码

模型的matlab代码

malmquist‐luenberger 模型的matlab代码Malmquist-Luenberger 模型是一个用于评估生产率变化的模型,它结合了Malmquist 指数和Luenberger 生产率指标。

以下是一个简单的MATLAB 代码示例,用于计算Malmquist-Luenberger 指数:matlab复制代码function[ml_index, g_index] =malmquist_luenberger(y_t, y_t1, y_t2)% 计算 Malmquist-Luenberger 指数f_t = y_t./y_t1; % 当前时期的技术前沿f_t1 = y_t1./y_t2; % 前一时期的技术前沿d_y = y_t - y_t1; % 当前时期与前一时期的生产距离函数d_y1 = y_t1 - y_t2; % 前一时期与前两时期的生产距离函数g_index = (f_t1.^2).*d_y1 ./ (f_t.^2).*d_y; % Luenberger 生产率指标ml_index = g_index ./ (1 + g_index); % Malmquist-Luenberger 指数end其中,y_t、y_t1和y_t2分别表示当前时期、前一时期和前两时期的产出向量。

该函数返回Malmquist-Luenberger 指数ml_index和Luenberger 生产率指标g_index。

需要注意的是,这只是一个简单的示例代码,实际应用中需要根据具体问题进行调整和改进。

另外,还需要注意数据的准确性和完整性,以及模型的适用范围和限制条件。

Malmquist-Luenberger 模型是一种用于评估生产率变化的模型,它可以用于分析不同时间或不同地区的经济绩效差异。

该模型结合了Malmquist 指数和Luenberger 生产率指标,可以分解为技术效率变化和技术进步两个主要部分。

伊辛模型MATLAB源代码

伊辛模型MATLAB源代码

伊辛(Ising)模型MATLAB源程序代码clear all;close all;clcN=16;Te=[0:0.2:5];imax=500000;MT=zeros(1,length(Te));Cv=zeros(1,length(Te));xx=zeros(1,length(Te));L=1;k=1;for T=TeMt=zeros(1,imax);step=100000;spin=round(rand(N)).*2-1;subplot(3,2,1)gspin=pcolor(spin);colormap(gray(2));title('自旋分布情况');Mt=zeros(1,imax);subplot(3,2,2)gMt=plot([1:imax],Mt);axis([0 imax -1.1 1.1]);title('固定温度下的磁矩');%xlabel('时间');%ylabel('Mt');%initial energyidx=repmat({':'},ndims(spin),1);idx{1}=[N 1:N-1];Etotal=spin(idx{:});idx=repmat({':'},ndims(spin),1);idx{2}=[N 1:N-1];Etotal=Etotal+spin(idx{:});idx=repmat({':'},ndims(spin),1);idx{1}=[2:N 1];Etotal=Etotal+spin(idx{:});idx=repmat({':'},ndims(spin),1);idx{2}=[2:N 1];Etotal=Etotal+spin(idx{:});Etotal=spin.*Etotal;Et0=zeros(1,imax);Et0(1)=-sum(sum(Etotal))/2; subplot(3,2,3);gEt=plot([1:imax],Et0);Ett=Et0(1);title('固定温度下的能量');%xlabel('时间');%ylabel('Et');%energy endsubplot(3,2,4);gMT=plot(Te,MT,'o-');axis([0 max(Te) 0 1.1]);title('磁矩随温度T的变化');%xlabel('T');%ylabel('MT');subplot(3,2,5)gCv=plot(Te,Cv,'o-');%axis([0 max(Te) 0 700]);title('热容随温度T的变化');%xlabel('T');%ylabel('Cv');subplot(3,2,6)gxx=plot(Te,xx,'o-');%axis([0 max(Te) 0 0.1]);title('磁导率随温度T的变化');%xlabel=('T');%ylabel=('X');for i=2:imaxa=ceil(rand(1,2).*N);if a(1)==1xl=N;elsexl=a(1)-1;endif a(1)==Nxr=1;elsexr=a(1)+1;endif a(2)==1yd=N;elseyd=a(2)-1;endif a(2)==Nyu=1;elseyu=a(2)+1;endds=spin(a(1),yu)+spin(a(1),yd)+spin(xl,a(2))+spin(xr,a(2)); de=2.*spin(a(1),a(2)).*ds;if rand<exp(-de/T)spin(a(1),a(2))=-spin(a(1),a(2));Ett=Ett+de;endMt(i)=sum(sum(spin))/(N^2);Et0(i)=Ett;if mod(i,step)==0set(gspin,'Cdata',double(spin));set(gMt,'Ydata',Mt);set(gEt,'Ydata',Et0);drawnowendendCv(L)=1/(T^2).*(mean(Et0.^2)-mean(Et0)^2);xx(L)=1/T.*(mean(Mt.^2)-mean(Mt)^2);MT(L)=mean(abs(Mt([0.9*imax:imax])));set(gMT,'Ydata',MT);drawnowset(gCv,'Ydata',Cv);drawnowset(gxx,'Ydata',xx);drawnowL=L+1;end经过本人亲自测试,运行正确。

Matlab产生IGES文件代码【matlab源码】

Matlab产生IGES文件代码【matlab源码】

Matlab产生IGES文件代码【matlab源码】毕业论文,设计,题目学院学院专业学生姓名学号年级级指导教师毕业教务处制表毕业Matlab产生IGES文件代码一、程序说明本团队长期从事matlab编程与仿真工作,擅长各类毕业设计、数据处理、图表绘制、理论分析等,程序代做、数据分析具体信息联系二、程序示例%产生数据并写入iges文件clcclearissearch=2;%0表示只写点,1表示只写线,2表示点线都写fprintf('正在写文件。

\n');%产生正弦波文件x=0:0.1:10;x=x';y=sin(x);Data=[x y];Data(1,3)=0;write_iges('iges_sin.igs',Data,issearch)%产生peaks文件[x,y,z]=peaks(30); xx=x(:);yy=y(:);zz=z(:);Data=[xx zz yy];write_iges('iges_peaks.igs',Data,issearch)%产生抛物线文件[x,y]=meshgrid(-1:0.1:1);z=x.^2+y.^2;xx=x(:);yy=y(:);zz=z(:);Data=[xx yy zz];write_iges('iges_paowu.igs',Data,issearch)%产生解释文件xx=[10;0];yy=[0;10];zz=[0;0];Data=[xx yy zz];write_iges('igesforexplain.igs',Data,0) %产生Matlab图标[x,y]=meshgrid(linspace(-1,1,16)); z=membrane;z=z(1:2:end,1:2:end,1:2:end); xx=x(:);yy=y(:);zz=z(:);Data=[xx zz yy];write_iges('iges_matlab.igs',Data,issearch)fprintf('写文件结束。

一维伊辛模型MATLAB

一维伊辛模型MATLAB

一维伊辛模型MATLAB源代码clear;clc;%本次模拟中,J=1,k=1kT=0.1:0.1:5;kT_len=length(kT);num=102;Ising=sign(2*rand(1,num)-1);%初始构型N=2*10^5;%------初始化,预分配内存------%E_total=zeros(1,kT_len);E_squ=zeros(1,kT_len);C=zeros(1,kT_len);M_total=zeros(1,kT_len);for i=1:kT_lenE_sum=0;E_squ_sum=0;M_sum=0;for turn=1:NE_old=E_calc(Ising);ri=randi([1,num],1);if ri==1ri=num-1;endif ri==numri=2;endIsing(ri)=-Ising(ri);%选择其中一个原子使其自旋态相反E_new=E_calc(Ising);E_diff=E_new-E_old;if E_diff>0 && exp(-E_diff/kT(i))<=rand()Ising(ri)=-Ising(ri);%保持原来的构型不变endif turn>N/2E=E_calc(Ising);M=sum(Ising)./num;E_sum=E_sum+E;E_squ_sum=E_squ_sum+E^2;M_sum=M_sum+M;endendE_total(i)=E_sum/(N/2);E_squ(i)=E_squ_sum/(N/2);C(i)=(E_squ(i)-E_total(i).^2)/(kT(i).^2);M_total(i)=M_sum/(N/2);end%-------画图部分--------%subplot(131),plot(kT,E_total,'--*b'),xlabel('kT'),ylabel('能量');title('温度与能量关系图');subplot(132),plot(kT,C,'--*r'),xlabel('kT'),ylabel('热容');title('温度与热容关系图');subplot(133),plot(kT,M_total,'--*k'),xlabel('kT'),ylabel('磁化强度'); title('温度与磁化强度关系图');%------函数文件部分,需另外新建函数文件,并保存------% function [E] = E_calc(Ising)%用于计算该构型下的能量len=length(Ising);E=0;for i=2:len-1E=E-Ising(i)*(Ising(i-1)+Ising(i+1))/2;end。

数学模型程序代码-Matlab-姜启源-第三章-简单的优化模型

数学模型程序代码-Matlab-姜启源-第三章-简单的优化模型

数学模型程序代码-Matlab-姜启源-第三章-简单的优化模型第3章简单的优化模型1. 生猪的出售时机p63~65目标函数(生猪出售纯利润,元):Q(t) = ( 8 – g t )( 80 + rt ) – 4t–640其中,t≥0为第几天出售,g为每天价格降低值(常数,元/公斤),r为每天生猪体重增加值(常数,公斤)。

求t使Q(t)最大。

1.1(求解)模型求解p63(1) 图解法绘制目标函数Q(t) = ( 8 – g t )( 80 + rt ) – 4t–640的图形(0 ≤t≤ 20)。

其中,g=0.1, r=2。

从图形上可看出曲线Q(t)的最大值。

(2) 代数法对目标函数Q(t) = ( 8 – g t )( 80 + rt ) – 4t–640用MATLAB求t使Q(t)最大。

其中,r, g是待定参数。

(先对Q(t)进行符号函数求导,对导函数进行符号代数方程求解)然后将代入g=0.1, r=2,计算最大值时的t和Q(t)。

要求:①编写程序绘制题(1)图形。

②编程求解题(2).③对照教材p63相关内容。

相关的MATLAB函数见提示。

★要求①的程序和运行结果:t=0:1:30;g=0.1;r=2;Q=(8-g.*t).*(80+r.*t)-4.*t-640;plot(t,Q)★要求②的程序和运行结果:程序:syms g t r ;Q=(8-g.*t).*(80+r.*t)-4.*t-640;q=diff(Q,t);q=solve(q);g=0.1;r=2;tm=eval(q)Q=(8-g.*tm).*(80+r.*tm)-4.*tm-640 运行结果:1.2(编程)模型解的的敏感性分析p63~64对1.1中(2)所求得的符号表达式t(r,g),分别对g和r进行敏感性分析。

(1) 取g=0.1,对t(r)在r=1.5:0.1:3上求r与t的关系数据,绘制r与t的关系图形(见教材p65)。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

伊辛(Ising)模型MATLAB源程序代码clear all;
close all;
clc
N=16;
Te=[0:0.2:5];
imax=500000;
MT=zeros(1,length(Te));
Cv=zeros(1,length(Te));
xx=zeros(1,length(Te));
L=1;k=1;
for T=Te
Mt=zeros(1,imax);
step=100000;
spin=round(rand(N)).*2-1;
subplot(3,2,1)
gspin=pcolor(spin);
colormap(gray(2));
title('自旋分布情况');
Mt=zeros(1,imax);
subplot(3,2,2)
gMt=plot([1:imax],Mt);
axis([0 imax -1.1 1.1]);
title('固定温度下的磁矩');
%xlabel('时间');
%ylabel('Mt');
%initial energy
idx=repmat({':'},ndims(spin),1);
idx{1}=[N 1:N-1];
Etotal=spin(idx{:});
idx=repmat({':'},ndims(spin),1);
idx{2}=[N 1:N-1];
Etotal=Etotal+spin(idx{:});
idx=repmat({':'},ndims(spin),1);
idx{1}=[2:N 1];
Etotal=Etotal+spin(idx{:});
idx=repmat({':'},ndims(spin),1);
idx{2}=[2:N 1];
Etotal=Etotal+spin(idx{:});
Etotal=spin.*Etotal;
Et0=zeros(1,imax);
Et0(1)=-sum(sum(Etotal))/2; subplot(3,2,3);
gEt=plot([1:imax],Et0);
Ett=Et0(1);
title('固定温度下的能量');
%xlabel('时间');
%ylabel('Et');
%energy end
subplot(3,2,4);
gMT=plot(Te,MT,'o-');
axis([0 max(Te) 0 1.1]);
title('磁矩随温度T的变化');
%xlabel('T');
%ylabel('MT');
subplot(3,2,5)
gCv=plot(Te,Cv,'o-');
%axis([0 max(Te) 0 700]);
title('热容随温度T的变化');
%xlabel('T');
%ylabel('Cv');
subplot(3,2,6)
gxx=plot(Te,xx,'o-');
%axis([0 max(Te) 0 0.1]);
title('磁导率随温度T的变化');
%xlabel=('T');
%ylabel=('X');
for i=2:imax
a=ceil(rand(1,2).*N);
if a(1)==1
xl=N;
else
xl=a(1)-1;
end
if a(1)==N
xr=1;
else
xr=a(1)+1;
end
if a(2)==1
yd=N;
else
yd=a(2)-1;
end
if a(2)==N
yu=1;
else
yu=a(2)+1;
end
ds=spin(a(1),yu)+spin(a(1),yd)+spin(xl,a(2))+spin(xr,a(2)); de=2.*spin(a(1),a(2)).*ds;
if rand<exp(-de/T)
spin(a(1),a(2))=-spin(a(1),a(2));
Ett=Ett+de;
end
Mt(i)=sum(sum(spin))/(N^2);
Et0(i)=Ett;
if mod(i,step)==0
set(gspin,'Cdata',double(spin));
set(gMt,'Ydata',Mt);
set(gEt,'Ydata',Et0);
drawnow
end
end
Cv(L)=1/(T^2).*(mean(Et0.^2)-mean(Et0)^2);
xx(L)=1/T.*(mean(Mt.^2)-mean(Mt)^2);
MT(L)=mean(abs(Mt([0.9*imax:imax])));
set(gMT,'Ydata',MT);
drawnow
set(gCv,'Ydata',Cv);
drawnow
set(gxx,'Ydata',xx);
drawnow
L=L+1;
end
经过本人亲自测试,运行正确。

相关文档
最新文档