题目写的看不懂。我也瞎答。
写了个MATLAB的小程序，用特征线法求解 偏微分方程组。
[XX,TT]=meshgrid(0:0.4:4,0:0.1:1);
N=size(XX,2);
T=size(XX,1);
u1=zeros(T,N);
u2=zeros(T,N);
dx=0.4;
dt=0.1;
ds=0.1;
u10=1:N;u20=(N:-1:1); % initial value
u10t=1:T;u20t=1:T; % value at a given node point
u1(1,:)=u10;
u2(1,:)=u20;
u1(:,1)=u10t;
u2(:,1)=u20t;
for tt=2:T % 特征线法求解线性偏微分方程组
for ii=2:N
u1(tt,ii)=(-400*u1(tt-1,ii-1)+400*u2(tt-1,ii-1))*ds;
u2(tt,ii)=(0.35*u1(tt-1,ii)-0.35*u2(tt-1,ii))*dt;
end
end
quiver(XX,TT,u1,u2)

1  引言
微分方程是描述一个系统的状态随时间和空间演化的最基本的数学工具之一，其在物理、经济、工程、社会等各方面都有及其重要的应用。然而，只有很少的微分方程可以解析求解，尤其对于偏微分方程，能解析求解的种类更是寥寥可数。更多的微分方程可以采用数值法进行求解，只要精度足够高，就可以满足科学和工程上的需求。
数值求解微分方程的基本思路是先把时间和空间离散化，然后将微分化为差分，建立递推关系，然后利用计算机强大的重复计算能力，快速得到任意格点处的值。Python的Numpy、Scipy工具包可以很好地实现此功能，Matplotlib工具包则可以将求解结果画为非常直观的图形。
接下来，我们先以常微分方程的数值求解为例，引入差分的思想，再将其推广到偏微分方程中。
2  常微分方程的差分求解
一般地，一阶常微分方程可以写为
首先，将连续的变量
 和 
 离散化，连续的函数 
 和 
 化为离散的序列 
 和 
 ，则上述微分方程可以化为差分方程【1】
从而我们得到递推关系
有了递推关系和初始条件以后，就可以利用Python的强大计算功能，得到任意的
 的值了，下面我们通过一个具体的例子来说明。
2.1  RC回路放电问题
对于一个
 回路，我们有
其中，
 分别为电流，电阻，电量和电容，利用 
 ，并定义 
，我们得到一个含初始条件的一阶常微分方程
这个方程当然可以解析求解，得到
 。另一方面按照差分法，可以得到递推关系
下面我们用Python进行数值求解，并和解析结果进行比较。
import numpy as np
import matplotlib.pyplot as plt
rc = 2.0 #设置常数
dt = 0.5 #设置步长
n = 1000 #设置分割段数
t = 0.0 #设置初始时间
q = 1.0 #设置初始电量
#先定义三个空列表
qt=[] #用来盛放差分得到的q值
qt0=[] #用来盛放解析得到的q值
time = [] #用来盛放时间值
for i in range(n):
t = t + dt
q1 = q - q*dt/rc #qn+1的近似值
q = q - 0.5*(q1*dt/rc + q*dt/rc) #差分递推关系
q0 = np.exp(-t/rc) #解析关系
qt.append(q) #差分得到的q值列表
qt0.append(q0) #解析得到的q值列表
time.append(t) #时间列表
plt.plot(time,qt,'o',label='Euler-Modify') #差分得到的电量随时间的变化
plt.plot(time,qt0,'r-',label='Analytical') #解析得到的电量随时间的变化
plt.xlabel('time')
plt.ylabel('charge')
plt.xlim(0,20)
plt.ylim(-0.2,1.0)
plt.legend(loc='upper right')
plt.show()图1：用改进的欧拉法差分得到的结果和解析结果的比较
在图1中给出了差分法得到的结果与解析法得到结果的比较，发现两者符合得很好，这说明对于这个问题，改进的欧拉法已经可以给出足够精确的结果。需要注意的是，这个微分方程本身比较简单，可以解析求解，而对于复杂得多的微分方程，没法解析求解，但是上述数值求解差分方法仍然是适用的。
3  偏微分方程的差分求解
有了差分代替微分的思想，接下来我们将其推广到偏微分方程的求解中。以一般二阶抛物型偏微分方程为例，一般的可以写为
仍然是将时间和空间离散化，将微分化为差分，即
其中
 和 
 分别为空间步长和时间步长， 
 和 
 分别标记空间指标和时间指标，则我们得到差分方程
