写点什么吧，希望坚持

写点什么吧，希望坚持

---

一点随想吧

• 一直想着坚持谢谢自己的博客，但是种种原因却没坚持下来，今天突然想起，对自己没能坚持感到很遗憾。之前一直希望自己每日能落笔千字，这样在博士毕业时，也能有一定的积累了，可惜的是没能坚持住。
• 每天写点什么呢？除了定期会更新自己对于计算和编程实现方面学习心得，还想更新一些自己每天的收获吧，我想记录一下自己每天都在干什么，害怕自己荒废时光。
• 最近一直在写开题报告，写写自己的一些想法和思考吧，为什么开题报告这么重要？如果对于所要研究的问题没有一个从现象到后背科学规律发掘的过程，那么研究就一定不会很全面、客观地反映出问题，清晰的研究构思，不是指必须准确、不能改变的定下研究要走的路，而是展现出对于问题思考的过程，随着研究的进行很可能出现我们预想不到的问题，但是对于背后科学规律的把握不能在一开始就出现逻辑上的错误。因此，从撰写开题报告来看，实际上是在完成一个认识、思考的过程，首先，在背景、意义中交代清楚目前出现实际问题，这也是所要研究问题的基础，以及价值。在研究目标中提出自己所做的工作能为解决、客服目前的问题带来哪些帮助。研究内容就是要写出为了完成、达到研究目标，我需要拿哪些工作，当然对于每一项工作，一定要写出我为什要做它！而不是做别的！在后续的研究方法中，要针对前面提出的研究内容，做一点怎么进行的思考工作，大体上写出这个内容如何进行思考，解决。我觉上述这些也就是开题报告撰写的思路和核心所在。

python解决一个抛物线PDE方程

之前一直用Matlab搞一点数值计算，后来发现了python，自己比较喜欢，所以也在尝试着用python解决问题，下面给一个我自己写的一个python解决扩散方程的代码。实际上，抛物线方程一般涉及到的问题就是热传导问题（非稳态），以及渗流固结。我主要解决固结理论，关于一维的太沙基固结控制方程的推到就不做说明了，这一点应该每个专业人士都掌握了，重点在于PDE的数值解。
* 问题的描述
u的一阶导数采用两点向前差分表示，u的二阶导数利用三点中心差分表示，固结方程如下：

$$C_v \cdot \dfrac{\partial ^2u}{\partial x^2} = \dfrac{\partial u}{\partial t}$$

边界条件：

$$u(0,t) = 0 \ kPa\\ u(1,t),x = 0 \ kPa$$
初始条件：
$$u(x,0) = 40 \ kPa$$

差分法：场函数u(i,j) ， x轴为i控制，y轴为t，采用j控制
对于t的一阶导数表示为：

$$\dfrac{\partial u}{\partial t} = \dfrac{u(i, j+1) - u(i, j)}{\Delta t}$$

$$\dfrac{\partial ^2 u}{\partial x^2} = \dfrac{u(i-1, j)-2 \cdot u(i, j) + u(i+1, j)}{\Delta x^2}$$
带入原始的抛物线方程后，转化为差分方程.
$$r = \dfrac{C_v \cdot \Delta t }{\Delta x^2}$$
为了简化差分的计算，采用r=0.5，则若 $\Delta x$确定，那么 $\Delta t$自然确定，得到最终的差分方程为：
$$u(i, j+1) = 0.5 \times (u(i-1, j) + u(i+1, j))$$
编程求解思路

首先对于场函数来说，必须确定计算区域，即按照实际需求确定一个t的范围，这样整个区域为一个矩形，然后进行差分点的确定，主要确定边的份数，为了使得差分公式第一次可以使用，必须要先处理整个区域的边界条件，即这个矩阵的周围点的值，然后才可以考虑差分计算中间点。

import numpy as np
import matplotlib.pyplot as plt

#初始参数的自定义
cv = 3e-6    #固结系数值定义
xL = 1    #x方程的长度
xn = 20   #x方向划分的份数
T = 5e6     #t总时间

u0 = 40000    #初始孔压值，这里先确定为一个均匀分布的值40kPa

