活动介绍
file-type

MATLAB中三角函数误差及符号工具箱解决方案

ZIP文件

下载需积分: 50 | 4KB | 更新于2025-04-22 | 119 浏览量 | 2 下载量 举报 2 收藏
download 立即下载
在讨论MATLAB中三角函数在pi或pi/2等特殊点的误差问题之前,我们需要先了解几个关键概念:圆周率π,浮点数表示,以及MATLAB中的符号计算。 首先,圆周率π是数学中的一个常数,表示圆的周长与直径的比例,其值约为3.14159。π是无理数,意味着它是一个不能准确表示为两个整数比例的实数,并且它的小数表示是无限且不循环的。在计算机程序中,我们无法存储一个精确的π值,只能使用π的一个近似值。这个近似值通常是一个浮点数,它是在计算机科学中用于表示实数的一种方法。浮点数由三部分组成:符号位、指数部分和尾数部分。由于有限的存储空间,浮点数只能提供数值的一个近似表示,这在数值计算中可能引入舍入误差。 在MATLAB中,pi是一个内置常量,用于代表圆周率π的近似值。由于它是一个浮点数,当我们在计算过程中使用这个值进行三角函数计算时,就可能遇到由于数值精度限制而导致的误差。例如,当计算sin(pi)或cos(pi/2)时,MATLAB中的浮点近似值不能完全准确地表示这些三角函数的理论值,因此计算结果可能非常接近于零,但又不完全为零。 这里的关键问题在于浮点数的表示误差。对于像π/2这样精确值为整数倍π/2的情况,计算得到的结果为一个非常小的数(例如cos(pi/2)的计算结果大约为6.12323399573677e-17),这反映了浮点数的精度限制。由于浮点数不能精确表示π,当其值被用在任何基于π值的计算中时,都可能引入误差。这种误差通常在数学上可以接受,但在一些对精度要求极高的应用场景中就显得不可忽视。 为了解决这种由浮点数近似值引入的误差问题,MATLAB提供了一种符号计算工具箱。使用符号工具箱中的函数可以进行符号计算,也就是说,可以在没有数值精度损失的情况下计算精确的数学表达式。例如,使用sin(sym(pi))代替sin(pi)可以在MATLAB中得到精确的0,因为sym(pi)代表的是一个精确的π值,而不是浮点数近似值。符号计算工具箱扩展了MATLAB的功能,允许用户进行没有数值近似的精确计算,这在很多工程和科研领域都是非常有用的。 MATLAB中的符号计算工具箱提供了一组丰富的函数,可以在数学表达式中使用符号变量,进行代数运算、微分、积分、极限计算等。这些计算结果是精确的,能够提供数学意义上的精确解,而不受计算机浮点数精度的限制。 总结一下,在MATLAB中进行三角函数计算时,如果涉及到了圆周率π的精确值,可能会遇到由浮点数近似值引起的误差。在大部分情况下,这种误差是可以接受的,但是在需要极高精度的场合,可以使用MATLAB的符号计算工具箱来避免这种误差。通过使用符号变量,用户可以得到准确的数学结果,而不必担心数值精度带来的影响。这在工程计算、科学分析以及教育演示等领域具有重要的应用价值。

相关推荐

filetype

