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)
当时出门没有电脑也没有键盘,用安卓移动端调试了好久才得数诶: