MATLAB基础示例

matlab中函数文件的语法规范

1、函数文件的名称即为函数名称，文件名命名规则与变量命名规则一样：

a.函数名必须以字母开头；

b.函数名区分大小写；

c.函数名可由字母、数字、下划线组成。

2、文件内第一行通过关键词function来指明函数文件的格式如下：

function[输出参数列表]=函数名(输入参数列表)

不同参数之间用逗号分隔。注意：

a.函数的文件名一般与程序中定义的函数名相同；二者不一致时，文件名为函数名。

b.函数名要有一定的实际含义。

if 示例

nrows = 4;
ncols = 6;
A = ones(nrows,ncols);
for c = 1:ncols
    for r = 1:nrows
        if r == c
            A(r,c)=2;
        elseif abs(r-c) == 1
            A(r,c) = -1;
        else
            A(r,c) = 0;
        end
    end
end
A

switch 示例

n = input('Enter a number:');
switch n
    case -1
        disp('negative one')
    case 0
        disp('zero')
    case 1
        disp('positive one')
    otherwise
        disp('other value')
end

for 示例：创建一个10阶Hilbert矩阵

s = 10;
H = zeros(s);
for c = 1:s
    for r = 1:s
        H(r,c) = 1/(r+c-1)
    end
end

while 示例：计算n！

n = 10;
f = n;
while n>1
    n = n-1;
    f = f*n;
end
disp(f)

break 示例：求随机数序列之和，直到下一随机数大于上限为止

limit = 0.8;
s = 0;
while 1
    tmp = rand;
    if tmp>limit
        break
    end
    s = s+tmp;
end
disp(s)

break 示例：统计文件magic.m中的代码行数，使用continue语句跳过空白行和注释

fid = fopen('magic.m','r');
count = 0;
while~feof(fid)
    line = fgetl(fid);
    if isempty(line)||strncmp(line,'%',1)||~ischar(line)
        continue
    end
    count = count+1;
end
count

nargin,nargout的含义及用法

%nargin返回在调用当前执行的函数时指定的输入参数个数。仅在函数体中使用此nargin语法。
function c = addme(a,b)
%addme(a) 调用本函数时，nargin被赋值为1
%addme(a,b) 调用本函数时，nargin被赋值为2
switch nargin
    case 2
        c = a + b
    case 1
        c = a + a
    otherwise
        c = 0;
end

nargout返回在调用当前执行的函数时指定的输出参数个数。仅在函数体中使用此nargout语法。

function[dif,absdif] = subtract(y,x)
%dif=subtract(y,x) 调用本函数时，nargout被赋值为1
%[dif,absdif]=subtract(y,x) 调用本函数时，natgour被赋值为2
dif = y - x;
if nargout > 1
    disp('Calculating absolute value')
    absdif = abs(dif);
end

求矩阵的最小值

%方法一
x=[1,2,3;4,5,6;7,8,9];
min(min(x))
%方法二
x=[1,2,3;4,5,6;7,8,9];
min(x(:))

判断是否为闰年

%方法一
function flg = isLeapYear(y)
%isLeapYear 判断是否为闰年
%如果y能被4整除且不能被100整除，或者能被400整除，则y为闰年
if mod(y,4) == 0 & mod(y,100) ~= 0 | mod(y,400) == 0
    flg = 1;
else
    flg = 0;
end

%方法二
function flg = isLeapYear(y)
%isLeapYear 判断是否为闰年
%如果y能被4整除且不能被100整除，或者能被400整除，则y为闰年
flg = mod(y,4) == 0 & ...
    mod(y,100) ~= 0 |...
    mod(y,400) == 0
end

创建一个5行4列的矩阵，使矩阵元素m（i，j）=100*i+j

n1 = 5;
n2 = 4;
for i = 1:n1
    for j = 1:n2
    M(i,j) = 100*i+j;
    disp(M);
    end
end

编程计算斐波那契数列的元素并输出前n项

function fibo = fiboSeq(n)
fibo = [1,1];
for k = 3:n
    fibo(k) = fibo(k-1)+fibo(k-2);
