问题及解答

在化工生产中常常需要知道丙烷在各种温度 $T$ 和压强 $P$ 下的导热系数 $K$. 下面是实验得到的一组数据:

试求 $T=99$°C 和 $P=10.3\times 10^3$ kPa 下的 $K$.

[分析]

这里可以将 $K$ 看成 $T$ 和 $P$ 的函数, 记为 $K=K(T,P)$.

在有限的数据下, 只能通过(合理的)插值去求出在特定点处的值.

已知 $K(T[j],P[j])=K_j$, 求 $K(x,y)$. 利用散乱数据插值法. 令

\[
r_j=\sqrt{(x-x_j)^2+(y-y_j)^2}
\]

$K(T,P)$ 定义为:

\[
 K(T,P)=\begin{cases}
K_j, & \text{当}\ r_j=0 \text{时},\\
\sum_{k=1}^{N}W_k(x,y)z_k & \text{当}\ r_j\neq 0 \text{时},\\
\end{cases}
\]

其中

\[W_k(x,y)=\dfrac{\dfrac{1}{r_k^2}}{\sum_{k=1}^{N}(1/r_k^2)}.\]

为此, 建立下面的函数fun1, 保存为 fun1.m 文件.

function z= fun1(T,P,K,x,y)
% 已知 K(T[j],P[j])=K_j, 求 K(x,y).
% 利用散乱数据插值法
% r_j=sqrt((x-x_j)^2+(y-y_j)^2)
% K(T,P) 定义为:
% K(T,P)=K_j, 当 r_j=0 时;
% 当 r_j ~=0 时,
% K(T,P)=\sum_{k=1}^{N}W_k(x,y)z_k
% where W_k(x,y)=\frac{1}{r_k^2}/\sum_{k=1}^{N}(1/r_k^2).
%-------------------------------------
n1=size(T,2);
n2=size(P,2);
n3=size(K,2);

if n1==n2 && n2==n3
    bool_flag=0;
    R=zeros(1,n1);
    W=R;
    sum=0;
    for j=1:n1
        r2=(x-T(j))^2+(y-P(j))^2;
        if r2==0
            bool_flag=1;
            z=K(j);
            break;
        end
        R(j)=r2;
        sum=sum+1/r2;
    end
    
    if bool_flag==0
        for i=1:n1
            W(i)=1/(R(i)*sum);
        end
        z=W*K';        
    end
else
    z=inf;
    
end

end

然后编写主程序 Ex7_6_1.m

% T （单位：摄氏度）
% P （单位：1000kPa）
% K=K(T,P), （K 导热系数）
%------------------------------------
% 在化工生产中常常需要知道丙烷在各种温度(T)和压强(P)下的热导系数K.
% 下面是实验得到的一组数据. 试求 T=99℃, P=10.3x10^3 kPa 下的热导系数K.
%
% 考虑使用散乱数据插值（参见书本P.158）
%--------------------------------------
clear
clc

T=[68 68 87 87 106 106 140 140];
P=[9.7981 13.324 9.0078 13.355 9.7918 14.277 9.6563 12.463];
K=[0.0848 0.0897 0.0762 0.0807 0.0696 0.0753 0.0611 0.0651];

x=99;
y=10.3;

z=fun1(T,P,K,x,y);
z

'---test---'
z=fun1(T,P,K,68,13.324);
z

figure
plot(T,P,'*r',x,y,'+')

figure
plot3(T,P,K,'-.or')