利用蒙特卡罗法求解热传导方程

作者: 深雨露 | 来源:发表于2020-01-17 21:09 被阅读0次

利用蒙特卡罗法求解热传导方程
使用scipy库的root和fsolve函数求解方程
热传导方程
Matlab使用ode45求解高阶微分方程
scipy一维热传导方程的求解
2020-01-16
非线性方程求解、幂法
第23课微分方程和exp(At)
数学工具
学子情缘

蒙特卡罗方法简介

蒙特卡罗（Monte Carlo, abbr. MC）方法是利用独立重复的统计实验来对物理及数学问题求解的方法。一个简单的关于MC方法的应用即求解图形的面积。

例：求解 $(x-\frac{1}{2})^2+(y-\frac{1}{2})^2=\frac{1}{2}$ 的面积？

%投点次数
M = 1000;
A = rand(M,2);
r = sqrt((A(:,1) - 0.5).^2 + (A(:,2) - 0.5).^2);
total = sum(r < 0.5);
pr = total / M;
disp(pr)

%绘图
xt = @(t) 0.5*cos(t) +0.5;
yt = @(t) 0.5*sin(t) + 0.5;
fplot(xt,yt)
hold on
plot(A(:,1),A(:,2),'.');
grid on
axis equal 
axis([0,1,0,1])

计算结果：0.7630

计算结果

MC求解面积的数学原理便是弱大数定理的方法

设 $X_1,X_2,\cdots$ 相互独立，服从同一分布的随机变量序列，且具有数学期望 $E(X_k) = \mu \quad (k=1,2,\cdots)$ ，作前 $n$ 个变量的算术平均 $\frac{1}{n}\sum_{k=1}^n X_k$ ，则对任意的 $\varepsilon> 0$ ，有
$\lim_{n \rightarrow \infty} \Pr \left\lbrace \left| \frac{1}{n} \sum_{k=1}^n X_k -\mu \right|<\varepsilon \right\rbrace =1.$

可以简单的理解为若统计的数据足够大，则事物的频率 $\frac{1}{n} \sum_{k=1}^n X_k$ 可无限接近于其期望 $\mu$ 。

利用MC方法求解一维热传导方程的边值问题

下面以2018年数学建模国赛A题为例，利用MC方法进行求解其第一问。

题目：热防护服的设计

人们在温度较高的环境下工作时，通常都要穿上专门的服装来避免灼伤。一般的避免灼伤专用服装是由三层的织物材料构成，分布标记为 I、II、III 层，而其中的 I 层与外界环境接触，III 层与皮肤之间还存在着空隙，所以我们将此空隙标记为 IV 层。为了降低研发专用服装成本、缩短研发设计专用服装的周期，首先把体内温度控制在 37ºC 的假人放置在实验室的高温环境中，然后建立数学模型来测量假人皮肤外侧的温度，并且解决下列问题：

（1）附件1给出了专用服装材料的一些参数值，附件 2 给出了在环境温度为 75 ℃、II 层的厚度为 6 毫米、IV 层的厚度为 5 毫米、工作时间为 90 分钟的条件下假人皮肤外侧的温度的变化情况。请建立适当的数学模型来计算温度的分布，最后把计算好的温度生成为 Excel 文件。

附件1. 专用服装材料的参数值

分层	密度( $kg/m^3$ )	比热 ( $J/(kg\cdot$ ºC))	热传导率 (W/(m·℃))	厚度 ( $mm$ )
I层	300	1377	0.082	0.6
II层	862	2100	0.37	0.6-25
III层	74.2	1726	0.045	3.6
IV层	1.18	1005	0.028	0.6-6.4

附件2. 假人皮肤外侧的测量温度

时间 (s)	温度 (℃)	时间 (s)	温度 (℃)
0	37.00	5396	48.08
1	37.00	5397	48.08
2	37.00	5398	48.08
3	37.00	5399	48.08
4	37.00	5400	48.08

解析题目

