熟话说要站在巨人的肩膀上看世界，该篇建立的数学物理方程的定解问题是张同学提出的，感谢大佬的付出

该篇文章分阶段编写现在是未完成状态

2023-05-20

做完了，我负责的第一题和第三题

2023-04-24

好，总算把常识部分都给借鉴完了，困了，晚安(￣o￣) . z Z

2023-04-23

  1. 按照报告要求，确定任务要求，以及设计原理

一、任务要求

  1. 设计合理，结论正确
  2. 内容规范，过程清晰
  3. 格式规范

(⊙﹏⊙)，等我查一下这个任务要求再补写这一部分！！！

二、设计原理

2.1 定解问题的提出

小提琴弦线设计分析问题可以转化为实际数理问题，即典型弦振动问题。由于主要研究小提琴弦线的E弦，因此问题可以描述为：一根分布均匀柔软的弦弦，两端分别在$x=0$以及$x=L$处固定，设初始时刻弦的形状为一抛物线，抛物线的顶点为$(L/2,h)$，讨论其弦振动的具体问题。<注：初值设定方式不收限制，但是不能全为零，否则无运动，在此以此初值是便于进一步的公式推导，L为弦长>

因此我们可得出定解问题如下：

2.2 求解定解问题

利用分离变量法对两端固定的一维弦振动问题进行求解。

基于两端固定长弦振动问题的结论：$u{(x,t)} = \sum{n=1}^{\infty}{(C_ncos(\frac{n \pi a}{L}t)+D_nsin(\frac{n \pi a}{L}t))sin(\frac{n\pi}{L}x)},n=1,2,3,\cdots$。

由正弦级数的正交性有：$Cn = \frac{2}{L}\int{0}^{L}{\frac{4h}{L^2}(L-x)xsin(\frac{n\pi}{L}x)}dx$，$Dn = \frac{2}{n\pi a}\int{0}^{L}{0}dx = 0$。

求导 $Lx-x^2$ $L-2x$ $-2$ $0$
积分 $sin(\frac{n\pi}{L}x)$ $-\frac{L}{n\pi}cos(\frac{n\pi}{L}x)$ $-(\frac{L}{n\pi})^2 sin(\frac{n\pi}{L}x)$ $(\frac{L}{n\pi})^3 cos(\frac{n\pi}{L}x)$

$C_n = \frac{8h}{L^3}{-(Lx-x^2)\frac{L}{n\pi}cos(\frac{n\pi}{L}x)+(L-2x)[\frac{L^2}{n^2 {\pi} ^2}sin(\frac{n\pi}{L}x)] - 2\frac{L^3}{n^3 {\pi}^3}cos(\frac{n\pi}{L}x)}|^L_0 = \frac{8h}{L^3}[0+0-2\frac{L^3}{n^3 {\pi}^3}cos(n\pi) - (0+0-2\frac{L^3}{n^3 {\pi}^3})]$

综上所述：$C_n = \frac{16h[1-(-1)^n]}{n^3 {\pi}^3}，n=0,1,2,3,\cdots$，$D_n = 0$

至此，定解问题的求解完毕，得出了弦振动的函数表达式：$u{(x,t)} = \frac{32h}{\pi^3}\sum{n=1}^{\infty}{\frac{1}{(2n+1)^3}cos(\frac{(2n+1)\pi at}{L})sin(\frac{(2n+1) \pi x}{L})},n=0,1,2,3,\cdots$

定解问题手推答案过程 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
% 定义常数
L = 0.68; % 弦长
h = 0.6; % 弦顶点高度
v = 380; % 波速
t_max = 0.025; % 模拟时间
dx = 0.001; % 空间步长
dt = dx/v; % 时间步长

% 定义空间网格
x = 0:dx:L;
n_points = length(x);

% 定义初始条件
u = zeros(1, n_points);
for i = 1:n_points
    if x(i) <= L/2
        u(i) = h*(x(i)/((L/2))^2)*(L-x(i));
    end
end

% 定义边界条件
u(1) = 0;
u(end) = 0;

% 迭代求解
u_new = u;
u_old = u;
for time_step = 1:t_max/dt
    for i = 2:n_points-1
        u_new(i) = 2*(1-((v*dt)/dx)^2)*u_old(i) ...
            - u_new(i) + ((v*dt)/dx)^2*(u_old(i+1) ...
            + u_old(i-1));
    end
    u_old = u;
    u = u_new;
end

% 绘制弦的振动图像
plot(x,u);
xlabel('x (m)');
ylabel('u (m)');
title('Small Vibrations of E String of Violin');

