python简单的计算方法_用python实现简单的有限元方法(二)
华中师范大学 hahakity
第一部分介绍了加权残差法(Weighted Residual Method)并引出有限元算法中经常出现的伽辽金(Galerkin)法。简单说明了有限元与有限差分法的区别:有限差分是对微分算子做差分近似,有限元是对待求函数做基函数展开近似。这一节继续介绍伽辽金有限元算法。
这里特意强调伽辽金有限元算法,是因为还有另一种基于泛函和变分的 Ritz 有限元算法。跳过 Ritz 变分法,直接讲解伽辽金法入门比较简单。
学习目标Galerkin 有限元法
分片基底函数展开
二阶微分算子的降阶
一维泊松方程:已知电荷密度求电势
预备知识加权残差法 (Weighted Residual Method)
Galerkin 有限元法
假设要求解的微分方程有统一的格式,
其中
是微分算子,
是要求解的场,
是源。
将
做函数基底展开 (使用爱因斯坦求和规则,省略
求和符号),
注意此时
为基底函数,人工设定,为已知量,
为待定系数。
选取基底函数
作为伽辽金法中的检验函数,则对第 i 个检验函数,加权残差为,
其中
是全域,
,
可以将
写成矩阵乘积形式,此处用
表示列向量,
是
构成的矩阵,
是
构成的列向量,
是
构成的列向量。
Galerkin 有限元法令残差
为0,通过解
,得到待定系数
。
分片基底函数展开
将全域离散化为 n 个子域,每个子域称作一个单元(“element”)。连续的函数
近似为有限个单元上的分片函数。如果单元很小,每个单元内
近似线性变化。
先考虑简单的一维问题,将区间分成 n 等份,每份长
。
对第 e 个单元区域
,线性插值函数为
其中
,
当
时,
,
当
时,
,
根据
的定义,得到第 j 个基底函数(注意它横跨两个单元),
的函数图像就是下图中间那个大三角形,覆盖第 j 和第 j+1 个单元。它跟
在第 j 个单元有重叠,跟
在第 j+1 个单元有重叠。此时内积都是局部进行,不重叠区域积分为0。
根据
可知第 j 行只有
,
,
三个非零矩阵元,K 矩阵是稀疏的 3 对角阵。
二阶微分算子降阶
算
时,很快遇到第一个问题,线性插值函数一阶微分在每个单元内为常数,
二阶或更高阶导数等于0。如果微分算子里有二阶导数,可以用分部积分法降阶,
二阶微分算子降阶之后变为,
只与全域边界条件有关,对于大部分不在边界上的点,这项为0。刚好在边界上的那些点,
或它的导数已知,因此可以将这一项单独列出来,令
在新的定义下,方程变为,
多出来的那项
是分部积分项在边界处的贡献
构成的列向量。
此时矩阵
第 j 行的非零矩阵元为,
每个矩阵元的积分区间由前面的示意图给出。
注意此时
或者等于 1/h,或者等于 -1/h,积分很好算。
简单的泊松方程求解
我们已经有了足够的知识,从头构造一个简单的一维 Galerkin 有限元算法。这里先演示一维有限元算法。步骤总结如下,将全域离散化为有限个子域(又称单元)
选择基底函数,对每个子域插值
用 Galerkin 法改写微分方程,函数内积计算
的矩阵元
解
这里拿一个简单的一维泊松方程演示伽辽金有限元算法。
物理问题:根据 Maxwell 方程,电场的散度正比于电荷密度
,而电场本身又可写成标量势
的负梯度
, 其中
。将第二个方程代入第一个得到给定电荷分布下静电势满足的柏松方程,
,二维情况下方程为,
,
一维情况下
其中
是真空电极化常数。
选择求解区域
,边界条件
。
选择一个简单的电荷密度分布,
此时,泊松方程有简单的解析解:
, 我们先假装不知道。
计算
的矩阵元,
这个积分用手积挺容易出错,写一小段 python 代码,将积分区域划分为
与
。在 jupyter-notebook 里运行,
import sympy
xjm1, xj, xjp1 = sympy.symbols(['x_{j-1}', 'x_{j}', 'x_{j+1}'])
x, h, eps = sympy.symbols(['x', 'h', '\epsilon'])
rho = lambda x: sympy.sin(sympy.pi * x)
A = - sympy.integrate(rho(x) * (x - xjm1)/h, (x, xjm1, xj))
B = - sympy.integrate(rho(x) * (xjp1 - x)/h, (x, xj, xjp1))
sympy.simplify(A + B)
输出结果:
Jackson 电动力学书本中介绍了一种简化方法,假设
在每个单元里是常数,只对
积分,则积分结果为,
先不管
, 矩阵方程的解为,
。写一段简单的 python 代码,构造矩阵
和列向量
,使用 scipy.sparse.linalg.inv 函数求
的逆,进而得到结果。
from scipy.sparse import dia_matrix
from scipy.sparse.linalg import inv
from numpy import pi
class FEM:
def __init__(self, nodes, xmin=0, xmax=1):
self.nodes = nodes
x = np.linspace(xmin, xmax, nodes)
self.x = x
self.h = x[1] - x[0]
def Kmatrix(self):
n = self.nodes
m = 1/self.h * np.ones(n)
data = [m, -2*m, m]
offsets = [-1, 0, 1]
# 使用 scipy.sparse 稀疏矩阵库,构造 3 对角稀疏矩阵 K
K = dia_matrix((data, offsets), shape=(n, n)).tocsc()
return K
def bvec(self):
'''假设 rho(x) 在每个单元内值为常数,仅对 u_j(x) 做积分'''
return - np.sin(pi * self.x) * self.h
def solve(self):
K = self.Kmatrix()
b = self.bvec()
return inv(K) * b
def compare(self):
ground_truth = 1/(pi**2) * np.sin(pi * self.x)
fem_res = self.solve()
plt.plot(self.x, ground_truth, label="ground truth")
plt.plot(self.x, fem_res, label="finite element")
plt.title("number of nodes =%s"%self.nodes)
plt.xlabel("x")
plt.ylabel(r"$\phi(x)$")
plt.legend(loc='best')
plt.show()
如果用 1001 个节点,计算结果比较精确,101 或 11 个节点,计算结果误差较大。
fem = FEM(1001)
fem.compare()
二维或高维时的网格生成
上面的一维问题有限元与有限差分算法区别很小,将方程右边的
除到左边,则 K 矩阵变成了二阶导数的差分矩阵。但对于高维且有复杂边界的偏微分问题,有限元算法就显示出它的优势。这一篇变得太长,下一篇尝试构造通用的二维和三维有限元求解器,以及如何在二维圆形区域求解本文例子中的静电势。
为了简化任务,使用 gmsh 软件的 python 库生成有限元分析需要的二维和三维网格。安装方法如下,
pip install gmsh
pip install pygmsh
使用 pygmsh,可以很快生成一个二维的圆形区域网格,
import pygmsh
with pygmsh.geo.Geometry() as geom:
geom.add_circle([0.0, 0.0], 1.0, mesh_size=0.2)
mesh = geom.generate_mesh()
# 将网格文件保存为 vtk 格式
mesh.write("test.vtk")
安装 paraview,打开 test.vtk, 则可以看到如下二维圆形区域的网格图,
总结
介绍了伽辽金有限元算法,一维分片函数展开,二阶微分算子降阶,数值求解了简单的一维泊松方程,介绍了网格生成软件 gmsh。将二维和三维有限元算法的通用求解器推迟到下一节。
警告:本文的求解器非通用求解器,仅作原理展示使用。
参考文献:
【1】Tao Pang,An Introduction to computational physics (second edition)
python简单的计算方法_用python实现简单的有限元方法(二)相关推荐
- python pip国内源_【Python】设置pip源为国内源及简单操作
一.pip国内源镜像: 二.修改源方法: 1.临时修改 可以在使用pip的时候在后面加上-index参数,指定pip源: pip install --index https://pypi.tuna.t ...
- python查询mysql数据库_用python操作mysql数据库(之简单查询操作)
1.mysql安装 此处省略一万字....... 2.pip安装MySQLdb模块 sudo pip install mysql-python 3.简单代码#!/usr/bin/env python ...
- python简单实践作业_【Python】:简单爬虫作业
使用Python编写的图片爬虫作业: #coding=utf-8 import urllib import re def getPage(url): #urllib.urlopen(url[, dat ...
- 基于python的性能测试工具_基于 Python 的性能测试工具 locust 与 LR 的简单对比[转发]...
背景 最近自己开发了一个小的接口,功能测完了,突然想测下性能,原来做性能测试,我一直用的是HP的LoadRunner,前一段时间正好看过locust,想想就用这个来测测性能吧. 由于对LR比较熟,正好 ...
- python简单图形输出_基于 Python Matplotlib 模块的高质量图形输出
Matplotlib 是一个由 John Hunter 等开发的,用以绘制二维图形的 Python 模块.它利用了 Python 下的数值计算模块 Numeric 及 Numarray,克隆了许多 M ...
- 简单叙述python的编程规范_简明 Python 编程规范
注:之前发布一篇<简明 Python 编程规范>(见:http://blog.csdn.net/lanphaday/article/details/2834883),本是我给当时所在的公司 ...
- python 图片拼接成数字_用Python语言对任意图像进行m*n的均匀分块并拼接还原(思路非常清晰,步骤简单)...
目录 1.读取原始图像 2.网格划分,将图像划分为m*n块 2.1分块后图像的存储问题 2.2图像的裁剪 2.3图像长宽的整除问题 方法一:四舍五入法 方法二:图像缩放法 方法三:非均分方法 3.显示 ...
- python实习内容过程_对Python实习的简单分析,也许可以帮到你
昨天对实习僧抓取了90条Python实习的数据,今天做一个简单的分析. 因为实习僧网站上关于Python实习岗位太少,所以以下内容进仅供娱乐. 各个地点的平均日工资,南京竟然比上海还高.回去筛选一看, ...
- 用python爬虫下载视频_使用Python编写简单网络爬虫抓取视频下载资源
我第一次接触爬虫这东西是在今年的5月份,当时写了一个博客搜索引擎,所用到的爬虫也挺智能的,起码比电影来了这个站用到的爬虫水平高多了! 回到用Python写爬虫的话题. Python一直是我主要使用的脚 ...
- 用python编写最简单的记事本_利用Python制作一个“电子记事本”
案例内容 今天的挑战就是写一个"记事本"小程序.程序的功能分为三个部分: 1.把内容记录到文件. 2.显示记录的所有内容. 3.删除不再需要的内容. 正式的"记事本&qu ...
最新文章
- 2022-2028年中国液化石油气(LPG)行业投资分析及前景预测报告
- 重型车辆盲区行为检查Behaviours – Heavy Vehicle Blind Spots
- 产品新人如何快速成长?
- 每周进度条(第10周)
- 【模电】0001 实用运放电路分析
- TestCenter测试管理工具功能详解四(I)
- matlab 图片制作动画制作,MATLAB作图之制作动画:单摆运动仿真
- Kubectl常用命令(三)
- html------轮播图
- ceph-deploy osd activate xxx bluestore ERROR
- 论文对图片与什么要求呢?
- seastar介绍及源码分析
- 监控利器之使用JConsole轻松监控JVM运行情况
- 1230v3配服务器内存性能提升,有点小失望 E3-1280 v3性能实测
- zpl 预览html,在将它发送到Zebra打印机之前,使用.NET WinForm打印预览ZPL II命令
- matlab 拉格朗日LM检验,求问,做LM检验和LR检验的stata命令
- 360计算机面试题,360笔试题目2015
- 【python 字母索引】找到英文句子里面每个单词最后一个字母的索引
- 找不到靶机的IP怎么办,教你几步就会了
- acu风格是什么意思_ACU GENII Ecwcs Gore美军风衣(冲锋衣)测评
热门文章
- 创建AD9361的vivado工程并导入SDK(ZYNQ平台)
- 计算机专业有必要考软考吗,软考初级程序员有用吗_有必要考吗_上学吧
- PAgP协议与LACP协议
- 神经网络学习小记录53——TF2搭建孪生神经网络(Siamese network)比较图片相似性
- 快递100码json
- axure android 原型设计工具,知乎和Quora高分APP原型设计工具推荐
- 【linux内核分析与应用-陈莉君】字符设备驱动
- 微信开发者工具 wxmi修改模版颜色_小白变大师试试免费设计工具:adaptiff
- c语言画bode图程序,根据上位机测得的Bode图的幅频特性,就能确定系统(或环节)的相频特性,试问这在什么系统时才能实现?...
- Oracle索引的建立及优缺点