end
end

提取一个整数的各位数，存到一个变量里

function r = getDigits(n)
r = mod(n,10);
r = g;
n = (n-r)/10;
while n ~= 0
    g = mod(n,10);
    r = [g,r];
    n = (n-g)/10;
end

求方程x^2-2*y=1的全部正整数解，并保存在矩阵中，输出正整数解的个数

%finds solutions
% x^2 - 2*y = 1 (x,y>=1,y<=1000)
s = [];
for y = 1:1000
    x = sqrt(1+2*y);
    if fix(x) == x
        s = [s;[x,y]];
    end
end
s

JiaoGu猜想

function JiaoGu = JiaoGu(n)
s = [];
while n ~= 1
    if mod(n,2) == 0
        n = n/2;
    else
        n = n*3+1;
    end
    s = [s,n];
end
disp(s)

编程输出某范围内的素数

function p = p(n)
p=2;
for k = 3:2:n;
    %if k is prime number
    %then p = [p,k]
    if isPrime(k)
        p = [p,k]
    end
end

function flg = isPrime(k)
for f = 3:2:sqrt(k)
    if mod(k,f) == 0
        flg = 0;
        return
    end
end
flg = 1;

编写程序做如下判定：如果a>3且b>3则显示字符串‘handle 1’，否则显示‘handle 2’

a = 3;b = 7;
if a>3 & b>3,
    disp('handle 1')
else
    disp('handle 2')
end

编写程序对一个百分制分数打分，要求不低于90分的成绩显示“成绩优异”；对于低于90分且不低于80分的成绩，显示“成绩优秀“；对于低于80分且不低于60分的成绩，显示成绩中等，其他则显示”未及格“。

grade = input('please input your score');
if grade >= 90
    disp('成绩优异')
elseif grade >= 80 & grade < 90
    disp('成绩优秀')
elseif grade >= 60 & grade < 80
    disp('成绩中等')
else
    disp('未及格')
end

编写程序使用switch语句完成百分制打分功能

grade = input('please input your score')
grade = fix(grade/10)
switch grade
    case {9,10},
        disp('成绩优异')
    case {8},
        disp('成绩优秀')
    case {6,7},
        disp{'成绩中等'}
    otherwise,
        disp('未及格')
end

已知分段函数f(x)定义如下，输入自变量x的值，然后输出函数值，

f(x)=(x+1)^2 (x<=0) and f(x)=(x-1)^2 (x>0)

x = input('please input a number')
if x > 0,
    disp((x-1)^2)
elseif x <= 0,
    disp((x+1)^2)
end

编写程序构造一个5行4列的矩阵

for i = 1:5
    for j = 1:4
        m(i,j) = i*100+j;
    end
end
disp(m)

编写程序将输入的字符串反序

str = input('please input a character:','s')
tmpstr = str;
i = 1;
len = length(str);
while i <= len
    str(len-i+1) = tmpstr(i)
    i=i+1;
end
str

设银行的年利率为4.25%，将10000元存入银行，问多长时间会连本带利翻一番

a0 = 10000;
a = a0;
r = 0.0425;
for i = 1:200
    a = a*(1+r);
    if a >= 2*a0
        disp(sprintf('存了%d年终于翻番了',i))
        break;
    end
end

现有一个求解线性方程组的迭代解法，对于对角元为0的情形无法处理，残差向量模较小时推出迭代。

下面的程序只给出了迭代法的框架，没有给出具体的迭代计算式：

%迭代法
%1.初始化
n = 100;
maxit = 300;%最大迭代次数
x = rand(n,1);
%2.主循环
for i = 1:maxit
    x = x + myexp;%以myexp表示具体的计算式
    r = A*x + b;
    d = A(i,i);
    if d == 0
        error('我不能处理对角元为0的问题，不要逼我！')
    end
    if norm(r) < 1e-6 %收敛条件
        break;
    end
end

从100-10000的整数中，从小到大搜索，找出第3个、第5个除以12余1，且除以25余2的整数。