2.3 琴弦发音原理推导

经过多个文献的学习得出，小提琴声音是通过把琴弦的振动经过琴码传递到共鸣体，使得共鸣腔体产生共振，带动箱内外空气振动而产生的，而弦的振动则来源于琴弓与琴弦的相互作用。小提琴发声产生的振动大致分为弦的振动、琴码的振动和共鸣箱的振动三种，在此我们仅研究琴弦的振动，正是因此我们所仿真出来的发音效果会带来一定的偏差。琴弦振动分为横振动、纵振动、扭转振动和倍频振动，根据课题要求与实际运算难度，在此我们只研究横振动。
横振动，将弦挑离其平衡位置再放掉，弦就开始做一个扁纱锭型的振动，它的振幅限制在两条明确的曲线之内。弦的横振动频率，可以用泰勒公式来表达$f=\frac{1}{2L}\sqrt{\frac{F}{\rho} }$，其中$\sqrt{\frac{F}{\rho} } = 0$，式中$f$为弦的振动频率，$\rho、F、L$依次为弦的线密度、拉力和长度。且看出横向振动频率与弦长成反比，与张力的平方根成正比，与弦线的密度的平方根成反比。通俗的来讲，拉弦或拨弦可发出声音。琴身内部的空气藉此产生振动而发声，称为共鸣。四条琴弦自小提琴的底部延伸，跨过琴马，至指板尾端以弦轸固定。调音方法是扭动弦轸以调整琴弦的松紧。琴音的高低决定于琴弦的大小、粗细及张力。愈短、愈细、愈紧的琴弦，所产生的音愈高。

2.4 相关参数的确定

有效琴弦长度L：《小提琴制作艺术》【英】克里斯.约翰逊,罗伊.考特诺尔著一书中讲到安装琴颈时琴颈站位与琴身站位比例为2比3，即F孔两个小缺口下口325毫米处。在此我们定弦长L=325mm=0.325m<注：此处我们选用的是4/4M尺寸的小提琴，即成人常用小提琴>

琴弦拉力F:由下表可得H311型号E-Mi弦拉力为18.6磅。

型号 1/16M 1/8M 1/4M 1/2M 3/4M 4/4L 4/4M 4/4H
有效弦长 毫米 216 241 265 290 310 328 328 328
英寸 8 1/2 9 1/2 10 1/2 11 1/2 12 1/4 13 13 13
H311 E-Mi钢弦 10.4 13 13.3 16 16.5 16.8 18.6 20.4
H311W E-Mi钢弦缠铝 16.3 18.1 19.9
H313 A-La多股钢弦缠铝 8 9.3 10.2 11.4 12.2 11.7 12.7 13.6
H312T A-La多股钢弦缠钛 12.2 13.5
H313 D-Re多股钢弦缠钛 7.5 8.3 9.7 10.3 11.3 9.2 11.5 12.2
H314 G-Sol 多股钢弦缠钨银 9.2 10.2 11.4
H314分数琴 G-Sol 多股钢弦缠钨镍 7 7.8 8.9 10 10.3
H315 G-Do 多股钢弦缠钨银 11.1 12.4

因此H311型号E-Mi弦拉力为$18.6磅=8.436kg=82.6728N$。

琴弦线密度:已知E-Mi钢弦的密度为$7.8g/cm^3=7800kg/m^3$,E弦直径为0.276mm则半径为$0.138mm=0.000138m$,由公式线密度=密度x琴弦的横截面积得，E弦的线密度是$4.67\times 10^{-4}kg/m$
至此，我们得到参数：$L=0.325m，F=82.6728N，ρ_线=4.67\times 10^{-4kg}/m$

2.5 小提琴E弦的空弦频谱的音频

要输出小提请E弦的空弦频谱的音频，需要进行一下步骤：

  1. 定义时间和时间步长，以便生成频谱
  2. 生成空白的频谱数据矩阵，大小为采样率乘以时长
  3. 生成E弦的频谱数据，将其添加到音频数据矩阵中
  4. 可以使用soundsc()函数将频谱数据播放，并将其保存为wav文件
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
clc;clear;

% 定义采样率和时间步长
Fs = 44100; %采样率
time = 2; %时长,单位s

%生成空白的频谱数据矩阵，大小为采样率乘以时长
audio = zeros(round(Fs * time), 1);