对于这种方法能量传递类问题，一般应考虑利用热传导方程来进行求解
$\frac{\partial u}{ \partial t} = k \cdot \nabla^2u,$
其中 $k$ 是热传导系数。对于如题的一维问题，可简化上述方程为
$\frac{\partial u}{\partial t} = k_{\iota} \cdot \frac{\partial^2 u}{\partial x^2}, \quad x \in [0,L],\iota \in \{1,2,3,4\}$
其中 $k_i$ 代表不同层次的热传导系数（如第 I 层的热传导系数 $k_1 = 0.082$ ）。

完善初始条件和边界条件

初始条件
$u(x,0) = \varphi(x),$
边界条件
$u(0,t)=\mu_1(x), \quad u(L,t)=\mu_2(x).$
至于在两层的界面之间的联系则依据热流量密度 $k_\iota \frac{\partial u}{\partial x}$ 连续的原理进行定义
$(k_\iota \frac{\partial u}{\partial x})^- = (k_{\iota+1}\frac{\partial u}{\partial x})^+$
下面使用MC方法求解热传导方程 $^{[1]}$ 。取差分距离点：对距离 $x$ 进行差分
$x_j = j \Delta x, \quad j = 1,2,\cdots;$
对时间 $t$ 进行差分
$t_n = n \Delta t, \quad n=1,2,\cdots$
利用六点半隐式差分方法对原方程差分化
$\frac{{u_j^{n + 1} - u_j^n}}{{\Delta t}} = \frac{k}{{2(\Delta {x)^2}}}[(u_{j + 1}^{n + 1} - 2u_j^{n + 1} + u_{j - 1}^{n + 1}) + (u_{j + 1}^n - 2u_j^n + u_{j - 1}^n)]$

其中拉布拉斯算子采用 $(n+1)$ 层和 $n$ 层的中央差分和的平均

$(n+1)$ 层
$\frac{\partial u}{\partial x} \approx \frac{(u_{j + 1}^{n + 1} - 2u_j^{n + 1} + u_{j - 1}^{n + 1})}{(\Delta x)^2}$
$n$ 层
$\frac{\partial u}{\partial x} \approx \frac{(u_{j + 1}^{n } - 2u_j^{n } + u_{j - 1}^{n })}{(\Delta x)^2}$
合并两层

$\frac{\partial u}{\partial x} \approx \frac{1}{2} \cdot \left[ \frac{(u_{j + 1}^{n + 1} - 2u_j^{n + 1} + u_{j - 1}^{n + 1})}{(\Delta x)^2} + \frac{(u_{j + 1}^{n } - 2u_j^{n } + u_{j - 1}^{n })}{(\Delta x)^2} \right]$

令 $\lambda = k \frac{\Delta t}{(\Delta x)^2}$ 得

$- \frac{\lambda }{2}u_{j + 1}^{n + 1} + (1 + \lambda )u_j^{n + 1} - \frac{\lambda }{2}u_{j - 1}^{n + 1} = \frac{\lambda }{2}u_{j + 1}^n + (1 - \lambda )u_j^n + \frac{\lambda }{2}u_{j - 1}^n$

考虑 $\lambda \leqslant 1$
$\begin{align} u_j^{n + 1} &= \frac{1}{{1 + \lambda }}\left[ {\frac{\lambda }{2}u_{j + 1}^{n + 1} + u_{j - 1}^{n + 1} + \frac{\lambda }{2}u_{j + 1}^n + u_{j - 1}^n + 1{\rm{ - }}\lambda u_j^n} \right] \\ &={\alpha _1}u_{j + 1}^{n + 1} + {\alpha _2}u_{j - 1}^{n + 1} + {\alpha _3}u_{j + 1}^n + {\alpha _4}u_{j - 1}^n + {\alpha _5}u_j^n \end{align}$
其中 $\sum_{i=5}^5 \alpha_i =1$ .

