一根均匀弦两端固定在x=0， x=l处。假设初始时刻速度为零，而在初始时刻弦的形状是一条顶点为(l/2，h)的抛物线，求弦振动的位移。用matlab仿真将结果以图形显示出来，然后分析一下l，h以及T以及线密度变化时对弦振动的影响

要解决这个问题，我们需要应用波动方程，考虑到弦的振动是简谐的。由于弦两端固定，我们可以确定振动模式为驻波，其形式通常由傅立叶级数描述。在这种特定情况下，因为弦的初始形状是一条抛物线，这提示我们可能需要使用数值方法，例如有限差分法来逼近解析解。

首先，让我们用Matlab编写一个函数来计算弦的振动位移。我们将使用欧拉-贝特定理来近似弦的振动。考虑到弦的固定端点，只有奇次谐波（正弦波）才有非零振幅。这里，`h`决定了抛物线的峰值，而`l`影响了波长。

function u = solveHarmonicOscillator(l, h, dt, dx, T)
    % 初始化
    L = l;
    N = floor(L/dx);
    x = 0:dx:L;
    
    % 弹簧常数 k (假设均匀弦，k = density * tension)
    density = 1; % 设定线密度，单位通常是kg/m
    tension = 1; % 设定张力，单位通常是N
    
    k = density * tension;
    A = h / sqrt(k); % 抛物线顶点到原点的距离
    
    % 定义初始条件 (抛物线)
    u0 = A * (1 - x.^2 ./ L^2);
    
    % 数值解
    t = 0:dt:T;
    [t, U] = ode45(@(t,u) harmonic_oscillator(t, u, x, L, k), t, u0);
    u = reshape(U(:, end), length(x));
    
    % 函数体...
end

function du_dt = harmonic_oscillator(t, u, x, L, k)
    % 根据波动方程 du/dt = -sqrt(k/L^3) * sin(x/L) * u
    du_dt = -sqrt(k/L^3) .* sin(x/L) .* u;
end

% 示例运行并可视化
[lValues] = 0.1:0.1:1;
hValues = 0.5:0.1:1.5;
figure;
hold on;
for i = 1:length(lValues)
    for j = 1:length(hValues)
        l = lValues(i);
        h = hValues(j);
        T = 2*pi/l; % 最小周期
        [x, u] = solveHarmonicOscillator(l, h, 0.01, 0.01, T);
        plot(x, u, '-o', 'MarkerSize', 3, sprintf('l=%f, h=%f', l, h));
    end
    hold off;
    xlabel('Distance along the string');
    ylabel('Displacement');
    title('Vibrating String with Different l and h Values');
    drawnow;
end

当`l`增大（即弦变得更长），周期`T`增加，因此频率减小；`h`（初始高度）会影响振动的最大振幅，但不会直接影响周期。线密度`density`一般只影响波的速度（通过弹簧常数`k`），并不会显著改变振动模式本身。张力`tension`会同时影响速度和最大振幅。

相关问题：
1. 这种数值方法是否适用于所有类型的初始条件？
2. 当张力`tension`发生变化时，如何修改波动方程和计算？
3. 对于实际的弦乐器，这种模型有哪些局限性？