一、设计缘由本文采用的是汉宁窗来设计滤波器,窗函数的好处如下:窗函数可以控制滤波器的频率响应:窗函数通常用于在频域上调整FIR滤波器的频率响应,从而满足设计要求。不同的窗函数具有不同的频域特性,因此选择合适的窗函数可以产生所需的滤波器响应。窗函数可以控制滤波器的截止频率:FIR滤波器通常通过选择截止频率来滤除信号中的不需要的频率成分。窗函数可以控制滤波器的截止频率,使滤波器具有所需的通带和阻带特性。窗函数可以改善滤波器的性能:使用窗函数可以改善滤波器的幅度响应和相位响应的性能,从而减少滤波器的失真和波纹。二、设计步骤1)确定选用的窗函数选择窗函数-要根据过度带宽和阻带最小衰减来选取合适的窗函数,本文以汉宁窗为例。六种窗函数的基本参数窗函数类型 旁瓣峰值/dB 过度带宽 阻带最小衰减/dB矩形窗     -13 -21三角窗 -25 -25汉宁窗 -31 -44汉明窗 -41 -53布莱克曼窗 -57 -74凯赛窗 -57 -80本题目中的通带下截止频率为400Hz,阻带下截止频率为350Hz,因此,过度带宽B为(经过归一化处理)。通过计算可得到窗函数需要的阶数N为124阶。2)构造希望逼近的频率响应函数,即理想滤波器的上下截止频率近似位于最终设计的FIRBF的过渡带宽的中心频率点,幅度函数衰减一半(约-6dB)。截止频率一般取。3)计算。直接使用MATLAB中的fir1函数,hn=fir1(N-1,wc,hanning(N))注意:fir1函数默认选用hamming窗。三、MATLAB实现代码如下clcclear%% ****************************settings*****************************fs=2000; %采样频率 fpl=400;fpu=500; %通带频率fsl=350;fsu=550; %阻带频率wpl=2*pi*fpl/fs;wpu=2*pi*fpu/fs; %通带角频率(归一化)wsl=2*pi*fsl/fs;wsu=2*pi*fsu/fs; %阻带角频率Bt=wpl-wsl; %过渡带带宽N0=ceil(6.2*pi/Bt); %计算滤波器的阶数wc=[(wpl+wsl)/2/pi,(wpu+wsu)/2/pi]; %计算理想带通滤波器截止频率N=N0+mod(N0+1,2); %保证长度为奇数 f1=450;f2=600;T=1; %时宽1sL=round(T*fs); %采样点个数 %% *****************************filter settings**********************n=0:N-1; hn=fir1(N-1,wc,hanning(N)); %调用fir1计算带通滤波器的h(n)[H,w]=freqz(hn,1,L); %计算频率响应函数 figure(1)magH=20*log10(abs(H)/max(abs(H))); %幅频响应subplot(311);stem(n,hn,'.');xlabel('n');ylabel('h(n)');title('汉宁窗FIR-BPF的单位脉冲响应');subplot(312);plot(w/pi*fs/2,magH);title('汉宁窗FIR-BPF的幅频响应')xlabel('频率/Hz');ylabel('20lg|Hg(e^j^\omega)|/max(|Hg(e^j^\omega)|)');grid onsubplot(313)plot(w/pi*fs/2,unwrap(angle(H))); %相频响应title('汉宁窗FIR—BPF的相频响应')xlabel('频率/Hz');ylabel('相位/rad')grid on figure(2);z=roots(hn);zplane(z,[])title('汉宁窗FIR-BPF的零极点分布图');legend('零点','极点');grid on; %% ******************************filter signals (time)**********************dt=1/fs;%t=0:dt:(L-1)*dt;t=linspace(0,T,L);x1=sin(2*pi*f1*t);x2=sin(2*pi*f2*t);x=x1+x2;% y=filter(hn,1,x);y=filtfilt(hn,1,x);% 解决滤波器时延% y=conv(hn,x); figure(3);subplot(411);plot(t,x1);title('450Hz正弦信号x1(t)');xlabel('t/s');ylabel('x1(t)'); subplot(412);plot(t,x2);xlabel('t/s');ylabel('x2(t)');title('600Hz正弦信号x2(t)'); subplot(413);plot(t,x);title('混合信号x(t)');xlabel('t/s');ylabel('x(t)'); subplot(414);% t1=0:dt:(L+124-1)*dt;plot(t,y);xlabel('t/s');ylabel('y(t)');title('滤波信号y(t)'); %% ***************************filter signals (frequnecy)***********************M=length(x);X=abs(fftshift(fft(x./(L))));x_fs=linspace(-fs/2,fs/2-1,L);% x_fs=(0:M-1)*fs/M-fs/2; %频率向量Y=abs(fftshift(fft(y./(L)))); figure(4);subplot(211);plot(x_fs,X);xlabel('Frequency'); ylabel('Amplitude'); title('混合后信号频谱')subplot(212);plot(x_fs,Y);xlabel('Frequency'); ylabel('Amplitude'); title('滤波后信号频谱')四、结果分析下图为MATLAB所绘制的频率响应图,我们可以从相频响应图中看到相位在通带范围内是线性。这是FIR滤波器的一大特点。 下图为FIR滤波器的零极点分布图,可看到大部分零点是在单位圆上的,少部分是在单位圆外的;FIR滤波器在原点z=0处是有N-1个极点,因此FIR滤波器是绝对稳定的。这也是FIR滤波器的一大特点。 下图是对450Hz正弦信号和600Hz正弦信号混合后的信号进行滤波,滤除了600Hz的信号,留下了450Hz的正弦信号。 下图是滤波前和滤波后的频谱图,可以直观的看出450Hz的信号被保留了下来。 五、关于吉布斯效应的分析吉布斯效应:将具有不连续点的周期函数(如矩形脉冲)进行傅立叶级数展开后,选取有限项进行合成。当选取的项数越多,在所合成的波形中出现的峰起越靠近原信号的不连续点。当选取的项数很大时,该峰起值趋于一个常数,大约等于总跳变值的9%。这种现象称为吉布斯效应。(1)在理想特性不连续点ωc附近形成过渡带。过滤带的宽度近似等于 WR(θ)主瓣宽度,Δω=4π/N 。(2)通带内增加了波动,最大的峰值在ωc- 2π/N 处。阻带内产生了余振,最大的负峰在ωc+2π/N处。通带与阻带中波动的情况与窗函数的幅度谱有关。 WR(θ)波动愈快(加大时),通带与阻带内波动愈快, WR(θ)旁瓣的大小直接影响波动的大小。这些影响是对hd(n)加矩形窗引起的,称之为吉布斯效应。 增加矩形窗口的宽度N不能减少吉布斯效应的影响。N的改变只能改变ω坐标的比例和  的绝对大小,不能改变主瓣和旁瓣幅度相对值。加大N并不是减少吉布斯效应的有效方法。寻找合适的窗函数形状,使其谱函数的主瓣包含更多的能量,相应旁瓣幅度就变小了;旁瓣的减少可使通带与阻带波动减少,从而加大阻带的衰减。但这样总是以加宽过渡带为代价的。

