地壳中的热传导是经典的边界条件随时间变化的热传导问题，因为地表温度随季节变化 假设地表某一点的日平均温度变化如下 To(t) = A + Bsin(2rt),这里t=365天,A=10°C,B =12°C.在地表以下20米全年的温度可以近似为11°C(高于地表的平均温度10C---由于地心的热量影响温度随着深度的增加而增加)地壳的热扩散率因地而异为简化问题设热扩散率是常数a2=D=012mday-1. 写程序计算地壳温度变化,深度到20米时间上限是10年除了表面和离地表20米处的温度其它地方的初始温度均设为10°(选择栅格点和时间步长h,首先运行出9年的结果在第10年画出4个温度轮廓某个轮廓代表3个月的时间,这四个轮廓在一幅图中，显示温度随深度和时间的变化

这是一个较为复杂的问题，需要使用更加高级的数值方法对其进行求解。一种常用的方法是有限元法，下面给出一个基于有限元法的 Python 代码示例。

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags, csr_matrix

# 定义模型参数
D = 0.12  # 热扩散率
a = np.sqrt(D / np.pi)  # 热扩散长度
T_surface = 22.0  # 地表温度
T_average = 11.0  # 地下平均温度
A = 10.0  # 日平均温度中的常数项
B = 12.0  # 日平均温度中的振幅
r = np.pi / 182.5  # 日平均温度中的频率

# 定义模型空间和时间范围
L = 20.0  # 模型深度
N = 200  # 空间网格数
M = 3650  # 时间步数
dx = L / N
dt = 86400  # 时间步长，单位为秒

# 定义初始温度分布
T = np.ones(N+1) * T_average

# 定义时间变化的表面温度函数
def To(t):
    return A + B * np.sin(2 * r * t)

# 构建有限元矩阵
A = diags([-1, 2, -1], [-1, 0, 1], shape=(N-1, N-1))
A *= D * dt / dx**2
A = csr_matrix(A)

# 进行时间迭代
for n in range(M):
    # 更新表面温度
    T[0] = To(n * dt)

    # 计算温度变化
    Tn = np.copy(T[1:N])
    T[1:N] = A.dot(Tn) + T[1:N]

# 绘制温度轮廓图
depths = np.linspace(0, L, N+1)
times = np.linspace(0, M*dt, 5)
for t in times:
    T_surface = To(t)
    T_plot = np.hstack((T_surface, T))
    plt.plot(T_plot, depths)
plt.xlabel('Temperature (°C)')
plt.ylabel('Depth (m)')
plt.ylim(L, 0)
plt.show()

这段代码使用了有限元矩阵来求解地壳中的热传导问题，并绘制了温度轮廓图。注意，这里的时间步长为一天，总时间步数为 3650，总时间为 10 年。由于计算量较大，这段代码可能需要一些时间才能运行完毕。