由此得到递推关系
下面我们考察一个具体的例子，一维热传导方程的求解。
3.1  一维热传导方程的求解
一维热传导方程是一个典型的抛物型二阶偏微分方程。设
 表示在时间 
 ，空间 
 处的温度，则根据傅里叶定律(单位时间内流经单位面积的热量和该处温度的负梯度成正比)，可以导出热传导方程【2】
其中
 称为热扩散率， 
 分别为热导率，比热和质量密度，都是由系统本身确定的常量。
为了具体，设
 ，设边界条件为
设步长为：
 ，从而 
 ，所以递推关系为
图2：递推关系示意图
图2直观地给出了差分法求解偏微分方程的过程。先把时空坐标都离散化，根据递推关系，由下一行的三个蓝点的值可以给出上一行的一个红点的值，由于边界条件和初始条件(即最下方和两边的绿线)已知，所以按这个递推关系可以得到网格中的所有值。下面我们用Python代码来实现。
import numpy as np
import matplotlib.pyplot as plt
h = 0.1#空间步长
N =30#空间步数
dt = 0.0001#时间步长
M = 10000#时间的步数
A = dt/(h**2) #lambda*tau/h^2
U = zeros([N+1,M+1])#建立二维空数组
Space = arange(0,(N+1)*h,h)#建立空间等差数列，从0到3，公差是h
#边界条件
for k in arange(0,M+1):
U[0,k] = 0.0
U[N,k] = 0.0
#初始条件
for i in arange(0,N):
U[i,0]=4*i*h*(3-i*h)
#递推关系
for k in arange(0,M):
for i in arange(1,N):
U[i,k+1]=A*U[i+1,k]+(1-2*A)*U[i,k]+A*U[i-1,k]
上述代码中，我们首先把位于0-3中的空间等分为30份，位于0-1的时间等分为10000份，然后通过设置初始条件、边界条件和递推关系并借助for循环就得到了1个30*10000的二维数组，里面放着每个离散的时空点的温度值。
为了直观地展现温度随时空的变化关系，接下来开始画图，首先画出不同时刻温度随空间坐标的变化
#不同时刻的温度随空间坐标的变化
plt.plot(Space,U[:,0], 'g-', label='t=0',linewidth=1.0)
plt.plot(Space,U[:,3000], 'b-', label='t=3/10',linewidth=1.0)
plt.plot(Space,U[:,6000], 'k-', label='t=6/10',linewidth=1.0)
plt.plot(Space,U[:,9000], 'r-', label='t=9/10',linewidth=1.0)
plt.plot(Space,U[:,10000], 'y-', label='t=1',linewidth=1.0)
plt.ylabel('u(x,t)', fontsize=20)
plt.xlabel('x', fontsize=20)
plt.xlim(0,3)
plt.ylim(-2,10)
plt.legend(loc='upper right')
plt.show()图3：不同时刻的温度随空间坐标的变化
从图3可以看到，温度关于
 呈现轴对称分布，这是由初始条件造成的。另外，对每一点的空间坐标，随着时间的推移，温度越来越低。
接下来，我们来画出温度等高线来描述温度随任意时空点的变化
#温度等高线随时空坐标的变化，温度越高，颜色越偏红
extent = [0,1,0,3]#时间和空间的取值范围
levels = arange(0,10,0.1)#温度等高线的变化范围0-10，变化间隔为0.1
plt.contourf(U,levels,origin='lower',extent=extent,cmap=plt.cm.jet)
plt.ylabel('x', fontsize=20)
plt.xlabel('t', fontsize=20)
plt.show()图4：温度随时空坐标的变化，温度越高，颜色越红
在图4中，利用颜色的深浅来标记温度，温度越高，颜色越红。从中同样可以看到，温度随空间的分布关于
 轴对称，而且随着时间的推移，温度越来越低。