count = 0;
m = [];%存储满足条件的整数，初始化为空矩阵
for i = 100:10000
    if mod(i,12) == 1 && mod(i,25) == 2
        count = count + 1;
        if count == 3 || count == 5
            m=[m,i];
            if count >=5 %找到了，退出该循环
                break
            end
        end
    end
end
m

用二分法求方程exp(-x)-sin(pi*x/2)=0在区间[0,1]内的根，要求误差不超过2^-5

f = inline('exp(-x) - sin(pi*x/2)');
a = 0; b = 1; n = 0; er0 = 2^(-5);
fa = f(a); fb = f(b); er = b-a;
while er>er0
    c = (a+b)/2;
    fc = f(c);
    if fa*fc<0
        b = c;
    else
        a = c;
        fa = fc;
    end
    disp([a,b]);
    er = b - a;
    k = k + 1
end

用迭代法求方程x^3-x-1=0在区间 [1, 2]内的根，要求误差不超过10^-7。

f = inline('(x+1)^(1/3)');
x0 = 1.5; a = 1; b = 2; k = 0; er0 = 10^-7;
er = b-a;
while er>er0
    x = f(x0);
    er = abs(x - x0);
    k = k + 1;
    x0 = x;
end
x0,k


f1 = inline('(x+1)^(1/3)');
x0 = 1.5; a = 1; b = 5; 
ezplot(f1,[1,2]),grid on,hold on;
f = inline('x^3 - x - 1');
ezplot(f,[1,2]),grid on,hold on;
k = 0; er0 = 10^-7;
er = b-a;
while er>er0
    x = f1(x0);
    er = abs(x - x0);
    k = k + 1;
    x0 = x;
    plot(x,x0,'o'),hold on;
end
x0,k

f(x)=x*e^x-1，用牛顿迭代法求f(x)=0在x0=0.5附近的根。

f = inline('x*exp(x) - 1');
df = inline('exp(x) + x*exp(x)');
x0 = 0.5;
x1 = x0 - f(x0)/df(x0);
er0 = 10^(-7);
k = 0;
while abs(x1-x0)>er0
    x0 = x1;
    x1 = x0 - f(x0)/df(x0);
    k = k + 1;
end
x1,k

用牛顿迭代法求x^3-2*sin(x)=0的根。

f = inline('x^3 - 2*sin(x)');
df = inline('3*x^2 - 2*cos(x)');
x0 = input('x0 = ');
x1 = x0 - f(x0)/df(x0);
er0 = 10^(-7);
k = 0;
while abs(x1-x0)>er0
    x0 = x1;
    x1 = x0 - f(x0)/df(x0);
    k = k + 1;
end
x1,k

绘制地图网络

function earth
t = linspace(-pi/2,pi/2,20);
theta = linspace(0.2*pi,50);
[t1,t2] = meshgrid(t,theta);
x = cos(t1).*cos(t2);
y = cos(t1).*sin(t2);
z = sin(t1);
mesh(x,y,z)
colormap([0,0,1])
axis square off

旋转正方形

theta = 15*pi/180;
ratio = 1/(cos(theta)+sin(theta));
x = [-1,1,1,-1,-1];
y = [-1,-1,1,1,-1];
plot(x,y,'b'),hold on
axis equal off %去掉坐标轴
for k=1:20
    pause(0.5)
    z = x + y*i;
    z = z*ratio*exp(i*theta);
    x = real(z);
    y = imag(z);
    plot(x,y,'b');

旋转抛物面

%z = x^2 - y^2
t = linspace(0,2*pi,100);
r = linspace(0,4,50);
[T,R] = meshgrid(t,r);
X = R.*cos(T);
Y = R.*sin(T);
Z = X.^2 +Y.^2;
surf(X,Y,Z)

