Matlab UESTC 数学实验题目参考

 注：上传可能导致了代码中空格消失或多余，请自行斟酌哦~

数学实验常考点

1.矩阵的求逆inv、转置（B‘）以及简单的线性方程组求解（AX=b得到X=inv(A)*b）

2.简单的符号计算（求导diff、求极限limit以及求积分quad/int），符号计算字母的替换subs，double型的转换

3.fminserach、fminbnd、fzero、dsolve、solve、ode23等等的使用（这些函数一般会声明指定你使用）

4.简单的数组运算（数组行列的交换位置、删除、加减、拼接）

5.分段函数的书写（建议用我在文档里的条件判断的方式而不要用ifelse）

6.子函数的书写，知道怎么调用和书写

7.linprog用于线性规划（目标函数和约束关系均线性），fmincon用于非线性规划，ga不仅用于非线性规划，还可以设置某几个自变量恒为整数（三个函数的表达式考场可以用help看，但是要知道是什么意思，尤其注意非线性不等式和等式的子函数书写方式）

8.蒙特卡罗方法求解概率、面积的方法

9.polyfit和polyval的使用（记得之前有一道题建议你使用lsqcurvefit，结果死活运行不出来，一换成polyfit就可以运行了）

10.简单的for/while循环用于运算

11.简单的findstr/strfind/find的运用

12.简单的meshgrid、mesh的运用

13.简单的rand()、nornmd的使用

数学实验题目

function r=myfun
% x = fminsearch(fun,x0)
% 其中fun为函数，x0为初始值点 输出的是极小值点
f=@(x)6.*x.^2+10.*x+11+10.*sin(x);
r=fminsearch(f,5);

function A=myfun
u=[8 5 ; 11 12];
v=[7 2; 1 9];
w=[4 7; 6 1];
%元素或者矩阵的拼接：同一行用逗号或空格隔开，不同列用分号隔开
A=[u,v;w,u];

function s=myfun
%线代知识：已知AX=b，则X=A^-1*b
A=[1 1 1 1;1 1 -1 -1;1 -1 1 -1;1 -1 -1 1];
b=[2;4;2;2];
s=inv(A)*b;

function r=myfun
rand('seed',1523)
A=fix(20*rand(10,10));
c=fix(5*rand(1,10));
A=[A;c];%不同列拼接用分号
A(:,[1 3])=A(:,[3 1]);%交换
A(1:2,:)=[];%删除
r=A;

function [C,D,E]=myfun
rand('seed',1291)
B1=fix(20*rand(8,8));
B2=fix(20*rand(8,8));
%转置：A；即可；求逆：inv(A)
C=B1+B2';
D=B1*B2;
E=B1'-B2';

function y=myfun(x)
%本次采用条件判断的方式书写（也可以用if else书写）
y=(x>=8).*(6.*x.^3+9)+(x>=6).*(x<=8).*(4.*x+4)+(x<6).*(9.*x.^2+5);
%解释：比如若x=9，那么(x>=8)就会为真，结果为1，其余判断为0
%注意输入的x可能是向量、矩阵，需要使用点乘点除

function s=myfun
s=0;%初始值设定
for i=1:500
if mod(i,2)==0
s=s+i;
end
end

function x=myfun
x(1)=1;%初始值设定
for i=2:20
x(i)=2*x(i-1)+0.5;
end

function r=myfun
syms x
f=((4-2*x)/(6*x-10))^x;
%符号计算求极限：limit(函数，趋近于多少)
%inf是正无穷，负无穷-inf
r=limit(f,inf);

function [x,fval]=myfun
%题目中函数、约束条件都是线性的，故使用线性规划linprog
%help查询得格式：[x,fval] = linprog(f,A,b,Aeq,beq,lb,ub)
%其中f为函数的系数，A、b对应约束不等式，Aeq、beq对应约束等式
%lb和ub自变量的最大值和最小值
%注意：必须是求min（求max加负号），不等式符号是<=(若是>=加负号)
f=[-2 -1 -5 -6];
A=[2 0 1 1 ; 2 2 1 2]; b=[8;12];
Aeq=[]; beq=[];
lb=[0 0 0 0]; ub=[];
[x,fval] = linprog(f,A,b,Aeq,beq,lb,ub);

注：遇到ga，fmincon，linprog等输入参数多的，可使用help查询输入形式，不用死记

function bb=myfun
c=0.354274;
bb=0;%初始值设定
for i=3:30
bb=bb+c^i;
end

function e=myfun
rand('seed',1240)
u=fix(13*rand(1,35));
e=u(1,4:9);%赋值

