用python怎样解偏微分方程组_用Python数值求解偏微分方程
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】傅里叶定律告诉我们单位时间通过单位面积的热量和该处的温度负梯度成正比,即
,其中
是热流,即单位时间通过单位面积的热量,
为热导率,
为温度。能量守恒定律告诉我们单位时间流出某闭合曲面的热量等于其内部减少的能量,即
,其中
表示单位体积的热容,而
和
分别为质量密度和单位质量的热容(比热),联合散度定理
,我们得到
,再将傅里叶定律带入,就得到了热传导方程
,其中
称为热扩散率。在一维的情况下,热传导方程就退化到了正文中的形式,即
。
用python怎样解偏微分方程组_用Python数值求解偏微分方程相关推荐
- python牛顿法解非线性方程组_利用python求非线性方程
最近在做的东西中有一件任务,相当于一个函数已知y来求x,网上找了各种办法最终得以实现.在此说明方法,并记录一些坑. 要求的函数比如:*log(x) - log(1-x) + 2.2 * (1 -2x) ...
- python牛顿法解非线性方程组_用牛顿迭代法解非线性方程组
题目: 用牛顿迭代法解非线性方程组 有两个非线性方程,未知数是x1,x2: (15x1+10x2)/[(40-30x1-10x2)^2×(15-15x1)]=5e-4; (15x1+10x2)/[(4 ...
- 用python怎样解偏微分方程组_Python能解偏微分方程吗
fipy 菲克定律是指在不依靠宏观的混合作用发生的传质现象时,描述分子扩散过程中传质通量与浓度梯度之间关系的定律.菲克定律是阿道夫·菲克(Adolf Fick)于1855年提出. 由菲克第二定律可以得 ...
- python牛顿法解非线性方程组_萌新请教牛顿法求解三元非线性方程组
该楼层疑似违规已被系统折叠 隐藏此楼查看此楼 请问牛顿迭代法求解三元非线性方程组,不同迭代初值貌似有很多不同结果,如何求解到满足要求的解,0 FindRoot[{-6.565159793181527` ...
- python牛顿法解非线性方程组_科学网—求解多元非线性方程组F(x)=0的Newton-Raphson方法及其MATLAB实现 - 王福昌的博文...
科学网对公式支持不太好,在博客园有相同博文 牛顿迭代法可以推广到多元非线性方程组 $boldsymbol{F}(boldsymbol{x})=boldsymbol{0}$的情况,称为牛顿-- 拉夫逊方 ...
- python牛顿法解非线性方程组_牛顿迭代法解非线性方程组(MATLAB版)
牛顿迭代法,又名切线法,这里不详细介绍,简单说明每一次牛顿迭代的运算:首先将各个方程式在一个根的估计值处线性化(泰勒展开式忽略高阶余项),然后求解线性化后的方程组,最后再更新根的估计值.下面以求解最简 ...
- python制作解压工具_使用python制作一个解压缩软件
python实现解压缩的重要模块就是--zipfile,其次是os 安装zipfile模块 首先得安装zipfile模块,打开cmd输入一下命令即可安装 pip install zipfile os是 ...
- python牛顿法解非线性方程组_matlab实现牛顿迭代法求解非线性方程组.pdf
matlab实现牛顿迭代法求解非线性方程组.pdf matlab 实现牛顿迭代法求解非线性方程组实现牛顿迭代法求解非线性方程组 已知非线性方程组如下 3*x1-cosx2*x3-1/20 x12-81 ...
- python input与返回值-Python 详解基本语法_函数_返回值
Python 详解基本语法 概要: 函数的返回值是函数重要的组成部分.函数的根本在于实现程序的部分功能,所以很多时候我们需要将函数执行后的结果返回给程序再由程序作出进一步的操作.可以说是函数的返回值令 ...
- python解初中题_用python解一道数独小题
个人的第一篇博文,还请多多支持,不当之处,还请多多指教.(以后有精力还会写更多的文章) 本人是一名大一狗,目前为止学了半年python,对python也就有一点点的了解,没事爱编写一些小程序玩,不过都 ...
最新文章
- c++ 11 锁_国民技术面向智能锁市场提供全系芯片与开源安全解决方案
- Android开发技术周报 Issue#81
- GridView常用总结
- 【学习笔记】JS进阶语法一事件基础
- php 命令执行脚本文件路径,php命令行(cli)下执行PHP脚本文件的相对路径的问题解决方法...
- AIX卷管理介绍以及利用空闲PP来创建文件系统
- mysql 查找配置文件 my.ini 位置方法
- 微软公布 Visual Studio 2020 上半年路线图
- android 图片墙拼贴,三步搞定 用APP打造图片文字拼贴效果
- Ubuntu 8.04 Hardy LTS 软件源设置
- 2019年计算机学业水平测试填空题,2019年计算机学业水平模拟测试选择题80题Word(含参考答案)...
- php课设报告致谢_奇安信CERT发布1月安全监测报告:需警惕这19个高危漏洞
- 英文XP系统安装中文包
- ROS操作系统入门学习
- telnet 测试IP和端口命令
- ps盖印图层在哪里_PS盖印图层快捷键
- 小米WatchS2和小米WatchS1 区别 哪个值得入手
- 能被2、3、4、5、6、7、8、9等数整除的数的特征
- HDLBITS笔记29:移位寄存器(包括4位移位寄存器,创建100位左/右旋转器,算术偏移,线性反馈移位寄存器等)
- Llinux装逼命令大全