%生成E弦频谱数据
fE = 659.2551; %E弦的频率 e^2
t = linspace(0, time, round(time * Fs)); %生成横轴(时间轴)
signal = sin(2 * pi * fE * t); %生成正弦波信号，用来产生E弦空弦的音频
signal = awgn(signal, 10, 'measured'); %添加高斯白噪声
signal = reverberator(signal, Fs, 'Time', 0.5, 'WetDryMix', 0.4); % 添加混响效果
audio(:,1) = signal; %将生成的正弦波信号添加到音频数据中

%播放并保存音频
soundsc(audio ,Fs); %播放音频
audiowrite('violin_E4_open_string_spectrum.wav',audio,Fs); %将音频保存为wav文

下面讨论时间步长的选取依据

2.6 时间步长的选取依据

我们要讨论时间步长的选取依据，主要是通过观察瀑布图查看频谱随时间的变化情况，以此确定时间步长。通常我们用以下步骤分析时间步长的选取依据。

  1. 根据采样率和E弦空弦音频的频率分布确定一个初始时间步长。

  2. 计算出音频信号的频谱，并将频谱转换为dB（分贝）单位。

  3. 绘制出频谱随时间的变化曲线（即瀑布图），其中时间步长为横轴，频率为纵轴。颜色表示该时间步长下该信号频率的强度。
  4. 通过观察频谱随时间的变化曲线，根据曲线的变化趋势和频率分布，逐步调整时间步长。重复步骤二和步骤三，直到确定一个好的时间步长。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
clc;clear;
% 读取音频文件
[x, fs] = audioread('E4空弦音频.wav');
fprintf('%d', fs); % 44100

% 计算时间步长
time_step = 0.01; % 初始时间步长
while(time_step <= 0.1) % 设定时间步长的上限
    % 转化时间步长为采样点数
    step = round(time_step * fs); %采样点 441
    % 分帧计算频谱
    [S, f, t] = spectrogram(x, hann(step), floor(step/2), step, fs); % floor(step/2)确保是正数
    % 转化频谱为dB
    S_dB = 20*log10(abs(S));
    % 绘制瀑布图
    waterfall(t, f, S_dB);
    % 设置坐标轴
    view(0, 90);
    xlim([t(1) t(end)]);
    ylim([0 max(f)]);
    xlabel('Time (seconds)');
    ylabel('Frequency (Hz)');
    zlabel('Magnitude (dB)');
    % 等待一段时间
    pause(0.1);
    % 对时间步长进行调整
    time_step = time_step + 0.005; % 设定步长调整值
end

以上，我们在观察时发现极难选取正确的时间步长。但是合适的步长应该满足采样率要求，即时间步长不应该过大也不应该过大或过小，否则会因为过采样或欠采样导致频谱信息失真(上述的代码输出很明显过采样了，都挤在了一起)。

  1. 根据奈奎斯特采样定理，步长应该小于信号最高频率的2倍的倒数，即$dt < \frac {1}{2fc}$。那么在小提琴的情况下，最高频率约为1318Hz。因此可选取步长为$dt = \frac {1}{2659.2551138hz} = 0.0007584317$秒 。
  2. 时间步长应该满足分析需求，即时间步长应该足够小，使频谱图显示出所需的详细信息。

综上所述，选取步长为0.0007584317是相对合适的选择。

下面使用Matlab模拟奈奎斯特定理，对2.6节的输出音频进行分析：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
clc; clear;
% 读取小提琴E弦的音频信号
[x, fs] = audioread('E4空弦音频.wav');

% 生成时间序列
T = length(x)/fs;
t = linspace(0, T, length(x));

% 对信号进行傅里叶变换，求出其频谱
X = fft(x);

% 找到最高频率fmax
f = linspace(0, fs/2, length(x)/2+1);
[~, idx] = max(X(1:length(f)));
fmax = f(idx);

% 计算采样频率fs的最小值
fs_min = 2*fmax;

% 根据fs和T计算出采样点数
N = fs*T+1;

% 选取时间步长k，使得fs/k > fs_min，并且k为正整数
k = ceil(fs/fs_min);

% 对时间序列进行下采样
t_new = t(1:k:end);

% 对下采样后的信号进行重构
x_new = interp1(t, x, t_new, 'spline');

% 对新采样的信号进行傅里叶变换，求出其频谱
X_new = fft(x_new);

% 绘制频谱，观察是否有混叠效应
f_new = linspace(0, fs/k/2, length(x_new)/2+1);
plot(f_new, abs(X_new(1:length(f_new))));

% 显示时间步长
disp(['时间步长为：', num2str(k/fs), '秒']);

可以看到这和理论值相差不多，还不错。

2.7 讨论弦拉力影响琴弦的音频频谱