# 由初始信息计算必须要的基础量
dx = xL/xn
x = np.arange(0,xL+dx,dx)
dt = 0.5*dx**2/cv
t = np.arange(0,T+dt,dt)
tn = t.size - 1    # t.size返回的是t轴方向存在的节点个数，因此t方向的分数需要减去1
u = np.zeros((tn+1,xn+1))    #初始化一个u矩阵，这个矩阵是依据节点个数建立的，

#处理边界条件
#由于 左边界和右边界均为0，zeros函数直接定义好了为0，即只做底边界的处理即可,注意边界的交点处取1/2即可

u[0,:] = u0*np.ones(xn+1)
u[0,0] = 0.5*u0
u[0,xn] = 0.5*u0

#计算u中的自由点

for rows in np.arange(tn-1):

    for cols in np.arange(1,xn):

        u[rows+1,cols] = 0.5 * ( u[rows,cols-1] + u[rows, cols+1] )

#解答的可视化

plt.figure(1)
cf1 = plt.contourf(x, t, u)
#plt.clabel(cf1)
plt.xlabel('x/ m')
plt.ylabel('T/ second')
plt.title(' one - dimmension consolidation equation')
plt.axis([0,1,0,150000])
plt.colorbar(cf1)
plt.savefig('contourf image.png',dpi = 150)
plt.show()

plt.figure(2)
cf2 = plt.contour(x, t, u)
plt.clabel(cf2, inline = 1.0, fontsize = 15)
plt.xlabel('x/ m')
plt.ylabel('T/ second')
plt.title(' one - dimmension consolidation equation')
plt.axis([0,1,0,150000])
plt.colorbar(cf2)
plt.savefig('contour image.png',dpi = 150)
plt.show()


plt.figure(3)
plt.plot(x,u[0,:], ls = '-', marker = '*', lw = 2, color = 'r', label = 'T_step=0' )
plt.plot(x,u[50,:], ls = '-', marker = 'o', lw = 2, color = 'y', label = 'T_step=50' )
plt.plot(x,u[100,:], ls = '-', marker = 's', lw = 2, color = 'm' , label = 'T_step=100')
plt.plot(x,u[150,:], ls = '-', marker = '^', lw = 2, color = 'c' , label = 'T_step=150')
plt.plot(x,u[200,:], ls = '-', marker = '*', lw = 2, color = 'k', label = 'T_step=200' )
plt.plot(x,u[250,:], ls = '-', marker = 'o', lw = 2, color = 'b', label = 'T_step=250' )
plt.plot(x,u[300,:], ls = '-', marker = 's', lw = 2, color = 'g' , label = 'T_step=300')
plt.plot(x,u[350,:], ls = '-', marker = '^', lw = 2, color = 'tomato' , label = 'T_step=350')
plt.title('curve of u with T_step')
plt.legend()
plt.savefig('curve1 image.png',dpi = 150)
plt.show()

plt.figure(3)
plt.plot(t,u[:,0],ls = '--', marker = 'o', lw = 1, color= 'c', label = 'x=0.0')
plt.plot(t,u[:,2],ls = '--', marker = 'o', lw = 1, color= 'k', label = 'x=0.1')
plt.plot(t,u[:,4],ls = '--', marker = 'o', lw = 1, color= 'r', label = 'x=0.2')
plt.plot(t,u[:,6],ls = '--', marker = 'o', lw = 1, color= 'g', label = 'x=0.3')
plt.plot(t,u[:,8],ls = '--', marker = 'o', lw = 1, color= 'b', label = 'x=0.4')
plt.plot(t,u[:,10],ls = '--', marker = 'o', lw = 1, color= 'y', label = 'x=0.5')
plt.plot(t,u[:,12],ls = '--', marker = 'o', lw = 1, color= 'm', label = 'x=0.6')
plt.plot(t,u[:,14],ls = '--', marker = 'o', lw = 1, color= 'orange', label = 'x=0.7')
plt.plot(t,u[:,16],ls = '--', marker = 'o', lw = 1, color= 'pink', label = 'x=0.8')
plt.plot(t,u[:,18],ls = '--', marker = 'o', lw = 1, color= 'tomato', label = 'x=0.9')
plt.plot(t,u[:,20],ls = '--', marker = 'o', lw = 1, color= 'c', label = 'x=1.0')
plt.title('curve of u with distance')
plt.axis([0, 150000, 0, 41000])
plt.legend()
plt.savefig('curve2 image.png',dpi = 150)
plt.show()