4  总结
在本文中，我们利用Python数值求解了常微分方程和偏微分方程，基本思想是先将连续的坐标离散化，然后将微分化为差分，由差分方程得到递推关系，然后利用计算机强大的重复计算能力得到任意格点处的函数值。虽然上面只算了两个例子，但是这种方法完全可以推广到任意偏微分方程的求解中。
在量子色动力学(QCD)中，由于强相互作用具有渐进自由的特性，所以在低能情况下没办法像QED那样使用微扰论计算，这时就要采用格点QCD的方法计算。其基本思想也是将时空离散化，然后从第一性原理的路径积分出发去计算。由于时空被离散化了，相当于人为地引入了一个最小时空距离
 ，在傅里叶变换到动量空间后相当于引入了一个最大的动量截断 
 ，所以计算结果不会出现紫外发散，从而可以算到很高的精度，在一些情况下，格点的计算结果甚至比实验更精确。所以，将连续参数离散化，把微分化为差分的思想，是极其重要的。
【1】不同的算法对于方程右边具体取什么形式并不一样，从而精度也不一样。例如，欧拉法右边取得是
 ；改进的欧拉法右边取的是 
 ；二阶Runge-Kutta法右边取的是 
 ；四阶Runge-Kutta法右边取的是 
 ，其中 
， 
 ， 
 ， 
 。这些算法的差别在于计算精度不同，并不改变差分的本质思想。为了具体，我们这里采用改进的欧拉法。
【2】傅里叶定律告诉我们单位时间通过单位面积的热量和该处的温度负梯度成正比，即
 ，其中 
 是热流，即单位时间通过单位面积的热量， 
 为热导率， 
 为温度。能量守恒定律告诉我们单位时间流出某闭合曲面的热量等于其内部减少的能量，即
 ，其中 
 表示单位体积的热容，而 
 和 
 分别为质量密度和单位质量的热容(比热)，联合散度定理 
 ，我们得到 
 ，再将傅里叶定律带入，就得到了热传导方程 
 ，其中 
 称为热扩散率。在一维的情况下，热传导方程就退化到了正文中的形式，即 
 。

fipy
菲克定律是指在不依靠宏观的混合作用发生的传质现象时，描述分子扩散过程中传质通量与浓度梯度之间关系的定律。菲克定律是阿道夫·菲克(Adolf Fick)于1855年提出。
由菲克第二定律可以得到动态扩散的偏微分方程。求解可以得到浓度分布和流出曲线。
不确定这个问题有没有解析解，不过数值求解是一种较为通用的解决方法。
fipy是目前难得的还活着的PDE求解python包，作者根据官方示例改写本程序。
相关推荐：《python视频教程》
问题
一个二维平板，顶端1摄氏度(100也可以，只是一个系数)，另外三个边缘0摄氏度，初始时刻整个板子都是0摄氏度，随之时间的推进，热量在板子上传递，最后达到平衡态，我们不仅希望知道平衡态的温度分布，也希望知道温度随时间是如何变化的。热量的传递由微分方程给出，可以简单地理解为热量按照温度降低最快的方向进行传递。
公式右边是温度的梯度，左边是温度随时间的变化
最后整个板子的温度分布大致呈现怎样
只有一个包需要导入import fipy as fp
确定求解区域，一个20*20的格点
#求解区域nx = 20ny = 20dx = 1.dy = dxL = dx * nxmesh = fp.Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
phi = fp.CellVariable(name = "solution variable",                      mesh = mesh,                      value = 0.)
创建微分方程
设立边界条件
创建画图
求解
建议在命令行里面运行，命令行里面可以获得动图，ipython里面只有最后一张图
spyder的ipython里面只有最后的一张图片

	

附：
微分方程的分类：常微分方程和偏微分方程。
1、常微分方程(ODE)是指微分方程的自变量只有一个的方程。最简单的常微分方程，未知数是一个实数或是复数的函数，但未知数也可能是一个向量函数或是矩阵函数，后者可对应一个由常微分方程组成的系统。
2、偏微分方程(PDE)是指微分方程的自变量有两个或以上，且方程式中有未知数对自变量的偏微分。偏微分方程的阶数定义类似常微分方程，但更细分为椭圆型、双曲线型及抛物线型的偏微分方程，尤其在二阶偏微分方程中上述的分类更是重要。
有些偏微分方程在整个自变量的值域中无法归类在上述任何一种型式中，这种偏微分方程则称为混合型。
常微分方程及偏微分方程都可以分为线性微分方程及非线性微分方程二类。
一个常微分方程是不是有特解呢？如果有，又有几个呢？这是微分方程论中一个基本的问题，数学家把它归纳成基本定理，叫做存在和唯一性定理。因为如果没有解，而我们要去求解，那是没有意义的；如果有解而又不是唯一的，那又不好确定。因此，存在和唯一性定理对于微分方程的求解是十分重要的。