在基于约束边界讨论弦拉力是如何影响琴弦的音频频谱时，我们可以将该问题拆解成两个基本问题：

  1. 讨论约束边界对小提琴弦拉力的影响。
  2. 讨论弦拉力是如何影响琴弦的音频频谱的。

琴弦的弦拉力可以用如下公式计算：$F = (π d^2 T) / 4L$。其中，F表示琴弦的弦拉力，d表示琴弦的直径，T表示琴弦的张力，L表示琴弦的长度。

根据胡克定律，琴弦的张力可以表示为：$T = (π^2kd^2 f^2) / L^2$。其中，k表示琴弦材料的杨氏模量，f表示琴弦的基波频率。将T代入F的公式中，可以得到：$F = (π^3 k d^4 f^2) / 4L^3$

然而，在实际情况下，约束边界会对琴弦的弦拉力大小产生影响。约束边界会限制琴弦的长度和振动形态，导致琴弦的弦拉力大小与未受约束时不同。考虑琴弦被放置在弦轴和尾块上时的约束边界情况。弦轴和尾块会对琴弦施加约束，导致琴弦的受力分布和振动形态发生变化，从而影响琴弦的弦拉力大小。根据受力平衡条件和振动方程，可以得到满足约束边界条件的琴弦弦拉力公式。然而，由于约束边界情况复杂，公式过于繁琐，难以进行简单推导。基于此，我们优先讨论弦拉力是如何影响琴弦的音频频谱的问题。

2.7.1 弦拉力影响琴弦的音频频谱

弦拉力会影响琴弦的音频频谱，主要表现为以下两方面：

  1. 频率变化：弦拉力越大，琴弦的自然频率就越高，产生的音频频率也就越高。这意味着琴弦音高会变高。相反，如果弦拉力减小，音高就会变低。
  2. 谐波加强：弦拉力对琴弦不同振动模式的谐波产生不同的影响。当弦拉力增加时，高阶谐波会加强，低阶谐波会减弱。这意味着琴弦的音质会更加明亮，富有谐波特性。相反，如果弦拉力减小，低阶谐波会加强，高阶谐波会减弱，产生的音质将更加柔和。

下面使用Matlab进行仿真

  1. 数值模拟准备

首先，我们需要对问题进行数值模拟。我们可以使用波动方程来模拟小提琴弦的振动。具体来说，我们可以根据以下波动方程来模拟弦的振动：$\frac{\partial^2 u}{\partial t^2}=c^2 \frac{\partial^2 u}{\partial x^2}$其中，$u(x,t)$是弦在位置$x$和时刻$t$的位移，$c$是波速。为了简化问题，我们假设弦密度为$\rho$，横截面积为$A$，则波速为$c=\sqrt{T/\rho A}$，其中$T$为弦张力。根据小提琴的实际参数，根据第2.2节相关参数的确定，我们假设小提琴e弦的密度为$4.67\times 10^{-4}$ g/cm，直径约为$0.276$ mm，故横截面积为$A=\pi (0.276/2)^2=0.0491$ mm$^2$，$L=0.325$m。在实验中，我们可以尝试不同的弦拉力$T$，以研究不同弦拉力下的音频频谱变化。

为了进行数值模拟，我们需要对弦进行离散化。我们可以将弦的长度$L$等分为$N$段，每段长度为$h=L/N$。我们将弦在每个节点上的位移$ui=u(x_i,t)$表示为一个向量$\boldsymbol{u}=(u_1,u_2,\cdots,u_N)$。类似地，我们将时间$t$离散化，时间步长为$\Delta t$，时间总长为$T$，则我们可以将弦在每个时间步长$t_n=n\Delta t$的位移表示为一个矩阵$\boldsymbol{U}=(\boldsymbol{u}^n){N\times M}$，其中$\boldsymbol{u}^n$表示弦在时间步长$tn$时各个节点的位移。我们可以使用有限差分法来求解弦的振动，具体来说，我们可以使用中心差分法来求解时间和空间上的导数：$\frac{\partial^2 u}{\partial t^2} \approx \frac{u_i^{n+1}-2u_i^n+u_i^{n-1}}{(\Delta t)^2}$，$\frac{\partial^2 u}{\partial x^2} \approx \frac{u{i+1}^{n}-2ui^n+u{i-1}^{n}}{h^2}$

