固体介质的一维热传导方程

Question

这是一段含有一维热传导方程的Matlab代码

function [loss] = compute_loss(rho, c, k, measured_times, t, dt, dx, nt, nx)
    alpha = k / (rho * c);
    r = alpha * dt / dx^2;
    T = ones(nx,nt) * 100; % 初始温度条件

    if r > 0.5
        error('选择的dt和dx值会导致解不稳定。尝试减小dt或减小dx');
    end

    % 求解热传导方程
    for i = 1:nt-1
        L_boundary_derivative = 0.01 * t(i)^2;
        for j = 2:nx-1
            T(j, i+1) = r * T(j + 1, i) + (1 - 2*r) * T(j, i) + r * T(j - 1, i);
        end
        T(end, i+1) = T(end-1, i+1) + L_boundary_derivative * dx;
    end

    % 初始化模拟的超声波传播时间
    start_times = 1:10; % 1秒、2秒、...、10秒
    echo_time = zeros(1, length(start_times));
    
    % 根据指定的起始时间点计算超声波传播时间
    for idx = 1:length(start_times)
        start_timestep = round(start_times(idx) / dt); % 转换秒到时间步
        current_T = T(:, start_timestep);
        V_initial = -0.4251 * current_T + 3260;
        one_way_time = sum(dx ./ V_initial); 
        if start_timestep + one_way_time/dt <= nt
            current_T = T(:, round(start_timestep + one_way_time/dt));
        else
            current_T = T(:, end); % 如果超出了时间范围，则使用最后的温度分布
        end
        V_return = -0.4251 * current_T + 3260;
        return_time = sum(dx ./ V_return);
        echo_time(idx) = one_way_time + return_time;
    end

    % 计算误差
    loss = sum((echo_time - measured_times).^2);
end

可以帮我看一下这个函数的问题吗？我将传入的参数rho从2000遍历到20000，但最后出来的结果都是一样的

%主程序
clc;
clear; 
% 设置参数
dt = 0.1;
dx = 0.1;
x = 0:dx:50;
t = 0:dt:10;
nt = length(t);
nx = length(x);
k = 50; % 导热系数
c = 400; % 比热
measured_times =[0.003138678215, 0.003138974298, 0.003139260027, 0.003139544182, 0.00313982763, ...
                 0.003140110528, 0.003140392974, 0.03140675219, 0.003140957191, 0.003141239013];

for i=2000:20000
    loss = compute_loss(i, c, k, measured_times, t, dt, dx, nt, nx);
    disp(loss)
end

Answer 1

引用 皆我百晓生 小程序回复内容作答：
这段代码中的问题是，你在计算热传导方程时，初始温度条件设置为 T = ones(nx,nt) * 100，这使得你的温度分布在整个传播时间内都是相同的，因此最终的结果也会是相同的。

你需要根据实际情况设置初始温度条件，这样才能得到不同的温度分布。一种常见的设置方法是在初始温度条件中考虑温度梯度，即根据物体的初始温度和边界条件来设置初始温度分布。

另外，你的 for 循环中只输出了误差的值，没有保存每个 rho 对应的误差值。你可以使用一个数据结构（如数组或矩阵）来保存误差值，以便最后进行分析或可视化。