不知是我想复杂了,还是你说的比较简单
一般f表述函数,你指的f是一个数吧,如果是别的,可以继续探讨
你试试这样行不行:
clear
all;clc;
f=input('please
input
f:');
n=(1/2)f
你这个问题似乎没说太清楚。
目标函数是(x-10)^2+(y-7)^2,这句应该是指“终点是(10,7)”吧?目标函数和终点应该并不是一个意思。
约束条件是-((x-5)^2+(y-5)^2)小于等于-4,这句应该是指“避开圆心在(5,5),半径为2的圆”吧?
步长02,精度0001,这个是指每一步移动的长度是固定的02(±0001)?
现在最大的问题是,怎样评价路径是不是最佳——按最短吗?你的题中并没有说明。
我们学校的数学建模上机课也有Mathlab程序,看看下面有没有你要找的。
一 基本运算
1 求
输入(12+2(7-4))/3^2执行
2 输入x = (52+13-08)10^2/25执行
再输入y= 2x+1执行
3 执行clear命令。观察结果
4计算圆面积Area = ,半径r = 2,则可键入
r=2;area=pir^2; area
问:语句末尾加分号与不加分号有何区别?请试验之
5常用函数
名称 含义 名称 含义
sin 正弦 exp E为底的指数
cos 余弦 log 自然对数
tan 正切 log10 10为底的对数
cot 余切 log2 2为底的对数
asin 反正弦 abs 绝对值
acos 反余弦
例:1)执行y = sin(10)exp(-034^2)
2) 想计算 的值
输入y1=2sin(03pi)/(1+sqrt(5))执行之
若又想计算 ,可以简便地用 *** 作:先按á键则会出现上面输入过的指令 y1=2sin(03pi)/(1+sqrt(5)) ;然后移动光标,把y1改成y2;把 sin 改成 cos 便可。即得
y2=2cos(03pi)/(1+sqrt(5))然后执行之。
系统默认4位有效数字,若想提高精度则可如下:
digits(10);sym(y2,'d') 执行就可精确到小数点后10位,还可将10改为其它数字试验
二 矩阵运算
1要得到矩阵 ,
可输入A = [1,2,3; 4,5,6; 7,8,9] 执行,观察结果
还可分行输入
A=[1,2,3
4,5,6
7,8,9]
效果相同
2 注意 %号后的语句为注释,练习时不必输入
>>a=[1,4,6,8,10] %一维矩阵
>>a(3) % a的第三个元素
ans =
6
»x =[1 2 3 4 5 6 7 8
4 5 6 7 8 9 10 11]; %二维2x8 矩阵
执行后双击左边Workspace里的x,观察之
» x(3) % x的第三个元素
ans =
2
» x([1 2 5]) % x的第一、二、五个元素
ans =
1 4 3
如需要还可定义b=x([1 2 5])执行后结果为
b =
1 4 3
>> x(2,3) % x的第二行第三列的元素
ans =
6
x(1:5) % x的第前五个元素
ans =
1 4 2 5 3
» x(10:end) % x的第十个元素后的元素
ans =
8 6 9 7 10 8 11
执行后双击左边Workspace里的x,观察是哪十个元素
» x(10:-1:2) % x的第十个元素和第二个元素的倒排
ans =
8 5 7 4 6 3 5 2 4
» x(find(x>5)) % x中大于5的元素
ans =
6 7 8 6 9 7 10 8 11
» x(4)=100 %给x的第四个元素重新给值
x =
1 2 3 4 5 6 7 8
4 100 6 7 8 9 10 11
» x(3)=[] % 删除第三个元素(不是二维数组)
x =
Columns 1 through 12
1 4 100 3 6 4 7 5 8 6 9 7
Columns 13 through 15
10 8 11
» x(16)=1 % 加入第十六个元素
x =
Columns 1 through 12
1 4 100 3 6 4 7 5 8 6 9 7
Columns 13 through 16
10 8 11 1
3 如不需要以前的变量时,为不干扰以后计算,可执行clear清除以前的变量
当元素很多的时候,则须采用以下的方式:
» x=(1:2:121); % 以起始值为1,增量值为2,终止值为121的矩阵
» x=linspace(0,1,100); % 利用linspace,生成以0为起始值,1为终止值,元素数目为100的矩阵
»a=[] %空矩阵
a =
[]
» zeros(2,2) %全为0的矩阵
ans =
0 0
0 0
» ones(3,3) %全为1的矩阵
ans =
1 1 1
1 1 1
1 1 1
» rand(2,4); %随机矩阵
4另外一种定义矩阵的方式
»a=1:7; b=11:2:23;
»c=[b a]; %利用上面建立的阵列 a 及阵列 b ,组成新阵列c
»d=[b ; a]; %利用a及b,组成新矩阵d
执行后双击左边Workspace里的c与d,比较之
再如 已知y=[-1,6,15,7,31,2,4,5];
x=y(3:5) %x为y的第三到第五个元素组成的新向量
或 x=[y(5),y(3),y(7)] %x为y的第五、第三、第七个元素组成的新向量
或这样更简单 x=y([5,3,7])
5 输入矩阵x=[4,8,12,10,23;6,3,15,13,19;9,1,2,18,14;11,7,5,21,17]
依次输入下列命令并执行,观察结果,各命令分别有什么作用?
max(x)
min(x) (问:如何得到整个矩阵的最小值与最大值?)
[m,n]=size(x)
L=length(x)
y=x’
a=x( :,2)
b=x( :,2)’
c=x(3, :)
d=x(1 :3,3 :5)
y(2,3)=y(2,3)/2
y(2, :)=y(2, :)/2
y( :,4)=y( :,4)+y( :,2)
6 点运算 执行下列命令,指出点运算的作用
x=1 :8 (或对另外的向量或矩阵来作)
y=2^x
z=x/y
w=x^2
u=sin(x)
常用命令
min 最小值 max 最大值
mean 平均值 std 标准差
sort 排序 diff 相邻元素的差
length 个数 sum 总和
dot 内积 cross 外积
三 画图
二维图形
命 令 含 义 plot绘图函数的叁数
plot 建立向量或矩阵各队队向量的图形 字元 颜色 字元 图线型态
loglog x、y轴都取对数标度建立图形 y ** 点
semilogx x轴用于对数标度,y轴线性标度绘制图形 k 黑色 o 圆
semilogy y轴用于对数标度,x轴线性标度绘制图形 w 白色 x x
title 给图形加标题 b 蓝色 + +
xlabel 给x轴加标记 g 绿色
ylabel 给y轴加标记 r 红色 - 实线
text 在图形指定的位置上加文本字符串 c 亮青色 : 点线
gtext 在鼠标的位置上加文本字符串 m 锰紫色 - 点虚线
grid 打开网格线 -- 虚线
hold on 命令用于在已画好的图形上添加新的图形
1 x=0:0001:10; % 0到10的1000个点(每隔0001画一个点)的x座标
y=sin(x); % 对应的y座标
plot(x,y); % 绘图
注:matlab画图实际上就是描点连线,因此如果点取得不密,画出来就成了折线图,请试验之
2 Y=sin(10x);
plot(x,y,'r:',x,Y,'b') % 同时画两个函数
3 若要改变颜色,在座标对后面加上相关字串即可:
x=0:001:10;
plot(x,sin(x),'r')
4 若要同时改变颜色及图线型态(Line style),也是在坐标对后面加上相关字串即可:
plot(x,sin(x),'r')
5 用axis([xmin,xmax,ymin,ymax])函数来调整图轴的范围
axis([0,6,-15,1])
6 MATLAB也可对图形加上各种注解与处理:(见上表)
xlabel('x轴'); % x轴注解
ylabel('y轴'); % y轴注解
title('余弦函数'); % 图形标题
legend('y = cos(x)'); % 图形注解
gtext('y = cos(x)'); % 图形注解 ,用鼠标定位注解位置
grid on; % 显示格线
7画椭圆
a = [0:pi/50:2pi]'; %角度
X = cos(a)3; %参数方程
Y = sin(a)2;
plot(X,Y);
xlabel('x'), ylabel('y');
title('椭圆')
8 绘制函数 在0 ≤ x ≤ 1时的曲线。
x=0:01:1
y=xexp(-x) %为什么用点运算?若不用会怎样
plot(x,y),xlabel('x'),ylabel('y'),title('y=xexp(-x)')
9 画出衰减振荡曲线 与它的包络线 及 。t 的取值范围是[0, 4π] 。
t=0:pi/50:4pi;
y0=exp(-t/3);
y=exp(-t/3)sin(3t);
plot(t,y,'-r',t,y0,':b',t,-y0,':b') % -r表示红色实线,:b表示蓝色点线,看上表
grid
10 在同一个画面上建立几个坐标系, 用subplot(m,n,p)命令;把一个画面分成m×n个图形区域, p代表当前的区域号,在每个区域中分别画一个图,如
x=linspace(0,2pi,30); y=sin(x); z=cos(x);
u=2sin(x)cos(x); v=sin(x)/cos(x);
subplot(2,2,1),plot(x,y),axis([0 2pi -1 1]),title('sin(x)')
subplot(2,2,2),plot(x,z),axis([0 2pi -1 1]),title('cos(x)')
subplot(2,2,3),plot(x,u),axis([0 2pi -1 1]),title('2sin(x)cos(x)')
subplot(2,2,4),plot(x,v),axis([0 2pi -20 20]),title('sin(x)/cos(x)')
三维图形
11三维螺旋线:
t=0:pi/50:10pi;
plot3(sin(t),cos(t),t) %参数方程
grid %添加网格
12 t=linspace(0,20pi, 501);
plot3(tsin(t), tcos(t), t); %注意点乘
也可以同时画出两条曲线,格式与二维情况类似,兹不举例。
13用mesh命令画曲面
画出由函数 形成的立体网状图:
a=linspace(-2, 2, 25); % 在x轴上从(-2,2)取25点
b=linspace(-2, 2, 25); % 在y轴上取25点
[x,y]=meshgrid(a, b); % x和y都是21x21的矩阵
z=xexp(-x^2-y^2); % 计算函数值,z也是21x21的矩阵
mesh(x, y, z); % 画出立体网状图
14 surf和mesh的用法类似:
a=linspace(-2, 2, 25); % 在x轴上取25点
b=linspace(-2, 2, 25); % 在y轴上取25点
[x,y]=meshgrid(a, b); % x和y都是21x21的矩阵
z=xexp(-x^2-y^2); % 计算函数值,z也是21x21的矩阵
surf(x, y, z); % 画出立体曲面图
四 程序设计
1 M-文件: 上面所做的运算都是在命令窗口中输入一条或两三条命令,然后执行,再输入,再执行,以这样交谈式的方式进行。如果为了解决某一问题需要很多命令,这样做就很不方便了。这时我们把解决某一问题的所有命令集中放在一个文档里,命名、保存。然后只要在命令窗口中输入文档名,执行即可。
例:(1)编写文档:点击MATLAB指令窗口上面最左端的图标 ,即新建文件,就可打开MATLAB文件编辑器。用户即可在空白窗口中编写程序。例如输入下面的程序:
x=linspace(0,2pi,20);
y=sin(x);
plot(x,y,'r+')
title('2D plot')
(2)点击文件编辑器上面工具条中的保存 ,命名(例如将上面的程序命名为picture),然后保存。像这样在MATLAB文件编辑器中编写的文件叫M-文件(M-file)。
(3)运行:i)在命令窗口中输入文件名(如上面的picture),然后执行。
ii)或直接在文件编辑器上面的工具条中找到debug(即调试),点击,再找到run(即运行),再点击即可。
同学们可以把前面画图的一些问题放在文件编辑器里再做一下。
2 自己编写函数:我们经常用到的像sin、cos、exp这样的一些函数都是MATLAB软件自身所带的函数,因此直接应用即可,但有时我们为了解决一些问题需要自己编写函数。自己编写函数有两个基本要求i)必须在MATLAB文件编辑器中编写。ii)函数名和文件名必须相同。 例: 编写函数 , 计算f(1)f(2)+f2(3)
(1)打开MATLAB文件编辑器,输入
function Y= fun1(x) % 表示Y是x的函数,x是自变量,fun1是函数名
Y=(x^3 - 2x^2 + x - 63)/(x^2 + 005x - 314);
然后保存。
注:在自己编写的函数前都要写上function,表示这是自己定义的函数。fun1表示函数名,那么最后文件名也应命名为fun1。
(2)这样在命令窗口中就可以像应用sin、cos那样来使用函数fun1,如:在命令窗口中输入 >> fun1(1)fun1(2)+fun1(3)fun1(3) 结果为:
ans =
-126023
3 for循环语句(这里的for语句与C语言中的for语句不同,要更简单一些)
例:一个简单的for循环示例。
for i=1:10; % i依次取1,2,…10,
x(i)=2i; % 对每个i值,重复执行该指令
end; % 表示循环结束,每一个for要对应一个end
x % 要求显示运行后数组x的值。
输入后观察结果,体会for语句的作用。
注:在MATLAB里(在C语言中也一样), 的作用表示把等号右边的值送给左边的变量,这和数学中相等的意思不同。下面的例子中都要这样理解,否则就不能明白程序的含义。
4 while循环语句
例: Fibonacci 数列:1,1,2,3,5,8,… 即: ,( 1,2,3…)现要求该数列中第一个大于10000 的元素。
a(1)=1;a(2)=1;i=2;
while a(i)<=10000
a(i+1)=a(i-1)+a(i);
i=i+1;
end;
i,a(i),
5(1)if-end语句,例:
cost=10;number=12;
if number>8
sums=number095cost;
end,
sums
(2)if-else-end语句,例:
cost=10;number=5; % 改变number的初值,看结果有何不同
if number>8
sums=number095cost;
else sums=number05cost;
end,
sums
6 例:用for 循环语句来寻找Fibonacc 数列中第一个大于10000 的元素。
n=100;a=ones(1,n); % a是一个一行,n列的所有元素为1的矩阵
for i=3:n
a(i)=a(i-1)+a(i-2);
if a(i)>=10000
a(i),
break; % 表示跳出循环
end;
end, i
7 练习:课本264页,参考例4右边的流程图114,编程序求解例4,自己设置误差,并与书上的结果比较。
五 拟合与插值
曲 线 拟 合 和 插 值 函 数
polyfit(x, y, n) 对描述n阶多项式y=f(x)的数据进行最小二乘曲线拟合
interp1(x, y, xo) 1维线性插值
interp1(x, y, xo, ' spline ') 1维3次样条插值
interp1(x, y, xo, ' cubic ') 1维3次插值
interp2(x, y, Z, xi, yi) 2维线性插值
interp2(x, y, Z, xi, yi, ' cubic') 2维3次插值
1 插值
看课本266页§112第一段,了解什么是插值。
例:考虑下列问题,12小时内,一小时测量一次室外温度。数据如下:
时间:1,2,3, 4, 5, 6, 7, 8, 9,10,11,12
温度:5,8,9,15,25,29,31,30,22,25,27,24
现在根据以上数据估计32,47等时刻的温度
hours=1:12;
temps=[5 8 9 15 25 29 31 30 22 25 27 24];
t=interp1(hours, temps, [32,47]) % 一阶线性插值,如果只估计一个点的值,则无须加方括号
改为t=interp1(hours, temps, [32,47], 'spline') 则为三次样条插值
如果输入如下程序,则画出插值曲线
hours=1:12;
temps=[5 8 9 15 25 29 31 30 22 25 27 24];
h=1:01:12;
t=interp1(hours, temps, h) ; % h后加上'spline'则为三次样条插值
plot(hours, temps, ' + ' , h, t)
用一阶线性插值和三次样条插值做课本268页例2,与书上之结果比较,然后挑课后题做一两道。
2 拟合
看课本270页§113,曲线拟合,比较拟合与插值有什么区别。
例:两组数据如下:
x=[0 1 2 3 4 5 6 7 8 9 1];
y=[-447 1978 328 616 708 734 766 956 948 930 112];
n=8;
p1=polyfit(x,y,n); % n表示用n阶多项式拟合,n=1为线性拟合,即通常所说最小二乘法
poly2sym(p1) % 前面的拟合命令只给出多项式的系数,用此命令则将结果转化为真正的多项式。或用 vpa(poly2sym(p1),10) 即取数值形式,取10位有效数字
x1=0:01:1; % 由此以后三句是画出拟合曲线的图像
y1=polyval(p1,x1); %此句是在x1这些点处求出多项式的值,送给y1
plot(x,y,'o',x1,y1)
改变n的数字,即用不同的多项式拟合,看看哪个结果好。
当n=10时,数据点之间出现大的波动。当企图进行高阶曲线拟合时,这种波动现象经常发生,并不利于我们认识两组数据之间的规律,因此并不是阶数越高越好,实际问题当中,适当选一个即可。
用上面的指令做课本271页例1及例2,将结果与书上之结果比较。
if u<=1是对u这个矩阵在进行判断,所以其判定结果是0,直接进行下一步else运算,及对u中的元素开方。
对一个矩阵中的具体元素做判断,不能用if u<=1来进行
%%%%没有用到内部函数!!!
PrimeNumber=[2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97];
k=1;
for i=2:100
temp=i(i-1)-1;
for j=1:length(PrimeNumber)
if temp==PrimeNumber(j)
Result(k,1)=i; %亲密数对中的第一个数
Result(k,2)=i-1; %亲密数对中的第二个数
Result(k,3)=temp; %亲密素数
k=k+1;
break;
end
end
end
% 计算亲密素数和
SumPrime=0;
for i=1:k-1
SumPrime=SumPrime+Result(i,3);
end
% 数对个数
disp('数对个数');disp(k-1);
% 亲密素数和
disp('亲密素数和');disp(SumPrime);
1修改了一下,变动不大。
function x=Gauss_(A,b)
A=[10,-7,0;5,-1,5;-3,2,6];
b=[7;6;4]; %%%%%应该是列矩阵啊
n=length(b);
for k=1:n-1 %……消去过程……
if A(k,k)==0
disp('错误信息')
break;
end
index = [k+1:n];%用向量比用循环好
m = -A(index,k)/A(k,k); %提前定义一个较短的变量,相当于你的A(i,k)=A(i,k)/A(k,k);
A(index,index) = A(index,index) + mA(k,index);%矩阵A消元
b(index) = b(index) + mb(k);%矩阵b消元
end
x = zeros(n,1);%提前预置一个变量的位置,以提高运算速度
x(n) = b(n)/A(n,n);%……回代过程……
for i = n-1:-1:1
x(i) = ( b(i) - A(i,[i+1:n])x([i+1:n]) )/A(i,i);%还是写在一起吧,比较简捷
end
x
%%%%%%%%%结束啦
运行结果:
x =
00000
-10000
10000
2验证
>> A=[10,-7,0;5,-1,5;-3,2,6];
>> b=[7;6;4];
>> X=A\b
X =
0
-1
1
3验证
>> Gauss_(A,b)
x =
00000
-10000
10000
以上就是关于如何用matlab编一段小程序 实现公式1/2*f=n 输入一个f 可以得到n全部的内容,包括:如何用matlab编一段小程序 实现公式1/2*f=n 输入一个f 可以得到n、matlab编个小程序、帮忙写个Matlab小程序等相关内容解答,如果想了解更多相关内容,可以关注我们,你们的支持是我们更新的动力!
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)