Matlab语法常见小问题解答（问答随笔第5期）

文章目录

• 1. Matlab坐标点基准线 2. isequal()返回总是0 3. 计算定积分报错

1. Matlab坐标点基准线

问：请问这两条线什么画啊? 也不知道该怎么称呼它，所以也没查到。


   答：坐标点的竖线是直接用的火柴杆图stem()，只不过坐标点花式了些；水平横线是用的line()，对坐标点取水平基准线：

clc,clear,close all;
x = 0:5;
y = 2.*x;
stem(x,y)
hold on
x0 = zeros(1,length(y));
line([0 1],[2,2])  %表示(0,2)与(1,2)之间有水平连线，以此类推画出其他水平横线

   line()函数具体用法参考：https://ww2.mathworks.cn/help/matlab/ref/line.html?s_tid=doc_ta

stem()函数具体用法参考：https://ww2.mathworks.cn/help/matlab/ref/stem.html?searchHighlight=stem&s_tid=srchtitle_stem_1

2. isequal()返回总是0

问：matlab相同的数字判断为不等



   复数的模的平方应该等于复数乘以其共轭，但是为什么matlab显示逻辑不正确呢？

   我也试了很多数据，都不对，后来又用了isequal()，也不对，连十六位精度都是一样的也判断不等？这是为什么？

答：isequal()只要是有一点不相等也是不行的，

z = 2+2i
z =
   2.0000 + 2.0000i
>> z1 = abs(z)^2
z1 =
    8.0000
>> z2 = z*conj(z)
z2 =
     8

这两个8虽然都是double类型，但其实精度不同，z2的8才是真正的整数8；

   z1:

而且如果让他俩相减的话就可以看出具体差多少：

>> z1-z2
ans =
   1.7764e-15

在这里将z1向下取整一下变为整数的double，两者就相等了：

>> isequal(floor(z1),z2)
ans =
  logical
   1

3. 计算定积分报错

最后一个Bug同样来自于未被采纳的问题，但仍具参考价值（今天 问答  被采纳破100了，

由此发篇随笔纪念一下）——利用Matlab求解复杂积分时遇到了这样一个报错：

Warning: Non-finite result. The integration was unsuccessful. Singularity likely

意思是积分域内部的奇点警告，对此，Mathworks给出了官方解决文档：Portal：积分域内部的奇点

但我的被积函数有点复杂…一时半会不知道该如何修改积分区域：

其中

积分化简：

我的Bug来源程序：

clc,clear;
format long
fo = 4;  %输入更改Fo
fo = sqrt(fo);
fun = '(x2*x3)/(x1^2)*exp(-(x2^2+x3^2)/(4*x1^2))*exp((x2*x3*cos(x4))/(2*x1^2))';
up = {num2str(fo),'0*x1+4','0*x2+1','0*x3+1'};
low = {'0','0*x1','0*x2','0*x3'};
f = nIntergrate(fun,low,up);
f = f/pi
function f = nIntergrate(fun,low,up)
% f,返回值，n重积分积分结果
% fun,是被积函数字符串形式；
% low 存储从外到内各重积分的积分下限函数；
% up 存储从外到内各重积分的积分上限函数（都是字符串形式）
% 为了保证函数正常运行，low和up内的函数遵循如下原则：设积分重数为m，最内层
% 到最外层积分变量依次以xm,...x2,x1来表示（规定只能以这样的顺序而且只能字母x代表变量）
% 最内层的变量积分范围不管是否和x(m-1)相关都要求出现x(m-1)，不是的可以在积分范围
% 后加上‘+0*x(m-1)’等形式的字符串，依次类推，次内层要求x(m-2)出现，一直到最
% 外层，最外层才是常数
% 
% 参考书籍：MATLAB高效编程技巧与应用：25个案例分析 91页
N = length(low);
fun = vectorize(fun);
if verLessThan('MATLAB','7.8')
    expr = GenerateExpr_quadl(N);
else
    if mod(N,2) == 0
        expr = GenerateExpr_quad2d(N);
    else
        expr = ['quadl(@(x1) arrayfun(@(x1)',GenerateExpr_quad2d(N-1),',x1),',low{1},',',up{1},')'];
    end
end
%============================================================
% 主要利用quadl函数递归构造求解表达式，适用于R2009a以前的版本
%============================================================
    function expr = GenerateExpr_quadl(n)
        if n == 1
            expr = ['quadl(@x',num2str(N),')',fun,',',low{N},',',up{N},')'];
        else
            expr = ['quadl(@(x',num2str(N-n+1),') arrayfun(@(x',num2str(N-n+1),')',GenerateExpr_quadl(n-1),',x',num2str(N-n+1),'),',low{N-n+1},',',up{N-n+1},')'];
        end
    end
%============================================================
% 主要利用quad2d函数递归构造求解表达式，适用于R2009a以及以后的版本
%============================================================
    function expr = GenerateExpr_quad2d(n)
        if n == 2
            expr = ['quad2d(@(x',num2str(N-1),',x',num2str(N),')',fun,',',low{N-1},',',up{N-1},',@(x',num2str(N-1),')',low{N},',@(x',num2str(N-1),')',up{N},')'];
        else
            expr = ['quad2d(@(x',num2str(N-n+1),',x',num2str(N-n+2),')','arrayfun(@(x',num2str(N-n+1),',x',num2str(N-n+2),')',GenerateExpr_quad2d(n-2),',x',num2str(N-n+1),',x',num2str(N-n+2),'),',low{N-n+1},',',up{N-n+1},',@(x',num2str(N-n+1),')',low{N-n+2},',@(x',num2str(N-n+1),')',up{N-n+2},')'];
        end
    end
f = eval(expr);
end

答：检查了很久，一直计算不出结果，后来还是拆成四个简单积分，用符号积分一步步算出来了

clc,clear;
 
syms x1 x2 x3 x4
f1 = exp((x2*x3/(2*x1^2))*cos(x4));
v1 = int(f1,x4,0,pi) 
f2 = exp(-(x2^2+x3^2)/(4*x1^2))*v1*x3;
v2 = int(f2,x3,0,1)
f3 = x2*v2;
v3 = int(f3,x2,0,1)
f4 = v3/(x1^2)/pi;
v4 = int(f4,x1,0,2)  %最后一个2代表Fo
f = vpa(v4,3)

当时出门没有电脑也没有键盘，用 安卓 移动端调试了好久才得数诶：

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删