Matlab 语法小问题 【问答随笔·5】

发布于:2023-01-22 ⋅ 阅读:(439) ⋅ 点赞:(0)

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)

当时出门没有电脑也没有键盘,用安卓移动端调试了好久才得数诶:
在这里插入图片描述

本文含有隐藏内容,请 开通VIP 后查看