function M=myfun
M=[];%初始化矩阵
%采用简单得for循环+条件判断解决
%由不等式可以知道x1~x3最大值为4，7，6
for x1=0:4
    for x2=0:7
        for x3=0:6
            if 4.8*x1+3*x2+3.3*x3<=22
            M=[M;x1,x2,x3];
            end
        end
    end
end

function B=myfun
rand('seed',2213)
B=fix(20*rand(9,8));
B=[B(:,1:2),B(:,3).^2,B(:,4:8)];     %注意考察元素点乘
B(:,[2 8])=B(:,[8 2]);                %2 8列互换位置
B=[B(:,1:6),B(:,7)+B(:,3),B(:,8)];
B=[B(1:7,:);zeros(1,8);B(9,:)];      %注意这里是第8行而不是第8列
%本题暂未考察元素的删除
%假设要删除第8行:B(8,:)=[];
%假设删除第3-5列:B(:,3:5)=[];
%假设删除第2、9行:B([2 9],:)=[];

function M=myfun
M=[];%初始化矩阵
for i=1:6
    m=[];%用于存储同一行的元素
    for j=1:8
        m=[m,i*j];
    end
    M=[M;m];
    m=[];%清空临时向量，用于下一次循环
end

function r=myfun
rand('seed',2477)
pp=97+fix(26*rand(1,31500));
pp=char(pp);
%查找k = strfind(str,pat) pat是要查找的部分
r=strfind(pp,'bab');

function v=testmain
x=-5:0.5:10;
v=fun(x);

function r=fun(x)
%采用条件判断的形式书写,也可用if else分段书写
%注意使用点乘点除，方便输入矩阵等形式的参数
r=(x+2<0).*(2.*x-3)+(x+2>=0).*(x-2<=0).*(x.^2+x)+(x>2).*(2.*x.^2+3.*x);

function p=myfun
%产生正态分布的X，不妨取10万个（这个大小可改，最后概率是一个范围答案）
X=normrnd(40,3,1,100000);
%离散型随机变量的生成
%步骤：产生0~1随机数，0~0.2的属于23，0.2~0.5的属于24，以此类推
r=rand(1,100000);
Y1=(r<0.2).*23+(r>=0.2&r<0.5).*24+(r>=0.5&r<0.8).*25+(r>=0.8).*26;
%找到满足X+Y1<=60的个数
%find函数：k = find(条件) 返回的是下标，即第几个
II=length(find(Y1+X<=60));
p=II/100000;

function Z=myfun
% 生成网格区域：[X,Y] = meshgrid(x,y)
x=-2:0.05:2; y=-2:0.05:2;
[X,Y] = meshgrid(x,y);
%生成矩阵
Z=cos(1+exp(X.*Y)).*sin(X.^2+3.*Y.^2);

function [val,x]=myfun
rand('seed',1321)
%蒙特卡洛：生成很多随机点去做模拟
val=inf;%表达式求min，故初始值设为正无穷
x=[];%初始化最大值点
%进行100000次实验，这个次数可以调整，合理即可
for i=1:100000
    x1=2*rand();
    x2=2*rand();
    if x1^2-x2>=0 && -x1-x2^2+2<=0
        f=x1^2+x2^2+8;
        if f<val %如果函数值比目前的最小值还要小
            val=f;
            x=[x1,x2];
        end
    end
end

function [D1,D2]=myfun
%D1：一阶导数 D2：二阶导数
syms a b x
f=3*a*x^3+a*b*exp(cos(2*x))+2*cos(x);
D1=diff(f,x,1);
D2=diff(f,x,2);

function xp=myfun
%用法： x = fzero(fun,x0)
fun=@(x)x+sin(x)-3;
xp=fzero(fun,2);

function idx=myfun(v)
idx=find(v<=0.6);

function G=myfun
rand('seed',2452)
MB=fix(15*rand(20,35));
G=MB(:,2:8);

function I=myfun
fun=@(x)1./(exp(x)+exp(-x));
I=quad(fun,0,1);

function r=myfun
t=[2.04 3.7 3.96 6.27 8.66 9.57 10.7 11.25 11.28 11.97];
R=[41.898 99.07 111.498 242.084 435.068 520.894 642.71 705.492 711.092 794.855];
% step1:多项式拟合——polyfit
% p = polyfit(x,y,n) n是拟合的阶数，p是多项式系数，降幂排列
p=polyfit(t,R,2);
t_new=3:0.5:40;
% step2:用拟合曲线去计算对应区间函数值——polyval
% y = polyval(p,x) p是polyfit求得的多项式系数
r=polyval(p,t_new);

function r=myfun
%用法:[x,fval] = fsolve(fun,x0)或者 x = fsolve(fun,x0)
%fun是函数fun=0的左半边，x0是初始点；x是解，fval是对应解处函数值
r=fsolve(@fun,[1 1 1]);