将上式代入波动方程中，得到：$ui^{n+1}=2u_i^n-u_i^{n-1}+\frac{c^2 \Delta t^2}{h^2}(u{i+1}^n-2ui^n+u{i-1}^n)$其中，$i=1,\cdots,N$，$n=1,\cdots,M$，$M=T/\Delta t$。然后，我们可以根据弦的初始形状、两端的边界条件以及所选择的弦拉力来确定初始条件和边界条件。具体来说，我们可以假设弦的初始形状为一条抛物线，即：$u_i^0=h-\frac{4h}{L^2}(i-\frac{L}{2})^2$其中，$i=1,\cdots,N$。我们可以将弦两端的边界条件设为固定边界条件，即：$u_1^n=u_N^n=0$。最后，我们可以根据所选择的弦拉力$T$计算波速$c$，从而求解出弦在不同时间步长的位移$\boldsymbol{U}$。由于我们希望研究弦振动的频谱变化，因此我们需要对弦的位移进行傅里叶变换，以得到弦振动的频谱分布。具体来说，我们可以将弦的位移$u_i$看作是一个复数$u_i=x_i+jy_i$，其中$x_i$表示实部，$y_i$表示虚部。然后，我们可以对实部和虚部分别进行傅里叶变换，得到弦振动在频域中的频谱分布。最后，我们可以将频谱分布绘制成音频频谱图，以研究不同弦拉力下的音频频谱变化。

  1. MATLAB程序实现

以下为使用MATLAB实现上述数值模拟和分析的程序：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
clc;clear;
% 参数设置
Length = 0.325; % 弦的长度，单位：米
heigh = 0.01; % 抛物线顶点的高度，单位：米
rho = 4.67e-4; % 弦的质量密度，单位：kg/m
Tensions = [10, 20, 30, 40, 50]; % 弦拉力的不同值，单位：牛顿
Nx = 100; % 空间节点数
Nt = 1000; % 时间节点数
Tmax = 0.01; % 最大时间，单位：秒
% 空间和时间步长
dx = Length / (Nx - 1);
dt = Tmax / (Nt - 1);
% 初始化空间和时间网格
x_space = linspace(0, Length, Nx);
t_stride = linspace(0, Tmax, Nt);
% 新建一个窗口
figure;
% 循环求解不同弦拉力下的音频频谱变化并绘制频谱图
for T = 1:length(Tensions)
    c = sqrt(Tensions(T) / rho); % 计算波速
    r = (c * dt / dx)^2;
    % 设置初始条件
    u = zeros(Nx, Nt);
    u(:, 1) = heigh * (1 - (4 * (x_space - Length/2).^2) / Length^2);
    u(:, 2) = u(:, 1);
    % 使用有限差分法求解
    for n = 2:Nt-1
        for i = 2:Nx-1
            u(i, n+1) = 2 * (1 - r) * u(i, n) - u(i, n-1) + r * (u(i+1, n) + u(i-1, n));
        end
    end
    % 计算频谱
    U = abs(fft(u, [], 2));
    freq = (0:size(U, 2)-1) * (1 / t_stride(end));
    % 在子图中绘制频谱图
    subplot(length(Tensions),1,T);
    plot(freq, U);
    xlim([0, 1000]);
    title(['Tension = ', num2str(Tensions(T)), ' N']);
    xlabel('Frequency (Hz)');
    ylabel('Amplitude');
end

得到不同拉力下小提琴e弦音频频谱的频谱图如下：

以上输出的频谱图，经分析得出如下结论：随着弦拉力的增加，频谱的高频部分变得更加明显，而低频部分则变得更加模糊。这是因为弦的振动模态随着拉力的变化而改变，导致不同频率的分量在频谱图中的贡献发生变化，且弦拉力越大，小提琴琴弦的音频频谱最高点呈向上增加的趋势，也就说明，弦拉力越大，琴弦的自然频率就越高，产生的音频频率也就越高。这意味着琴弦音高会变高。相反，如果弦拉力减小，音高就会变低。弦拉力对琴弦不同振动模式的谐波产生不同的影响。当弦拉力增加时，高阶谐波会加强，低阶谐波会减弱。这意味着琴弦的音质会更加明亮，富有谐波特性。相反，如果弦拉力减小，低阶谐波会加强，高阶谐波会减弱，产生的音质将更加柔和。

下面我们根据不同弦拉力生成相对应的小提琴e弦音频文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
clc;clear;

% 参数设置e
Length = 0.325; % 弦长，单位：m
heigh = 0.138; % 弦半径，单位：m
rho = 4.67e-4; % 线密度，单位：kg/m
Fs = 44100; % 采样频率，单位：Hz
duration = 2; % 音频持续时间，单位：s
n_modes = 10; % 考虑的振动模态数量

