用Matlab求水仙花数,完美数,回文数和亲合数
数值分析 - 第1章 MATLAB入门
（2）>> u=[x,y] u = 7 4 3 -1 -2 -3
3、解： sum=0;a=[4 -1 2 -8 4 5 -3 -1 6 for i=1 : length(a) if a(i)>0, sum=sum+a(i); end end sum sum = 21
-7];
7 5 3 （4）>> d=b(4:-1:1) d = 3 5 7 9 （5）>> e=sort(b) e =
1水仙花数是一个3位自然数应当满足在100999之间2由于还需要满足各位数的立方和等于该数本身因此需要求出该数的每位数字另外我们还要用到取模函数mod和取整函数fix
1.3
习题课
1、通过使用 help 命令来学习使用 date、clock 命令，了解它们的显示格式。 解答：本题的主要目的让学生学习使用 help 命令来学习一些命令的使用方法。 在 MATLAB 中输入 help date 或者 help clock，会得到 date 和 clock 命令的一些说明。结果如下： DATE Current date as date string. S = DATE returns a string containing the date in dd-mmm-yyyy format. 它表示 date 命令的显示结果是按照：天、月、年的格式显示的。 CLOCK Current date and time as date vector. CLOCK returns a six element date vector vector containing the current time and date in decimal form: CLOCK = [year month day hour minute seconds] The first five elements are integers. The seconds element is accurate to several digits beyond the decimal point. FIX(CLOCK) rounds to integer display format. 这是 help clock 的结果，它提示了 clock 命令的显示格式为：年、月、日、小时、分钟、秒。 然后我们分两次输入 date 和 clock 命令后就会得到各自的结果。通过 help 命令的提示我们就可以轻松看懂结果 的意思了。 >> date ans = 22-Oct-2007 >> clock ans = 1.0e+003 * 2.0070 0.0100 0.0220 0.0170 0.0130 0.0313 2、编写一个计算公式，在 MATLAB 中计算一个半径为 3 的圆的周长。 解答：先写出周长计算公式：c=2πr，然后将里边的数学运算法转换成 MATLAB 的运算符为：c=2*pi*r。 下面我们介绍两组语句，它们都可以输出周长的值。 （1） >> r=3; >> c=2*pi*r; >> c c= 18.8496 （2） >> r=3; >> 2*pi*r ans =
matlab计算公式
matlab计算公式MATLAB是一种强大的数学计算软件，广泛应用于科学与工程领域。
它提供了丰富的计算功能，包括数值计算、矩阵运算、绘图以及数据分析等。
本文将介绍MATLAB中常用的计算公式及其使用方法。
一、数值计算1. 四则运算在MATLAB中，进行加减乘除等基本运算非常简单。
例如，要计算两个数的和，可以使用如下公式：```matlabsum = a + b;```其中，`a`和`b`是待计算的数值，`sum`则为计算结果。
2. 幂运算如果需要进行幂运算，可以使用`^`操作符，例如：```matlabresult = x ^ n;```这个公式表示将`x`的`n`次方赋值给`result`。
3. 开方运算开方运算可以通过调用`sqrt`函数来实现：```matlabroot = sqrt(x);````x`为待开方的数值，`root`为计算结果。
二、矩阵运算在MATLAB中，矩阵运算也是常见且重要的功能之一。
1. 矩阵相加要计算两个矩阵的和，可以使用如下公式：```matlabresult = matrix1 + matrix2;```其中，`matrix1`和`matrix2`是待相加的矩阵，`result`为计算结果。
2. 矩阵相乘矩阵相乘可以使用`*`操作符，例如：```matlabresult = matrix1 * matrix2;```这个公式表示将`matrix1`和`matrix2`进行乘法运算，并将结果赋值给`result`。
3. 转置矩阵如果需要计算矩阵的转置，可以使用`'`符号，例如：```matlabtranspose = matrix';```这个公式表示将`matrix`进行转置操作，并将结果赋值给`transpose`。
三、绘图功能MATLAB提供了丰富的绘图功能，可以用于生成各种类型的图形，如折线图、散点图、柱状图等。
1. 折线图要绘制折线图，可以使用`plot`函数，例如：```matlabx = 0:0.1:2*pi;y = sin(x);plot(x, y);```这个公式表示生成从0到2π的一系列数据，并计算其对应的正弦值，然后通过`plot`函数将结果绘制成折线图。
第二章 MATLAB数值计算
{矩阵（Matrix）：是一个矩形的数组，即二维 数组，其中向量和标量都是矩阵的特例，0×0矩 阵为空矩阵([])。
{数组（Array）：是指n维的数组，为矩阵的延 伸，其中矩阵和向量都是数组的特例。
复数（Complex）
以m×n的矩阵a为例，若元素a(i,j)则 对应的“单下标”为s= (j-1)×m＋i。
a(1,2) a(4)
a(2,3) a(8)
2. 子矩阵块的产生
子矩阵是从对应矩阵中取出一部 分元素构成的。 （1）全下标方式
例 设矩阵A是已知的66矩阵
A(1:3,2:5) A的前3行，及第2到5列元素构成的子矩阵
数值采用十进制表示，可以用带小 数点的形式直接表示，也可以用科学计 数法，数值的表示范围是10-309～10309 。 例如：
-2、5.67、2.56e-56(表示2.56×10-56)、 4.68e204(表示4.68×10204)
矩阵和数组的概念
{标量（Scalar）：是指1×1的矩阵，即只含一个 数的矩阵。
表示数组A和B中的对应元素相乘。A和B数组必 须大小相同，除非其中有一个是标量。
(3) 矩阵和数组的除法(division) 矩阵的除法运算表达式有两种：“A\B”和“A/B”， 运算符为“\”和“/”分别表示左除和右除。 X=B/A是X*A=B的解，B/A=B*A-1。 X=A\B是方程A*X=B的解， A\B=A-1*B。 数组的除法运算表达式有两种：“A.\B”和“A./B”， 表示数组相应元素左除和右除。 A和B数组必须大小相同，除非其中有一个是标量。
产生均匀分布的随机矩阵，元素取值 范围0.0～1.0。 产生正态分布的随机矩阵
matlab编程题整理
作业一1输出x,y两个中值较大的一个值x=input(‘x’);y=input(‘y’);if x>yxelseyend2输入x，计算y的值。
计算函数的值y=x+1,x<0,y=2x-1,x≧0x=input错误！未指定书签。
(‘x);if x<0y=x+1elsey=2*x-1end3输入一学生成绩，评定其等级，方法是：90~100分为“优秀”，80~89分为“良好”，70~79分为“中等”，60~69分为“及格”，60分为“不合格”x=input(‘x’)if x>100|x<0y=’输入错误’elseif x>=90y=’优秀’elseif x>=80y=’良好’elseif x>=70y=’中等’elseif x>=60y=’及格’elsey=’不合格’emd4某超市节日期间举办购物打折的促销活动，优惠办法是：每位顾客当天一次性购物在100元以上者，按九五折优惠；在200元以上者，按九折优惠；在300元以上者，按八五折优惠；在500元以上者，按八折优惠。
x=input(‘x’);if x>=500y=x*0.8elseif x>=300y=x*0.85elseif x>=200y=x*0.9elseif x>=100y=x*0.95else y=xend 5编程计算：s=1+2+3+…+100sum=0;for i=1:100sum=sum+i;endsum引申1!-2!+3!-4!+5!- (99)sum=0;for i=1:99pdr=1;for k=1:ipdr=pdr*k;endsum=sum+pdr*(-1)^(i-1);endsum引申1*2*3*4*……*100sum=1；for i=1:100sum=sum*iendsum6计算1~100的奇数和sum=0;for i=1:2:100sum=sum+i;endsum7百元买百鸡问题。
(完整版)第一讲Matlab基本数值计算
第一讲Matlab 基本数值计算一、矩阵在Matlab中，一个矩阵可以使数学意义上的矩阵，也可以是标量或者向量。
对于一个标量（一个数）可以将之作为11⨯的矩阵，而向量（一行或一列）则可以认为是1n⨯⨯或者1n⨯的矩阵。
另外，一个00矩阵在Matlab中被认为是空矩阵，用“[]”表示。
1、矩阵的创建矩阵的创建可以有以下几种形式⑴直接输入>> A=[1 2 3;4 3 7;2 4 1]注意：每行间的元素用逗号或空格分开，行与行之间用分号或回车分开，矩阵标示是一对中括号[ ]。
也可以采用数组编辑器（Array Editor）像在Excel电子表格中据那样输入数据。
⑵通过语句和函数产生常用的特殊矩阵：zeros:全零矩阵,ones：全1矩阵，eye：单位矩阵，rand:随机矩阵，diag:对角阵等。
例：>> A=ones(3,4)>> E=eye(3)>> D=diag([3 5 2])⑶对矩阵进行裁剪或拼接⑷从外部文件装入数据外部数据文件可以是以保存的Matlab工作空间，也可以是文本（.txt）文件，或者是电子表格创建的文件（.xls）.例：已知一个文本格式的数据文件E:\Mathmodel\data1.txt>> load e:\Mathmodel\data1.txt得到一个变量名与文件名相同的矩阵（data1）。
注意：文件的扩展名不能省略。
例：已知一个Excel文件的路径为E:\Mathmodel\data2.xlsa. 缺省操作:>> NUMBER=xlsread('E:\Mathmodel\data2.xls')>>[NUMBER,TXT]=xlsread('E:\Mathmodel\data2.xls')默认操作是从第一个工作表（sheet1）中提取数据。
b. 从指定的工作表（而不是第一个）中提取数据：>> NUMBER=xlsread('E:\Mathmodel\data2.xls','S2')或者>> NUMBER=xlsread('E:\Mathmodel\data2.xls',2)c.从指定的工作表中读取指定区域的数据：>> NUMBER=xlsread('E:\Mathmodel\data2.xls',2,'g3:i8')2、Matlab的矩阵运算⑴基本运算矩阵的加（+）、减（-）、乘（*）、乘方（^）运算法则与代数中的定义完全一致。
MATLAB软件应用 解题代码 线性代数
2%求全部4位水仙花数m=1000:9999%生成所有四位数的m向量m1=rem(m,10)%求m除以10的余数赋给m1m2=rem(fix(m/10),10)%先调用fix，对m/10的结果取整，再求取整后的数/10的余数m3=rem(fix(m/100),10)%先调用fix，对m/100的结果取整，再求取整后的数/10的余数m4=fix(m/1000)%对m/1000取整赋予m4k=find(m==m1.*m1.*m1.*m1+m2.*m2.*m2.*m2+m3.*m3.*m3.*m3+m4.*m4.*m4.*m4) %调用find函数，在向量m中找到水仙花数的序号赋给变量ks=m(k)%输出水仙花数%反序x='Kat is beautiful'%生成字符串result=strrep(x,'is',x(6:-1:5))%调用strrep函数将x中的子字符串is替换为利用冒号表达式表示的反序，初项为6，末项为5，步长为-13A=magic(7)%产生一个7阶魔方矩阵A([2 6],:)=A([6 2],:)%交换矩阵A的第一行和第二行B=A%得到矩阵Bsum(B(1,:))%求B矩阵第一行元素之和sum(B(1,:))==sum(B(2,:))%验证B矩阵各行元素之和是否相等sum(B(2,:))==sum(B(3,:))sum(B(3,:))==sum(B(4,:))sum(B(4,:))==sum(B(5,:))sum(B(5,:))==sum(B(6,:))sum(B(6,:))==sum(B(7,:))sum(B(:,1))%求B矩阵第一列元素之和sum(B(:,1))==sum(B(:,2))%验证B矩阵各列元素之和是否相等sum(B(:,2))==sum(B(:,3))sum(B(:,3))==sum(B(:,4))sum(B(:,4))==sum(B(:,5))sum(B(:,5))==sum(B(:,6))sum(B(:,6))==sum(B(:,7))sum(diag(B))%求B矩阵主对角线元素之和B1=flipud(B)%将B矩阵上下翻转，B矩阵的副对角线就变成了B1矩阵的主对角线sum(diag(B1))%求B1矩阵主对角线元素之和即求B矩阵副对角线元素之和sum(B(1,:))==sum(B(:,1))%验证B矩阵每行、每列元素之和是否相等sum(diag(B))==sum(diag(B1))%验证B矩阵主、副对角线元素之和是否相等sum(B(1,:))==sum(diag(B))%验证B矩阵每行、每列元素之以及主、副对角线元素之和是否相等B(:,[2 6])=A(:,[6 2])%交换矩阵B的第一列和第二列C=B%得到矩阵Csum(C(1,:))%求C矩阵第一行元素之和sum(C(1,:))==sum(C(2,:))%验证C矩阵各行元素之和是否相等sum(C(2,:))==sum(C(3,:))sum(C(3,:))==sum(C(4,:))sum(C(4,:))==sum(C(5,:))sum(C(5,:))==sum(C(6,:))sum(C(6,:))==sum(C(7,:))sum(C(:,1))%求C矩阵第一列元素之和sum(C(:,1))==sum(C(:,2))%验证C矩阵各列元素之和是否相等sum(C(:,2))==sum(C(:,3))sum(C(:,3))==sum(C(:,4))sum(C(:,4))==sum(C(:,5))sum(C(:,5))==sum(C(:,6))sum(C(:,6))==sum(C(:,7))sum(diag(C))%求C矩阵主对角线元素之和C1=flipud(C)%将C矩阵上下翻转，C矩阵的副对角线就变成了C1矩阵的主对角线sum(diag(C1))%求C1矩阵主对角线元素之和即求C矩阵副对角线元素之和sum(C(1,:))==sum(C(:,1))%验证C矩阵每行、每列元素之和是否相等sum(diag(C))==sum(diag(C1))%验证C矩阵主、副对角线元素之和是否相等sum(C(1,:))==sum(diag(C))%验证C矩阵每行、每列元素之以及主、副对角线元素之和是否相等4x=[2,4,8;0,-6,-4;8,1,7]%自定义一个非奇异的3阶方阵xn=norm(x) %直接求x的2-范数max(eig(x.'*x))^(1/2)%求矩阵x的转置乘矩阵x本身所得的矩阵的最大特征值的平方根%两者的结果相等x=[0,6,6,1,1,6,6,1,1,6,6,0;0,0,1,1,5,5,6,6,10,10,11,11]%生成结点坐标矩阵x，第一行为横坐标，第二行为纵坐标A1=[2,1;1,1]%定义变换矩阵A1A2=[0,-1;1,0]%定义变换矩阵A2A3=[2,0;0,2]%定义变换矩阵A3y1=A1*x%再利用A1对x进行变换，得到y1矩阵y2=A2*x%再利用A2对x进行变换，得到y2矩阵y3=A3*x%再利用A3对x进行变换，得到y3矩阵%分别绘制变换前后的图形subplot(2,2,1)fill(x(1,:),x(2,:),'r')subplot(2,2,2)fill(y1(1,:),y1(2,:),'r')subplot(2,2,3)fill(y2(1,:),y2(2,:),'r')subplot(2,2,4)fill(y3(1,:),y3(2,:),'r')5%判断n是否为质数n1 = input('请输入一个不大于100的整数：','s');n = str2num(n1);%因为默认输入的是ASCII码值，所以要有这一步，转换成数字if rem(n,2)==0 & n/2~=1 %表示如果输入的数字被2整除，则输出'您所输入的数为合数'，且把2排除disp('您所输入的数为合数')elseif rem(n,3)==0 & n/3~=1%条件表示如果输入的数字是奇数，且被3整除，则输出'您所输入的数为合数'，且把3排除disp('您所输入的数为合数')elseif rem(n,3)~=0 & rem(n,5)==0 & n/5~=1%条件表示如果输入的数字不被2、3整除，但被5这整除，则输出'您所输入的数为合数'，且把5排除disp('您所输入的数为合数')elseif rem(n,3)~=0 & rem(n,5)~=0 & rem(n,7)==0 & n/7~=1%条件表示如果输入的数字不被2、3、5整除，但被7这整除，则输出'您所输入的数为合数'，且把7排除disp('您所输入的数为合数')else %以上条件都不满足，则输出num本身disp(n)end%判断n是否为两个质数的乘积n1 = input('请输入一个不大于100的整数：','s');n = str2num(n1);%因为默认输入的是ASCII码值，所以要有这一步，转换成数字if rem(n,2)==0 & isprime(n/2) %表示如果输入的数字被2整除且除以2后为质数，则判定输入的数为两个质数的乘积disp([num2str(n),'是2和',num2str(n/2),'的乘积']) %输出‘n是2和n/2的乘积’elseif rem(n,3)==0 & isprime(n/3) %表示如果输入的数字被3整除且除以3后为质数，则判定输入的数为两个质数的乘积disp([num2str(n),'是3和',num2str(n/3),'的乘积']) %输出‘n是3和n/3的乘积’elseif rem(n,3)~=0 & rem(n,5)==0 & isprime(n/5) %表示如果输入的数字被5整除且除以5后为质数，则判定输入的数为两个质数的乘积disp([num2str(n),'是5和',num2str(n/5),'的乘积']) %输出‘n是5和n/5的乘积’elseif rem(n,3)~=0 & rem(n,5)~=0 & rem(n,7)==0 & isprime(n/7) %表示如果输入的数字被7整除且除以7后为质数，则判定输入的数为两个质数的乘积disp([num2str(n),'是7和',num2str(n/7),'的乘积']) %输出‘n是7和n/7的乘积’else %以上条件都不满足，则输出'您输入的数不能分解为两个质数的乘积' disp('您输入的数不能分解为两个质数的乘积')end6%用辛普森法计算四分之一圆的面积，从而计算圆周率的近似值disp('以下是辛普森法求圆周率')a=0; %给a赋值为0b=1; %给b赋值为1n=input('n=?'); %输入n的值，即把区间分成n等份h=(b-a)/n; %计算每个区间的宽度，赋给变量hx=a:h:b %产生自变量向量x，x中包含了要取的n+1个自变量的值f=sqrt(1-x.*x);%得到一个与x对应的函数值向量赋给变量f,(这里一定要用点乘)s=[] %定义一个空矩阵sfor k =1:2:n %循环n次，步长为2s1=(f(k)+4*f(k+1)+f(k+2))*h/3; %求以抛物线弧段为曲边，以[Xk,Xk+2]为底的曲边梯形面积s=[s,s1]; %将s1添加到矩阵s中，循环n次后，s里面的元素分别时endpai=4*sum(s) %调用sum函数求s各元素之和，即整个曲边梯形的面积，再乘以4即为pi8figure(1)u1=linspace(-pi,pi,100);u2=linspace(-2*pi,2*pi,200);u3=linspace(0,2*pi,400);plot(u1,sin(2*u1)+1,'b-',u2,cos(u2),'r--',u3,sin(u3)+cos(2*u3),'k-.')figure(1)x=linspace(0,0.3);y=x.*sin(1./x);plot(x,y,'b')title('figure(1)')figure(2)r=1;t=0:6*pi;x1=r*(t-sin(t));y1=r*(1-cos(t));plot(x1,y1,'r')title('figure(2)')figure(3)syms x yezplot('4*x^2+9*y^2=36')grid ontitle('figure(3)')10figuresubplot(1,2,1);t=-1:0.1:1;[X,Y]=meshgrid(t);Z=X.^3+Y.^2;mesh(X,Y,Z);hold onsurf(X,zeros(size(Y)),Y.^3)title('{z=x^{3}+y^{2}}与x0z平面的相交图')subplot(1,2,2);t=-1:0.1:1;[X,Y]=meshgrid(t);Z=X.^2+Y.^5;mesh(X,Y,Z);hold onsurf(zeros(size(Y)),Y,X.^5)title('{z=x^{2}+y^{5}}与y0z平面的相交图')11%画月亮t=linspace(0,2*pi,100);subplot(2,2,1)x=sin(t);y=cos(t);p1=y>0.5;y(p1)=NaN;plot(x,y)hold onx11=sqrt(3)*sin(t)./2;y11=sqrt(3)*(0.5+cos(t))./2;p11=y11>0.5;y11(p11)=NaN;plot(x11,y11)axis([-1.1,1.1,-1.1,1.1])axis squaregrid onsubplot(2,2,2)x=sin(t);y=cos(t);p2=x>0.5;x(p2)=NaN;plot(x,y)hold onx11=sqrt(3)*(0.5+sin(t))./2;y11=sqrt(3)*cos(t)./2;p11=x11>0.5;x11(p11)=NaN;plot(x11,y11)axis([-1.1,1.1,-1.1,1.1])axis squaregrid onsubplot(2,2,3)x=sin(t);y=cos(t);p3=y<-0.5;y(p3)=NaN;plot(x,y)hold onx11=-sqrt(3)*sin(t)./2;y11=-sqrt(3)*(0.5+cos(t))./2; p11=y11<-0.5;y11(p11)=NaN;plot(x11,y11)axis([-1.1,1.1,-1.1,1.1])axis squaregrid onsubplot(2,2,4)x=sin(t);y=cos(t);p4=x<-0.5;x(p4)=NaN;plot(x,y)hold onx11=-sqrt(3)*(0.5+sin(t))./2; y11=-sqrt(3)*cos(t)./2;p11=x11<-0.5;x11(p11)=NaN;plot(x11,y11)axis([-1.1,1.1,-1.1,1.1])axis squaregrid on13% (a)x1=[1,3,5];p=poly(x1);p=p*3;px=poly2str(p,'x')x2=[2,4];q=poly(x2);q=q*-2;qx=poly2str(q,'x')% (b)A=[2,3;1,1];B=[2,1;1,1];pA=polyvalm(p,A)qA=polyvalm(q,A)pB=polyvalm(p,B)qB=polyvalm(q,B)%点乘的话，p(A)q(A)=q(B)p(A) pA.*qBqB.*pA%乘的话，p(A)q(A)和q(B)p(A)不相等pA*qBqB*pA% (c)a=0;b=10;f=diff(polyval(polyint(p),[a b]))p=[3,-27,69,-45]px=poly2str(p,'x')q=[-2,12,-16]qx=poly2str(q,'x')r=[1,1]rx=poly2str(r,'x')q1=[0,0,q]f=conv(p,r)+q1x=roots(f)14%算积分%【a】fun1=@(x) x+x.^2+x.^3;I1=rectangular(fun1,-1,1,1000)x=linspace(-1, 1, 100);y=x+x.^2+x.^3;I2=trapz(x,y)I3=quad(fun1,-1,1)%【b】fun2=@(x,y) sin(y.*(x+y)./(x.^2+4)) I=quad2d(fun2,1,10,1,10)function f=rectangular(fun,a,b,n)%矩形法:h=(b-a)/n;x=a:h:b;y=x;for i=2:n+1y(i)=fun((x(i)+x(i-1))/2);endf=h*sum(y(1:end));endsyms x y zeq1 = 3*x + 2*y - z == 4;eq2 = x - y + z == 1;eq3 = -x + 4*y + 5*z == 8;[x,y,z] = solve([eq1,eq2,eq3], [x,y,z])14%算积分%【a】fun1=@(x) x+x.^2+x.^3;I1=rectangular(fun1,-1,1,1000)x=linspace(-1, 1, 100);y=x+x.^2+x.^3;I2=trapz(x,y)I3=quad(fun1,-1,1)%【b】fun2=@(x,y) sin(y.*(x+y)./(x.^2+4)) I=quad2d(fun2,1,10,1,10)function f=rectangular(fun,a,b,n)%矩形法:h=(b-a)/n;x=a:h:b;y=x;for i=2:n+1y(i)=fun((x(i)+x(i-1))/2);endf=h*sum(y(1:end));endsyms x y zeq1 = 3*x + 2*y - z == 4;eq2 = x - y + z == 1;eq3 = -x + 4*y + 5*z == 8;[x,y,z] = solve([eq1,eq2,eq3], [x,y,z])。
水仙花数通用算法
break
if flag:
print(tempnumber)
matchlist = [0 for x in range(10)]
flag = 1
以上就是代码了。。。
如有错误或需要讨论。。。
请移步：
/valbendan/article/details/7461243
max = 10 ** weishuo
if weishuo is 1:
min = 0
else :
min = 10 ** (weishuo-1)
max9 = max // number[9][1]
min9 = min // number[9][1]
for s[9] in range(min9,max9e(0,weishuo-s[9]+1):
for s[7] in range(0,weishuo-s[9]-s[8]+1):
for s[6] in range(0,weishuo-s[9]-s[8]-s[7]+1):
for s[5] in range(0,weishuo-s[9]-s[8]-s[7]-s[6]+1):
number = [[x * i ** weishuo for x in range(weishuo+1)]for i in range(10)]
matchlist = [0 for x in range(10)]
s = [x for x in range(10)]
bignumber = 0
flag = 1
for s[1] in range(0,weishuo-s[9]-s[8]-s[7]-s[6]-s[5]-s[4]-s[3]-s[2]+1):
matlab实验题
实验一、Matlab 基本操作1、 计算9.248.26107sin 369.12÷⎪⎭⎫ ⎝⎛π+的值 2、 产生一个5阶魔术方阵，并执行如下操作：（1） 将矩阵的第2行3列元素赋值给变量c（2） 将由矩阵第2，3，4行第3，5列构成的子矩阵赋值给变量d3、给出区间[0，1]上的6个等分点数据。
4、建立如下矩阵 （1）1010200700020070002007⨯⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭ （2）1010010101001010100⨯⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭ 5、设有分块矩阵⎥⎦⎤⎢⎣⎡=⨯⨯⨯⨯22322333S O R E A ，其中E,R,O,S 分别为单位阵、随机阵、零阵和对角阵，试通过数值计算比较A*A 与A.*A 的区别。
实验二、Matlab 程序设计1、 已知函数⎪⎩⎪⎨⎧≤≤<≤<≤-+=2110,101,1)(2x x x x x x f 计算)5.1(),5.0(),1(f f f -. 2、 用for-end 循环语句求：100！和∑=1001i i 。
3、 用while-end 循环语句求不超过1000的偶数之和与奇数之和。
4、 建立一个命令M-文件：求所有的“水仙花数”，所谓“水仙花数”是指一个三位数，其各位数字的立方和等于该数本身。
例如，153是一个水仙花数，因为153=13+53+33。
5、建立如下矩阵（1）10101000120011100⨯⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭ （2）20072000002007300000200740000020075000002007⎛⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭实验三、使用Matlab 作图1、 在同一坐标系下面画出)5.0sin(2.01.0x e y x +=和)5.0cos(2.01.0x e y x +=在区间]2,0[π上的曲线图。
2、 绘制三维螺旋线：].,0[,5.0,sin 2,cos 2π∈===t t z t y t x3、 画出曲面)sin(xy z =的网线图。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
题目:对深圳市人口及医疗的预测编号: 10055队长李荣姓名: 李荣(09041117)王遂(09041124)李春晓(09055025)1.水仙花数:(1)题目若一个三位数自然数的各位数的各位数字的立方和等于该数本身,则称该数为水仙花数,例如,153=1^3+5^3+3^3,所以153为水仙花数,编程计算所有水仙花数.(2)分析问题:本题特点在于一个数的各个位数立方和与其本身的关系很明确且很简单,可以用if语句或find语句来完成.如果用if语句,则又有两种选择:1.命一个三位数,再表示其各个位数的数字;2.命各个位数的数字,再表示该三位数.如果用find语句,其思路与if语句类似.(3)问题求解:方法一:for x=100:1:999a=fix(x/100);b=fix(x/10-10*a);c=x-100*a-10*b;if x==a^3+b^3+c^3xendend方法二:for x=1:1:9for y=0:1:9for z=0:1:9if x^3+y^3+z^3==100*x+10*y+zm=100*x+10*y+zendendendend方法三:shui=100:999;i=floor(shui/100);j=floor(mod(shui,100)/10);k=floor(mod(shui,10));p=i.^3+j.^3+k.^3;shui(find(p==shui))(4)结论及分析:通过实验,结果正确,证明分析无误.(5)结果:153 370 371 4072.回文数(1)题目:对于一个自然数,若将各位数字倒序排出,加到原数字上,反复多次,若能得到一个从左到右读和从有到左读完全一样的数,则称该自然数能产生回文数.通过编程计算,你能找出多少个能产生回文数的数,又能找到多少不能产生回文数的数,二者的最小数是多少?(2)分析问题:本题关键在于如何将一个数的各个数位的数字倒序排出,并加到原数上.又由题目知必然要用到循环语句.(3)假设:由于将一个数的各个数位的数字倒序排出,并不断加到原数上比较困难,我们假设该数为四位数并假设倒序排出加到原数得到的新数不超过10000,于是有了下面的简单程序for a=0:9;for b=0:9;for c=0:9;for d=0:9;x=1000*a+100*b+10*c+d;while x<10000;if a==d&b==c;x=x;elsex=x+1000*d+100*c+10*b+a;endendendendxendend实际操作时，我们发现：该程序不仅具有假设中的缺陷,而且在实际操作中运行速度很慢,无法得到结果,所以必须将方法优化。
后来我们使用取余和取整操作，借助循环语句while，得出如下程序，突破了本题的关键。
for a=1:100b=0;c=0;x=a;while x>=10b=mod(x,10);x=fix(x/10);c=10*c+b;endc=10*c+x;此程序实现了一个数的倒序操作,并不用考虑位数问题.(4)问题求解：for a=1:100c=0;x=a;while x>=10b=mod(x,10);x=fix(x/10);c=10*c+b;endc=10*c+x;A=a+c;B=0;y=A;while y>=10C=mod(y,10);y=fix(y/10);B=10*B+C;endB=10*B+y;while B~=AA=A+B;m=0;z=A;while z>=10n=mod(z,10);z=fix(z/10);m=10*m+n;endm=10*m+z;B=m;end[a,A]end(5)结论及分析：修改后的程序通过四个while循环语句完成了（1）数字的倒序排出（2）倒序排出数与原数无上界的不断相加（3）相加之和倒序、正序排列所得数的比较，得出了正确的的结果。
运行速度较快.(6)结果：最小的能产生回文数的数是1,100个数中找到了98个能产生回文数的数.其中89，98运算结果为：1.0e+012 * 1.0e+012 *0.0000 8.8132 0.0000 8.8132它们是不能产生回文数的数，则得到两个这样的数，最小的为89..(7)结论推广：在命a的范围时,我们的范围是1到100,而事实上a可以任意取值,故可以通过国修改a的范围算出任意区间上的结果.3.完美数(1)题目:对于一个数而言，若它的真因子之和等于数本身，则称该数为完美数。
通过编程计算，你能找到多少个完美数。
(2)分析问题:问题的关键在于找出一个数的真因子，通过使用mode命令可以找到其所有因子。
然后利用for循环得到因子之和，最后用if语句使之与原数相等而导出结果。
(3)假设:由于暂时没有想到得到真因子的方法，我们忽略了“真”字，于是得出下列程序：for a=1:1000;c=0;for b=1:1000;if mod(a,b)==0;c=c+b;endendif c==aaendend结果只得了一个a=1.于是我们改换思路，将以上程序中a除以b的商当做真因子并作了部分修改，得到了结果。
在做亲和数实验时，为了将完美数从结果中去掉，我们添加了变量之间的不等命令。
从中我们得到了启示，对以上程序添加了a>b命令，再次得到结果。
(4)问题求解：方法一for a=1:1000;c=0;for b=1:1000;if mod(a,b)==0&a>b;c=c+b;endendif c==aaendend方法二for a=1:1000;s=0;for b=2:1000;if mod(a,b)==0;c=a/b;s=s+c;endendif s==asendend(5)结论及分析：通过对自己原来固有思维的大胆突破，我们得出了结果；通过对原来程序的修改，我们再次得到结果。
修改后的程序运行较快，结果正确。
(6)结果：一千以内的完美数有：6 ,28, 496(7)结论推广：程序中的a为1到1000，而a的取值可以时任意区间的自然数，可以通过修改a的范围得到其他结果。
4.亲和数(1)题目:亲和数又叫友好数，是指两个自然数，其中每个自然数的真因数之和等于另一个数。
通过编程计算，你能得到多少对亲和数。
(2)分析问题:此题实际上是在完美数的基础上的一个拓展，故其思想与完美数基本一样。
(3)假设:我们假设亲和数可以是两个相同的完美数，于是将完美数的方法二分别用于两类数上，并进行了适当的修改，得出以下程序：for a=1:1000;x=0;for b=2:1000;if mod(a,b)==0;c=a/b;x=x+c;endendfor l=1:1000y=0;for m=2:1000;if mod(l,m)==0;n=l/m;y=y+n;endendif y==a&l==x;[a,l]endendend此程序有两大缺点：1、运行速度慢；2、没有将完美数排除在外。
于是我们又做了些优化,将以上程序中的x直接取代l，于是得出了运行较快的程序，而这样无疑是将运算范围缩小了一半。
此外，我们添加了a~=l命令，解决了第二个缺陷。
同时这个命令使我们得到了启示，于是得出了完美数的方法一，也随即得出了亲和数的第三、四种方法。
(4)问题求解：方法一for a=1:10000;x=0;for b=2:10000;if mod(a,b)==0; c=a/b;x=x+c;endendfor l=1:10000y=0;for m=2:10000;if mod(l,m)==0; n=l/m;y=y+n;endendif y==a&l==x&a~=l; [a,l]endendend方法二for a=1:20000;x=0;y=0;for b=2:20000;if mod(a,b)==0; c=a/b;x=x+c;endendfor m=2:10000;if mod(x,m)==0; n=x/m;y=y+n;endendif y==a&a~=x;[a,x]endend方法三for a=1:10000;c=0;for b=1:10000;if mod(a,b)==0&a>b; c=c+b;endendfor x=1:10000;z=0;for y=1:10000;if mod(x,y)==0&x>y; z=z+y;endendif c==x&z==a&x~=a;[a,x]endendend方法四for a=1:10000;c=0;z=0;for b=1:10000;if mod(a,b)==0&a>b; c=c+b;endendfor y=1:10000;if mod(c,y)==0&c>y; z=z+y;endendif a==z&a~=c[a,c]endend(5)结论及分析：虽然4种方法都可得出结果，但通过优化的方法，明显缩短了运算时间。
(6)结果：我们得出了10000以内的所有亲和数：220 284284 2201184 12101210 11842620 29242924 26205020 55645564 50206232 63686368 6232(7)结论推广：程序中的运算范围我们定为1到10000，，但通过改变对变量的赋值范围，就可以得出其他更多的结果实验感想本次实验我们初步了解了MATLAB的使用方法，尤其对简单编程有了进一步的实践，并掌握了一些基本操作方法和编程思想。
通过实验，电脑软件的强大功能让我们认识到了其对数学研究的重要性和必要性。
由于第一次接触MATLAB，在实验过程中我们遇到了不少困难，但我们通过查阅书本、互相讨论、反复研究以及向学长、老师请教等途径将它们一一克服。
尤其是我们的团结协作，在此体现得非常明显，很多问题都是在大家不断的商量与完善中才得出的结果。
而对于同一个的问题，每个人通常都会有不同的看法，在不断的对比与优化中我们才逐渐得到比较好的程序。
另外在实验报告的完成过程中，大家各尽所长。
这次试验使我们更加清楚的认识到团结合作的重要性。
在今后的学习中，我们要更认真地学习MATLAB，使之成为我们学习数学和日后进行研究的好帮手；同时，我们要更加重视团队合作，取长补短，精益求精。
2009年12月23日。