%编写一个子函数（因为fun有三个函数组成）
function s=fun(x)
s(1)=3*x(1)-cos(x(2)*x(3))-0.5;
s(2)=x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;
s(3)=exp(-x(1)*x(2))+20*x(3)+(10*pi-3)/3;

function y=myfun
%用法：[t,y] = ode23(odefun,tspan,y0)
% odefun：微分方程，tspan：自变量范围，y0：初始值
% t:自变量（跟输入的tspan相当），y：函数值
x=0:0.01:1;
[t,y]=ode23(@odefun,x,1);
plot(t,y);
function dydt=odefun(x,y)
dydt=2*x+2*x/y;

function H=myfun
rand('seed',1755)
V=fix(49*rand(30,33));
H=V./43;

function D2f=myfun
%采用符号计算的方法
syms f(x)
%f(x)相当于声明f是以x为自变量的函数
f=(3*x+3)^8;
D2f=diff(f,2);
%如果前面没写f(x)，即未声明x为自变量时，应在求导时写diff(f,x,2)
%subs:将对应符号计算中的字母替换成新的字母/向量等等
%使用：snew = subs(目标符号计算式子,旧字母,新字母/矩阵等等)
v=2:0.3:12;
D2f=subs(D2f,x,v);
%题目要求以double形式存储，但是现在仍为符号计算形式，故需转double
D2f=double(D2f);

function a=myfun
x=-2*pi:0.01:2*pi;
a=5.*exp(sin(x))+2.*x.^2-x+1;
plot(x,a);

function [s,m]=myfun
x=0:0.2:50;
f=17.*sin(x).*cos(x)+30.*exp(-6.*x);
%求和：sum 最大值：max 最小值：min
s=sum(f);
m=max(abs(f));

预先准备：先使用matlab画出两个对应的函数图像，观察区域位置。绘图代码：

t=0:0.01:2*pi;
x1=cos(t);
%通过plot预估区域位置
y1=1+sin(t);
plot(x1,y1,'b');
hold on
grid on
%画出原型
x2=-2:0.01:2;
y2=abs(cos(x2));
plot(x2,y2,'r')
axis equal%使得最后的图像正方（不然圆就看着像椭圆了）

可见，包围区域有两个（总共就是一个圆）要求的实际上就是下半区域的面积。

最后的代码实现：

function r=myfun
%step1:在包围区域的对应的矩形区域撒下随机点，个数合理即可
%可知x在-1~1,y在0~2，该矩形区域面积为2*2=4
num=100000;%合理的个数
x=-1+2*rand(1,num);
y= 0+2*rand(1,num);
%step2:找到落入规定区域的随机点
II=length(find(y<=abs(cos(x)) & y>=(1-sqrt(1-x.^2))));
%在圆的曲线上方，在|cosx|曲线的下方
%step3:求得落入目标区域概率
p=II/num;
r=4*p;%概率*面积

function r=myfun
rand('seed',3539);
v=fix(10*rand(1,50));
%find找到符合条件的v的位置
p=find(v>=3&v<=7);
%使用for循环进行求和
r=0;%设置初始值
for j=1:length(p)
r=r+v(p(j));
end

function s=myfun
s=0; r=1.07;
for i=1:19
s=s+r^(2*i);
end

补档：

ode23的使用

function r=ode23fun_2
clear;
clc;
clear all
%ode23解出结果
[t,y]=ode23(@fun,[0:0.1:5],[0;2]);
r=y(:,2);
%存储函数式子
function dfun=fun(t,y)
dfun=[y(1);20*(1-y(2)^2)*y(1)+0.5*y(2)];

 

蒙特卡洛计算二重积分

字符串拼接：

ga的使用

clc;
clear all;


nvars=3;
A=[3 2 6;4 5 2;2 9 7];
b=[50;40;100];
lb=[0;0;0];
ub=[15;9;25];
Aeq=[]; beq=[];
[x,fval] = ga(@func,nvars,A,b,Aeq,beq,lb,ub,@con1,[3])



function myfun1=func(x)
myfun1=2*(x(1)-1)^2+3*(x(2)-1)^2+0.2*x(1)*x(2)+0.2*x(2)*x(3)+(2*x(3)-2)^2;
end

function [C,Ceq]=con1(x)
C(1)=-2*x(1)^2-x(2)^2-x(3)^2+5;
C(2)=x(1)^2+x(2)^2+2*x(3)^2-200;
Ceq=[];
end

 

exprnd使用方法

clear all; clc;
ar(1)=exprnd(5);
 
i=1;
while ar(i)<=30*60
i=i+1;
ar(i)=ar(i-1)+exprnd(5);
end
 
plot(1:i,ar,'r.');
title('到达时刻散点图');
xlabel('学生序号');
ylabel('到达时刻/秒');
grid on