filetype

%% 参数设置 clear all; close all; epsilon = 0.01; % 摄动参数 N = 100; % 总网格数(偶数) max_iter = 10000; % 最大迭代次数 tol = 1e-6; % 收敛容差 omega = 1.7; % SOR松弛因子 %% 生成Shishkin网格 sigma = min(0.5, 2*epsilon*log(N)); % 过渡点位置 M = N/2; % 每侧网格数 % 生成网格坐标 x_left = linspace(0, sigma, M+1); % 左侧细网格 x_right = linspace(sigma, 1, M+1); % 右侧粗网格 x = unique([x_left, x_right]); % 合并网格 n = length(x)-1; % 实际区间数 %% 计算步长 h = diff(x); %% 构造系数矩阵 A = sparse(n-1, n-1); b = zeros(n-1, 1); for i = 1:n-1 xi = x(i+1); % 二阶导数项 (非均匀网格中心差分) h_prev = h(i); h_next = h(min(i+1, length(h))); coeff = -epsilon * 2 / (h_prev + h_next); % 填充三对角部分 if i > 1 A(i, i-1) = coeff/h_prev; end A(i, i) = A(i, i) - coeff*(1/h_prev + 1/h_next); if i < n-1 A(i, i+1) = coeff/h_next; end % 对流项 (二阶迎风) b_coeff = -(1 + xi*(1 - xi)); if i < n-2 % 足够点进行二阶差分 h1 = h(i); h2 = h(i+1); denom = h1*(h1 + h2); A(i, i) = A(i, i) + b_coeff*(-(2*h1 + h2))/denom; A(i, i+1) = A(i, i+1) + b_coeff*((h1 + h2)^2)/(h2*denom); A(i, i+2) = A(i, i+2) + b_coeff*(-h1^2)/(h2*denom); else % 边界附近用一阶迎风 A(i, i) = A(i, i) + b_coeff/h(i); if i < n-1 A(i, i+1) = A(i, i+1) - b_coeff/h(i); end end % 反应项 A(i, i) = A(i, i) + 2; % 右端项计算 D = 1 - exp(-1/epsilon); term1 = (2*(1 - exp(-xi/epsilon)) - (xi*(1-xi)/epsilon)*exp(-xi/epsilon))/D; term2 = (1 + xi*(1-xi))*(pi/2)*cos(pi*xi/2); term3 = -(epsilon*(pi^2/4) + 2)*sin(pi*xi/2); b(i) = term1 + term2 + term3; end %% SOR迭代求解 u = zeros(n-1, 1); % 初始猜测 for iter = 1:max_iter u_old = u; for i = 1:n-1 sigma = A(i,:)*u - A(i,i)*u(i); u(i) = (1-omega)*u(i) + omega*(b(i) - sigma)/A(i,i); end if norm(u - u_old, inf) < tol break; end end %% 结果处理 u = [0; u; 0]; % 添加边界值 %% 精确解计算 u_exact = (1 - exp(-x'/epsilon))/(1 - exp(-1/epsilon)) - sin(pi*x'/2); %% 绘图比较 figure; plot(x, u_exact, 'r-', 'LineWidth', 2); hold on; plot(x, u, 'b--', 'LineWidth', 2); xlabel('x'); ylabel('u(x)'); legend('Exact Solution', 'Numerical Solution', 'Location', 'Best'); title(['\epsilon = ', num2str(epsilon), ', N = ', num2str(N)]); grid on; %% 误差分析 error = norm(u - u_exact, inf); disp(['最大误差: ', num2str(error)]); disp(['收敛迭代次数: ', num2str(iter)]);这个代码有问题,帮我改正一下

filetype

一切可以从自由波动方程开始说起: (其中v为波速) 考虑最理想的情况:忽略重力,忽略厚度的均质二维方形平面,只有其中心处被固定 以该方形平面为xy平面,固定中心为原点,垂直平面方向为z轴,建立空间直角坐标系 从而自由波动方程化为: 设z具有分离变量的形式解: 代入波动方程化简后有: 对于等式左边,引入分离常数,设 ( ) 整理后即 ,类似于单摆的动力学方程 从而时谐项对应平面上各点的简谐振动: 其中A和B为与初始条件相关的常数 对于原波动方程等式右侧: 继续分离变量,引入分离常数: ( ) 类似地,有 ( ) 显然三个分离常数之间满足一个约束条件: 同理,X和Y也有与T形式相似的通解: 现在开始讨论边界条件: 首先注意到,Chladni图案讨论的始终是二维平面上的驻波 一方面,由于整个平面只有中心点处是被约束的,因而要形成波节,则一定要经过中心点 另一方面,参考甩动绳子时形成的一维驻波,绳子自由的那一端通常应该是波腹 因而这里假定二维方形平面自由振动的边缘同样也全部都是波腹 题外话1:我并不太确定这样的假设是否合理,实际上网上更多文献中假定平面边缘为波节,但大多数的Chladni图案的制造操作中,并没有任何固定平面边缘的迹象,因而我觉得波节假设不太合理 这样一来限定了平面上可能的振动模式,量化了振动模 对于简单的正弦/余弦波,振动模的一维示意图如下(设正方形平板的边长为L): 正弦模 余弦模 题外话2:其实一开始建系时取方形平面中心为原点,感觉上好像更对称,但实际上使表达式更复杂了,大部分文献中采取将方形平板的一个顶点置于原点,整体置于第一象限的做法 量化之后的X和Y的表达式为(其中n和m均取非负整数): 其实上式中的待定乘数因子alpha与beta,以及一开始求出的时谐函数T,在这个问题中我们是并不关心的,因为对于驻波,比如Chladni图案,需要关注的是波动函数中始终为0,不随时间演化而改变的波节,乘子和时间演化并不影响波节的位置,因此现在最重要的是研究: 这一类隐函数对应的图像 但上式对应的一类图像和实验中观察的Chladni图案相似度并不高,它们只是简单地分别沿x方向和y方向将平板用直线分为几个矩形部分,每个矩形部分单独振动而已: m和n同时取1~5时的波节图案 这是因为对于同一种振动模式(m,n),x轴和y轴实际上是等价的,也就是说(n,m)和(m,n)对应了两个“简并”的振动模,类比于(其实离得有点远)量子力学中全同粒子构成的对称波函数和反对称波函数,这两种“简并”模(m,n)和(n,m)也可以构成“对称”模(m,n)+(n,m)和“反对称”模(m,n)-(n,m),对应的波动函数即: 其中 而 若以m和n分别为横纵指标在网格中分别画出对应各种模式(m,n)的波节图 显然无论是对称模还是反对称模,其在对称指标(m,n)和(n,m)下对应的图像都一样 另外注意到,当取n=m时,反对称模不存在 因此构造这样一个模的表格: 上三角部分(不包括主对角线)为反对称模 下三角部分(包括主对角线)为对称模 特别注意,由于原点固定为波节的边界条件限制,对称模中m和n均为偶数时,对应的振动模式不满足该边界条件,应除去(下图中均未除去) m和n分别取从0到20的整数,作图如下: 放大查看 网上经常能找到的正方形钢板上的Chladni图案记录[1]: 和之前大图中的局部对比: 再回过头来看题图: 在第一张大图中可以找到如下对应: 大图中3行5列,对应题图440Hz 大图中11行13列,对应题图880Hz 大图中11行19列,类似题图2400Hz 大图中15行19列,类似题图4000Hz 部分Mathematica代码如下(参考该网站[2]修改,取波速v和边长L均为单位1): \[Psi][x_, y_, p_, q_] := Cos[p \[Pi] (x - 1/2)] Cos[q \[Pi] (y - 1/2)]; With[{cols = RGBColor /@ {"#F66095", "#2BCDC1", "#393E46"}}, Table[ContourPlot[\[Psi][x, y, p, q] + Sign[p - q] \[Psi][x, y, q, p] == 0, {x, -0.5, 0.5}, {y, -0.5, 0.5}, Axes -> False, Frame -> False, ContourStyle -> Directive[CapForm["Round"], Thickness[.01], Blend[cols[[;; -2]]]], MaxRecursion -> 2, PlotRangePadding -> -0.01, ImageSize -> 54, Background -> cols[[-1]]], {p, 0, 10, 1}, {q, 0, 10, 1}]]这些是部分内容和代码,请你提取其中的信息并写出MATLAB代码

filetype

``` %% 悬臂梁参数设置(来自论文4.1节) lb = 0.447; % 梁长(m) wb = 0.0305; % 梁宽(m) tb = 0.7e-3; % 梁厚(m) Eb = 70e9; % 弹性模量(Pa) rho_b = 2770; % 密度(kg/m³) Cb = 0.01; % 阻尼系数 % 压电作动器参数(A/B在根部,C在中段) lp = 0.0526; % 压电片长度(m) wp = 0.01; % 压电片宽(m) tp = 3.5e-3; % 压电片厚(m) d31 = 1.7e-10; % 压电常数(m/V) Ep = 7.1e10; % 压电片弹性模量(Pa) Cp = 3e-9; % 电容(F) Ka = (wb*tb*d31*Ep)/(2*Cp); % 压电耦合系数 % 压电片位置(A/B在根部,C在中段) xA1 = 0.012; % A左边界位置(m) xA2 = xA1 + lp; % A右边界位置(m) xB1 = xA2 + 0.01; % B左边界位置(m) xB2 = xB1 + lp; % B右边界位置(m) xC1 = lb/2 - lp/2;% C左边界位置(m) xC2 = lb/2 + lp/2;% C右边界位置(m) %% 动力学模型(状态空间方程) % 模态参数(假设前两阶模态) w1 = 15.8; % 第一阶模态频率(rad/s) w2 = 98.4; % 第二阶模态频率(rad/s) xi1 = 0.02; % 阻尼比 xi2 = 0.02; % 状态矩阵A, B(来自论文方程5-6) A = [0 0 1 0; 0 0 0 1; -w1^2 0 -2*xi1*w1 0; 0 -w2^2 0 -2*xi2*w2]; % 输入矩阵B(双驱动:B和C作动器) phi11 = calculate_phi(xB1, xB2, 1); % 需定义calculate_phi函数 phi12 = calculate_phi(xC1, xC2, 1); phi21 = calculate_phi(xB1, xB2, 2); phi22 = calculate_phi(xC1, xC2, 2); B = Ka * [0 0; 0 0; phi11 phi12; phi21 phi22]; % 输出矩阵C(末端位移) C = [1 1 0 0]; % 假设模态叠加为末端位移 sys = ss(A, B, C, 0); % 状态空间模型 %% DRNN神经网络初始化 inputSize = 2; % 输入:控制误差[error1, error2] outputSize = 3; % 输出:PID参数[Kp, Ki, Kd] hiddenSize = 5; % 隐含层节点数 % DRNN权重初始化 Wx = randn(hiddenSize, inputSize)*0.01; Wy = randn(hiddenSize, hiddenSize)*0.01; Wh = randn(outputSize, hiddenSize)*0.01; bh = zeros(hiddenSize, 1); bo = zeros(outputSize, 1); % 训练参数 learningRate = 0.1; state = zeros(hiddenSize,1); % 隐含层状态 %% PID控制器初始化 Kp = 1.0; Ki = 0.1; Kd = 0.01; % 初始PID参数 integral = 0; previous_error = 0; %% 仿真参数 dt = 0.001; % 时间步长 T = 2; % 总时长 t = 0:dt:T; % 输入信号:作动器A的正弦+随机扰动 freq = 10; % 正弦频率(Hz) Amp = 50; % 幅值(V) UA = Amp*sin(2*pi*freq*t) + 5*randn(size(t)); % 初始化变量 X = zeros(4, length(t)); % 状态变量 Y = zeros(1, length(t)); % 输出位移 U_control = zeros(2, length(t)); % 控制输入(B和C) %% 主仿真循环 for k = 1:length(t)-1 % 当前误差(目标位移为0) error = -Y(k); % DRNN前向传播调整PID参数 input = [error; integral]; state = tanh(Wx*input + Wy*state + bh); pid_params = Wh*state + bo; Kp = pid_params(1); Ki = pid_params(2); Kd = pid_params(3); % PID控制计算 integral = integral + error*dt; derivative = (error - previous_error)/dt; U_pid = Kp*error + Ki*integral + Kd*derivative; % 解耦控制(动态调整输入分配) U_control(:,k) = [0.6*U_pid; 0.4*U_pid]; % 简化解耦 % 系统输入 = 扰动A + 控制B/C U_total = [UA(k); U_control(1,k); U_control(2,k)]; % 状态更新(使用ode45求解微分方程) [~, x] = ode45(@(t,x) A*x + B*U_total(2:3), [0 dt], X(:,k)); X(:,k+1) = x(end,:)'; Y(k+1) = C*X(:,k+1); % DRNN在线训练(误差反向传播) target = [1.2*error; 0.8*integral; 0.5*derivative]; % 目标PID参数 delta_o = (pid_params - target).*(1 - tanh(pid_params).^2); delta_h = (Wh'*delta_o) .* (1 - state.^2); % 权重更新 Wh = Wh - learningRate * delta_o * state'; Wx = Wx - learningRate * delta_h * input'; Wy = Wy - learningRate * delta_h * state'; end %% 可视化结果 figure; subplot(2,1,1); plot(t, Y); title('悬臂梁末端位移'); xlabel('时间(s)'); ylabel('位移(m)'); subplot(2,1,2); [Pxx,f] = pwelch(Y, [], [], [], 1/dt); semilogy(f, Pxx); title('功率谱密度'); xlabel('频率(Hz)'); ylabel('PSD'); % 辅助函数:计算振型函数导数差 function delta_phi = calculate_phi(x1, x2, mode) % 假设振型函数(实际需根据论文公式实现) if mode == 1 phi_prime = @(x) cos(pi*x/2); else phi_prime = @(x) -sin(3*pi*x/2); end delta_phi = phi_prime(x2) - phi_prime(x1); end```末端位移始终为0 ,并且他的解耦对吗

weixin_38620893
  • 粉丝: 4
上传资源 快速赚钱