依据上式可以得知， $u_j^{n+1}$ 的求得与 $u_{j + 1}^{n + 1} ,u_{j - 1}^{n + 1} ,u_{j + 1}^n ,u_{j - 1}^n ,u_j^n$ 这五项及其相应对 $u_j^{n+1}$ 的贡献权重/比例（即 $\alpha_i$ ）密切相关。

那么我们可设想人随机走动的模式：位于 $(j,n+1)$ 处的人分别以概率 $\alpha_1,\alpha_2,\alpha_3,\alpha_4,\alpha_5$ 向临近点移动。确定待移位置后，人到达待移位置，继续以类似的形式移动，直到人到达边界处（ $n=1$ 或者边界点 $Q_{k1}$ ）.将其赋值为 $f(Q_{k1})$ .

原理图

对一个 $u(x_j,t_{n+1})$ 的子样 $f(Q_{k1})$ 重复 $M$ 次掷点，从而得到 $f(Q_{k1}),f(Q_{k2}),\cdots,f(Q_{kM})$ 。那么可认为 $u(x_j,t_{n+1})$ 的无偏估计为
$\widehat{u_M}=\frac{1}{M} \sum_{i=1}^{M} f(Q_{ki})$
依据上面的思路，即可求解出整个的热量传输过程的温度变化

计算机实现步骤

为加快求解速度，采用MATLAB中MATLAB Coder应用生成mex文件后（使用JIT编译器）进行生成。上述应用只能对函数进行编译，因而下文均在函数——MC_PDE中进行书写。下文以第IV层计算为例予以说明。

运行实现

定义函数的参数和返回值

function u = MC_PDE(temp_90,M)
% temp_90 ：90分钟外界边界条件
% M = 10：投点次数

初始条件设置

%热传导方程系数
k = 0.028;
%内层边界条件（皮肤）
edge_first = repmat(37,1,5401);
%长度维起始行
edge_location = 1;
edge_location = int32(edge_location);
% 长度维终止点
l_end = 50;
l_end = int32(l_end);

%% x方向
% x方向划分间隔
h = 0.1;
% x方向：值域，x ∈ ( 0, x_max)
x_max = 15.2;
%x_initial = 0 : h :1;
% x方向迭代次数
x_times = x_max / h - 1;
x_times = int32(x_times);

%% t方向
% t方向划分间隔
tau =1;
% 时间值域：t ∈ ( 1, t_max)
t_max = 5401;
% 时间迭代次数
t_tims = t_max / tau;
t_tims = int32(t_tims);

temp_total = 0;

%系数求解
r = k * ( tau / (h^2) );
coef = 1 / (1 + r);

%%u初始化
u = zeros(x_times+1,t_tims) * nan;
% 初值条件
u(:,1) = 37;

% 边值条件
u(edge_location,:) = edge_first;
u(x_times+1, :) = temp_90(:,2)';

随机走动部分