% 计算系数B_n
x_space = linspace(0, Length, 1000);
B_matux = zeros(1, n_modes);
for n = 1:n_modes
    B_matux(n) = 2/Length * trapz(x_space, heigh * (1 - 4 * x_space.^2 / Length^2) .* sin(n * pi * x_space / Length));
end

% 计算不同拉力下的音频频谱变化
figure;
T_values = linspace(50, 150,5); % 弦拉力范围，单位：N
for i = 1:length(T_values)
    T = T_values(i);
    c = sqrt(T / rho); % 波速，单位：m/s
    
    % 计算振动模态的频率
    freqs = (1:n_modes) * c / (2 * Length);
    
    % 生成音频信号
    t = linspace(0, duration, Fs * duration);
    y = zeros(1, length(t));
    for n = 1:n_modes
        y = y + B_matux(n) * sin(n * pi * t / Length) .* cos(2 * pi * freqs(n) * t);
    end
    
    % 绘制频谱图
    subplot(1,5,i);
    Y = fft(y);
    f = linspace(0, Fs, length(Y));
    plot(f(1:end/2), abs(Y(1:end/2)));
    xlabel('Frequency (Hz)');
    ylabel('Magnitude');
    title(sprintf('E string spectrum at T = %.1f N', T));
    
    % 保存音频文件
    audiowrite(sprintf('E_string_T_%.1fN.wav', T), y, Fs);
end

以上生成音频文件，可以明显听出不同弦拉力下，音频的声音差异，具体表现为：弦拉力越大，琴弦的自然频率就越高，产生的音频频率也就越高，琴弦音高会变高。当弦拉力增加时，琴弦的音质会更加明亮，富有谐波特性。相反，如果弦拉力减小，低阶谐波会加强，高阶谐波会减弱，产生的音质将更加柔和。

2.7.2 约束边界对小提琴弦拉力的影响

以上我们讨论了弦拉力对小提琴音频频谱的影响，我们现在基于约束边界讨论边界条件对弦拉力的影响。

我们可以通过解弦的振动方程来分析不同约束边界对弦拉力的影响。弦的振动方程为:$\frac{\partial^2 u}{\partial t^2} = v^2\frac{\partial^2 u}{\partial x^2}$，其中，$u(x,t)$表示弦的位移，$v$表示波速。

根据初值条件，我们可以得到弦的初始形状为:

弦的两端分别在$x=0$和$x=L$处固定，因此我们可以得到边界条件为:$u(0,t) = 0,\ u(L,t) = 0$或者$\frac{\partial u}{\partial x}(0,t) = 0,\ \frac{\partial u}{\partial x}(L,t) = 0$

根据边界条件，我们可以得到弦的特征方程为:$sin(kL) = 0$，其中，$k$为波数，满足$k=\frac{n\pi}{L}$，$n$为整数。

不同的约束边界对应了不同的特征方程，从而影响了弦的振动模式和拉力。下面我们分别分析不同约束边界对弦拉力的影响。

  1. 自由边界

自由边界表示弦两端不受限制，即边界条件为$u(0,t) = u(L,t) = 0$。此时，特征方程为:$sin(kL) = 0$，解得$k=\frac{n\pi}{L}$，对应的振动模式为:$u_n(x,t) = A_nsin(\frac{n\pi}{L}x)cos(\frac{n\pi}{L}vt)$

弦拉力可以表示为:$T = \mu\frac{\partial u}{\partial x}(0,t) = \mu\frac{\partial u}{\partial x}(L,t)$，其中，$\mu$为弦的线密度。根据弦的振动方程和振动模式，我们可以得到:$T_n = \frac{4\mu h}{n\pi}$

显然，自由边界下弦的拉力与振动模式的波数$n$成反比，即波数越大，拉力越小。这是因为在自由边界下，波的传播是任意方向的，不同波数的波可以互相干涉抵消。

  1. 固定边界

固定边界表示弦两端完全固定，即边界条件为$u(0,t) = u(L,t) = 0$和$\frac{\partial u}{\partial x}(0,t) = \frac{\partial u}{\partial x}(L,t) = 0$。此时，特征方程为:$cos(kL) = 0$

解得$k=\frac{(2n+1)\pi}{2L}$，对应的振动模式为:$u_n(x,t) = A_nsin(\frac{(2n+1)\pi}{2L}x)cos(\frac{(2n+1)\pi}{2L}vt)$，弦拉力可以表示为:$T = \mu\frac{\partial u}{\partial x}(0,t) = \mu\frac{\partial u}{\partial x}(L,t)$