求y(x) = a*x/sqrt(a^2 + x^2）在5a处的二阶导数值

function r = myfun
syms a x y(x)
y(x) = a*x/sqrt(a^2 + x^2);
d2y = diff(y,x,2);
r = d2y(5*a);
disp(r)

非线性方程求根：fzero

f1 = @(x) (x + 2)^x - 2;
fplot(f1,[-1,2])
x0 = [-1,2];
r2 = fzero(f1,x0) %在x0附近搜索根

蒙特卡罗方法随机投点

f = @(x) x.^3 - x - 1;
r1 = roots([1 0 -1 -1])
x = rand(1,1e+6)*4 - 1;
y = abs(f(x));
id = find(y == min(y));
x(id)
r2 = x(id)

不动点迭代法的实现

%方法1
function [root,n] = fixpoint(phai,x0,tol)
if (nargin == 2)
    tol = 1.0e-5;
end
err = 1;
root = x0;
n = 0;
while(err>tol)
    n = n + 1; %迭代次数
    r1 = root;
    root = phai(r1); %计算函数值
    err = abs(root - r1);
end

%方法2

牛顿迭代法的实现

tol = 1e-6;
f = @(x) x.^3 - x - 1;
df = @(x) 3*x.^2 - 1;
x0 = 3.5;
y0 = f(x0);
x1 = x0 - y0/df(x0);
while abs(x0 - x1)>tol
    x0 = x1;
    y0 = f(x0);
    x1 = x0 - y0/df(x0);
end
x1

曲线拟合

%多项式拟合函数：polyfit
%调用格式：p = polyfit(x,y,n)
%说明：x，y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量p

%多项式求值函数：poltval
%调用格式：p = polyval(p,x)
%说明：y = polyval(p,x)为返回对应自变量x在给定系数p的多项式的值

x = [2 3 4 5 7 8 10 11 14 15 16 18 19];
y = [106.42 108.26 109.58 109.5 110 109.93 110.49 110.59 110.6 110.9 110.76 111 111.2];
figure,plot(x,y,'+')
v = polyfit(x,y,2); %多项式拟合
t = 1:0.1:19;
u = polyval(v,t); %多项式求值
hold on
plot(t,u,'r')

欧拉法

%dydt = t-2y,y(0) = 1
n = 100;
t = linspace(0,2,n);
h = t(2) - t(1);
y = zeros(1,n);
y(1) = 1;
dydt = @(t,y) t-2*y;
for k = 2:n
    y(k) = y(k - 1) + dydt(t(k - 1),y(k - 1))*h;
end
plot(t,y,'-')

Matlab求解微分方程

%用函数文件定义一阶微分方程右端函数
%用matlab命令ode23（）求数值解
%基本语法：[T,Y] = ode23(‘F’,Tspan，y0)
%'F'是包括函数文件名字的字符串或函数句柄
%Tspan = [t0,tN]是常微分方程的求解区域
%y0表示初始条件
%tspan可以具体指定节点

方法1：子函数法
%dydt = t-2y,y(0) = 1
function main
n = 200;
t = linspace(0,2,n);
%dydt = @(t,y)t-2*y;
[t,Y] = ode23(@dydt,t,1); % @dydt处可改为‘f’
plot(t,Y,'-')
end

function f = dydt(t,y)
f = t - 2*y;
end

方法2：
%马尔萨斯人口增长模型
%dN/dt = 0.015*N,N(1994) = 12
t = [1994 2021];
dNdt = @(t,N) 0.015*N;
[t,N] = ode23(dNdt,t,12);
plot(t,N,'-O')

ode求解常微分方程组

function L8fcz
y0 = [0,2];
tn = [0,3];
[t,y] = ode23(@fun2,tn,y0);
plot(t,y(:,1),'o',t,y(:,2),'+')
xlabel('t');
ylabel('y_1 & y_2');
legend('y_1','y_2')
end

function dydt = fun2(t,y)
dydt = [-y(1)*exp(1 - t) + 0.8*y(2);y(1) - y(2)^3]; %方程组
end

ode求解高阶常微分方程组

%对方程组进行降阶，将二阶变为一阶的导数

求解方程

syms x t
str = int(exp(-t^2),0,x) - 0.5;
f = @(x) (pi^(1/2)*erf(x))/2 - 1/2;
x1 = fzero(f,0)
x2 = fsolve(f,0)
x3 = double(solve(str))
x = linspace(0,1,1e+6);
y = abs(f(x));
id = find(y == min(y));
x4 = x(id)

加油站油库深度12米呈啤酒桶形状,顶部和底部半径为2米,中心半径为3米。纵截面右侧是抛物线 x=3−y.^2/36 (−6≤y≤6)，计算储油量

t = linspace(-6,6,30);
r = 3 - t.^2/36;
cylinder(r,30),axis off
colormap(’jet’)
syms y h
f = 3 - y*y/36;
V = pi*int(f*f,y,-6,h)
V12 = subs(V,h,6) %符号表达式中的变量替换方法——subs（）

数值积分

函数trapz的使用
%Q=trapz(Y) 通过梯形法计算Y的近似积分（采用单位间距）
%Q=trapz(X,Y) 根据X制定的坐标或标量间距对Y进行积分
Y = [1 4 9 16 25];
Q = trapz(Y)

函数quad的使用
%q=quad(fun,a,b) 尝试使用递归Simpson积分法求取函数fun从a到b的近似积分，误差小于1e-6，fun是函数句柄
%q=quad(fun,a,b,tol) 使用绝对误差容限tol代替默认值1e-6
F = @(x)1./(x.^3 - 2*x - 5);
Q = quad(F,0,2)

函数dblquad的使用
%矩形区域上的二重积分的数值计算
%q=dblquad(fun,xmin,xmax,ymin,ymax)
%q=dblquad(fun,xmin,xmax,ymin,ymax,tol) 使用容差tol代替默认值1e-6
F = @(x)1./(x.^3 - 2*x - 5);
Q = dblquad(F,1,2,0,1)

函数quad2d的使用
%计算二重数值积分——tiled方法
%q=quad2d(fun,a,b,c,d) 逼近fun(x,y)在平面区域a<=x<=b和c(x)<=y<=d(x)上的积分，边界c和d均可为标量和函数句柄

最优化

%线性规划问题
C = [500,600];
A = [-1,-1;1,0;0,1];
b = [-5000;5000;3000];
Aeq = [];
beq = [];
Lb = [0;0];
Ub = [inf,inf];
X = linprog(-C,A,b,Aeq,beq,Lb,Ub) %若改为[x,fval] = linprog(-C,A,b,Aeq.beq,Lb,Ub)
Z = C*X                           %则此时Z = -fval

C = [500,600];
A = [-1,-1];
b = -5000;
Aeq = [];
beq = [];
Lb = [0;0];
Ub = [5000,3000];
X = linprog(-C,A,b,Aeq,beq,Lb,Ub)
Z = C*X

随机试验

%高尔顿板模拟
%50层钢柱，1e5个钢珠
n = 50;
numBall = 1e5;
x = rand(n,numBall);
x(x>0.5) = 1;
x(x<=0.5) = 0;
Pos = sum(x);
hist(Pos,n+1)

%用蒙特卡罗方法估算圆周率
%设圆的半径为1，则其外切正方形的面积为4，圆的面积为pi，现在模拟产生在正方形ABCD中均匀分布的点n个，如果这n个点中有m个点在该圆内，则圆的面积与正方形ABCD的面积之比可近似为m/n，即pi/4≈m/n。
n = 1e6;
x = unifrnd(-1,1,n,1); %x = 2*rand(1,n)-1;
y = unifrnd(-1,1,n,1); %y = 2*rand(1,n)-1;
id = find(x.^2+y.^2<=1);
m = length(id); % m = sum(x.^2+y.^2<=1)
mypi = 4*m/n

%甲乙两船在24小时内独立地随机到达码头，设两船到达码头时刻分别为X和Y，且X ~ U(0 , 24), Y ~ U(0 , 24)。如果甲船到达码头后停留2小时,乙船到达码头后停留1小时.问两船相遇的概率有多大？
N = 2000;
P = 24*rand(2,N);
X = P(1,:);
Y = P(2,:);
I = find(X<=Y&Y<=X+2);
J = find(Y<=X&X<=Y+1);
F = (length(I)+length(J))/N