for j_value = (edge_location + 1) : l_end
    disp(j_value)
    for n = 2:5401
        % 投掷点
        temp = zeros(1,M);
        for point_num = 1 : M
            t = int32(n);
            l = int32(j_value);
            while(1)
                xi = rand();
                if xi > 0 && xi < coef * (1/2 * r)
                    l = l + 1;
                    
                    %检验函数
                    if l == 152
                        basicvalue = u(152,t);
                        temp(point_num) = basicvalue;
                        signal = 1;
                    elseif   l == 1
                        basicvalue = u(1,t);
                        temp(point_num) = basicvalue;
                        signal = 1; 
                    elseif t == 1
                        basicvalue = u(l,1);
                        temp(point_num) = basicvalue;
                        signal = 1;
                    else
                        temp(point_num) = nan;
                        signal = 0;
                    end

                    if signal == 1
                        break
                    end

                elseif xi >= coef * (1/2 * r) && xi < coef * ( r)
                    l = l - 1;

                    %检验函数
                    if l == 152 || l == 1 || t == 1 
                        basicvalue = u(l,t);

                        temp(point_num) = basicvalue;
                        signal = 1;
                    else
                        temp(point_num) = nan;
                        signal = 0;
                    end

                    if signal == 1
                        break
                    end
                elseif xi >= coef * ( r) && xi < coef
                    t = t - 1;

                    %检验函数
                    if l == 152 || l == 1 || t == 1 
                        basicvalue = u(l,t);

                        temp(point_num) = basicvalue;
                        signal = 1;
                    else
                        temp(point_num) = nan;
                        signal = 0;
                    end

                    if signal == 1
                        break
                    end
                elseif xi >= coef && xi < coef * (1 + (1/2)*r)
                    l = l + 1;
                    t = t - 1;

                    %检验函数
                    if l == 152 || l == 1 || t == 1 
                        basicvalue = u(l,t);

                        temp(point_num) = basicvalue;
                        signal = 1;
                    else
                        temp(point_num) = nan;
                        signal = 0;
                    end

                    if signal == 1
                        break
                    end
                elseif xi >= coef * (1 + (1/2)*r) && xi < 1
                    t = t - 1;

                    %检验函数
                    if l == 152 || l == 1 || t == 1 
                        basicvalue = u(l,t);

                        temp(point_num) = basicvalue;
                        signal = 1;
                    else
                        temp(point_num) = nan;
                        signal = 0;
                    end

                    if signal == 1
                        break
                    end
                end
            end    
        end
        
        %期望计算部分
        temp_total = sum(temp);
        u(j_value, n) =  temp_total / M;
    end
end

运行结果分析

下面结论依据掷点数为1000的条件得到。运行时间为10803.82s（约3h）。

三维图

不可否认，由于MC方法求解的特殊性，温度随着时间和厚度的变化并不均匀。这一点可以从上图不断凸起的低矮“小山”这一现象中可知。

二维热图

二维热图可清楚可见三道明显的分界线，进热源点随时间温度变高，符合构建模型预期。

总结

从计算复杂度看，复杂度略显偏大。求解速度呈现出求解区域靠近边界附近运算速度快，中间较为慢的特点，这个特点与这个算法的机制密切相关（因为中间的点难以足够快的速度移动到边界）。
从求解方法看，MC法相比隐式求解求解而言，对某个多维点而言，不必求解出各个层次的值，只需对这个多维点取足够多的子样即可求解。

$[1]$ 刘春光. 偏微分方程边值问题的蒙特卡罗解法[D]. 吉林大学, 2004.

利用蒙特卡罗法求解热传导方程
蒙特卡罗方法简介蒙特卡罗（Monte Carlo, abbr. MC）方法是利用独立重复的统计实验来对物理及数学...
使用scipy库的root和fsolve函数求解方程
求解非线性方程使用scipy库的 root, fsolve 函数求解非线性方程。求解传热方程（热辐射+热传导）...
热传导方程
热传导方程方程及其定解问题的导出齐次热传导方程：非齐次热传导方程：当物体体积很大，考虑时间很短和较小范围...
Matlab使用ode45求解高阶微分方程
原问题为利用ode45求解下面微分方程：初值为可以转为二阶微分方程利用Matlab求解，下面为源代码 fun...
scipy一维热传导方程的求解
1.源码实现 2.运行程序 3.执行结果
2020-01-16
机器学习概念统计学样本 —— 数据集待求解方程（待求解方程中的未知数为参数） —— 预训练模型求解参数后的N元方程...
非线性方程求解、幂法
实验六 test1 利用roots函数求解非线性方程：解： fzero函数样例运行测试画出以下非线性函数，利用f...
第23课微分方程和exp(At)
怎么求解一阶方程？怎么求解一阶导数？怎么求解常系数线性方程？将它们转换成线性代数问题的思路是常系数线性方程解...
数学工具
desmos - 函数图形绘制方程求解
学子情缘
记得刚上初二，数学进入方程求解的教学内容，如一元n次方程的求解，二元一次方程组的求解等。一元二次方程有公式法，有因...