根据弦的振动方程和振动模式，我们可以得到:$T_n = \frac{2\mu h}{L}$

显然，固定边界下弦的拉力与振动模式的波数$n$无关，即所有波数的波都有相同的拉力。这是因为固定边界将弦的振动限制在某些特定的模式下，不同波数的波不能互相干涉抵消。

  1. 半自由边界

半自由边界表示弦的一端固定，另一端自由，即边界条件为$u(0,t) = 0$和$\frac{\partial u}{\partial x}(L,t) = 0$。此时，特征方程为:$\tan(kL) = 0$

解得$k=\frac{n\pi}{L}$，对应的振动模式为:$u_n(x,t) = A_nsin(\frac{n\pi}{L}x)cos(\frac{n\pi}{L}vt)$

弦拉力可以表示为:$T = \mu\frac{\partial u}{\partial x}(0,t) = \mu\frac{\partial u}{\partial x}(L,t)$根据弦的振动方程和振动模式，我们可以得到:$T_n = \begin{cases}\frac{2\mu h}{L} & \text{}n\text{ 为奇数} \ \frac{4\mu h}{n\pi} & \text n\text{ 为偶数}\end{cases}$

显然，半自由边界下弦的拉力与振动模式的波数$n$的奇偶性有关。当$n$为偶数时，拉力与固定边界相同；当$n$为奇数时，拉力与自由边界相同。这是因为半自由边界将弦的振动限制在某些特定的模式下，但仍然允许不同波数的波在自由边界处干涉。当$n$为偶数时，干涉后的波是奇函数，不能在半自由边界处产生位移，从而弦的拉力与固定边界相同；当$n$为奇数时，干涉后的波是偶函数，可以在半自由边界处产生位移，从而弦的拉力与自由边界相同。

综上所述，不同约束边界对弦拉力的影响是显著的。对于分布均匀柔软的弦，弦的拉力与振动模式的波数、约束边界的类型都有关系。在实际应用中，我们需要根据具体情况选择合适的约束边界以及相应的振动模式，以满足设计要求。

下面是四种不同边界条件对线拉力的影响

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
% 定义
L = 0.6; % 弦的长度，单位：米
h = 0.01; % 抛物线顶点的高度，单位：米
rho = 0.01; % 弦的质量密度，单位：kg/m
Tensions = linspace(10, 50, 5); % 弦拉力不同值，单位：牛顿
Nx = 100; % 空间节点数
Nt = 1000; % 时间节点数
Tmax = 0.1; % 最大时间，单位：秒
% 计算间和时间步长
dx = L / (Nx - 1);
dt = Tmax / (Nt - 1);
% 初始化空间和时间网格
x = linspace(0, L, Nx);
t = linspace(0, Tmax, Nt);
% 循环求解不同弦拉力下的音频频谱变化
for T = Tensions
  % 计算波速
  c = sqrt(T / rho);
  r = (c * dt / dx)^2; % 修改了分母
  % 设置初始条件
  u = zeros(Nx, Nt);
  u(:, 1) = h * (1 - (4 * (x - L/2).^2) / L^2);
  u(:, 2) = u(:, 1);
  % 自由边界
  u(1, :) = 0;
  u(end, :) = 0;
  for n = 2:Nt-1
    for i = 2:Nx-1
      u(i, n+1) = 2 * (1 - r) * u(i, n) - u(i, n-1) + r * (u(i+1, n) + u(i-1, n));
    end
  end
  % 计算频率
  U = abs(fft(u, [], 2));
  freq = (0:Nt-1) * (1 / dt / Nt); % 修改了计算频率的方式
  freq1 = freq(2:end);
  U1 = U(:, 2:end);
  [pks, locs] = findpeaks(U1(:, 1));
  freq1 = freq1(locs);
  % 绘制频率-拉力图
  figure;
  plot(T, freq1(1), 'ro');
  hold on;
  plot(T, freq1(2), 'bo');
  plot(T, freq1(3), 'go');
  plot(T, freq1(4), 'yo');
  plot(T, freq1(5), 'mo');
  title('Free Boundary');
  xlabel('Tension (N)');
  ylabel('Frequency (Hz)');
  % 固定边界
  u = zeros(Nx, Nt); % 重新初始化u
  u(:, 1) = h * (1 - (4 * (x - L/2).^2) / L^2);
  u(:, 2) = u(:, 1);
  u(1, :) = 0;
  u(end, :) = 0;
  for n = 2:Nt-1
    u(2, n) = r * (u(3, n) - u(2, n-1)) + 2 * (1 - r) * u(2, n) - u(2, n-2); % 修改了公式
    u(end-1, n) = r * (u(end-2, n) - u(end-1, n-1)) + 2 * (1 - r) * u(end-1, n) - u(end-1, n-2); % 修改了公式
    for i = 3:Nx-2
      u(i, n+1) = 2 * (1 - r) * u(i, n) - u(i, n-1) + r * (u(i+1, n) + u(i-1, n));
    end
  end
  % 计算频率
  U = abs(fft(u, [], 2));
  freq = (0:Nt-1) * (1 / dt / Nt); % 修改了计算频率的方式
  freq2 = freq(2:end);
  U2 = U(:, 2:end);
  [pks, locs] = findpeaks(U2(:, 1));
  freq2 = freq2(locs);
end
    
  %绘制频率-拉力图
  figure;
  plot(T, freq1(1), 'ro');
  hold on;
  plot(T, freq1(2), 'bo');
  plot(T, freq1(3), 'go');
  plot(T, freq1(4), 'yo');
  plot(T, freq1(5), 'mo');
  title('Free Boundary');
  xlabel('Tension (N)');
  ylabel('Frequency (Hz)');
    
  %自由-固定边界
  u(1, :) = 0;
  u(end, :) = 0;
  for n = 2:Nt-1
      u(2, n) = r * (u(3, n) - u(2, n-1)) + 2 * (1 - r) * u(2, n) - u(2, n-1);
      u(end-1, n) = r * (u(end-2, n) - u(end-1, n-1)) + 2 * (1 - r) * u(end-1, n) - u(end-1, n-1);
      for i = 3:Nx-2
          u(i, n+1) = 2 * (1 - r) * u(i, n) - u(i, n-1) + r * (u(i+1, n) + u(i-1, n));
      end
  end
    
  % 计算频率
  U = abs(fft(u, [], 2));
  freq = (0:size(U, 2)-1) * (1 / t(end));
  freq3 = freq(2:end);
  U3 = U(:, 2:end);
  [pks, locs] = findpeaks(U3(:, 1));
  freq3 = freq3(locs);
  % 绘制频率-拉力图
  figure;
  plot(T, freq1(1), 'ro');
  hold on;
  plot(T, freq1(2), 'bo');
  plot(T, freq1(3), 'go');
  plot(T, freq1(4), 'yo');
  plot(T, freq1(5), 'mo');
  title('Free Boundary');
  xlabel('Tension (N)');
  ylabel('Frequency (Hz)');

三、附录

3.1 抛物线的求解

已知抛物线的定点和焦点，可以使用抛物线的标准方程或顶点法求出抛物线的方程。

  1. 标准方程

抛物线的标准方程为：

其中，a是抛物线的开口方向和大小，b是抛物线的对称轴位置，c是抛物线的顶点的纵坐标。

根据题意，已知定点和焦点，可以求出抛物线的开口方向和大小，以及顶点的横坐标。进而可以求出b和c。

以定点为抛物线的顶点，设焦点在y轴正半轴上，设焦距为2p，定点的坐标为(h,k)。

根据抛物线的定义，抛物线上任意一点到焦点的距离等于该点到定点的距离，即：

平方化化简可得：

其中p为焦距。根据定义，抛物线的开口方向和大小由参数a决定，a=1/4p。

因此抛物线的方程为：

其中，h，k，p为已知的定点和焦点的参数。

  1. 顶点法

以定点为抛物线的顶点，将焦点所在的直线定义为抛物线的对称轴。设焦距为2p，焦点的坐标为（h，k+p），则可得到焦点所在的直线方程为x=h。

根据抛物线的定义，顶点到对称轴的距离等于焦距的一半，因此抛物线的方程为：

其中，h，k，p为已知的定点和焦点的参数。

3.2 奈奎斯特采样定理

采样过程所应遵循的规律，又称取样定理、抽样定理。采样定理说明采样频率与信号频谱之间的关系，是连续信号离散化的基本依据。
在进行模拟/数字信号的转换过程中，当采样频率fs.max大于信号中最高频率fmax的2倍时(fs.max>2fmax)，采样之后的数字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的2.56～4倍；采样定理又称奈奎斯特定理。
如果对信号的其它约束是已知的，则当不满足采样率标准时，完美重建仍然是可能的。 在某些情况下（当不满足采样率标准时），利用附加的约束允许近似重建。 这些重建的保真度可以使用Bochner定理来验证和量化。具体见百度百科奈奎斯特抽样定理