有限元求解两点边值问题之二
问题:$$-\frac{d^2u}{dx^2}=f(x)$$
边值条件为$$u(a)=0,u(b)=1$$,
精确解为$$u=x^{\frac{3}{4}}$$
其Galerkin推导有限元方程的过程可参见《微分方程数值解法》李荣华第四版P226
代码如下
# -*- coding: utf-8 -*-
"""
Created on Wed Jun 6 21:04:25 2018@author: shaowuGalerkin有限元法求解两边值问题,方程形式:
p*u''=f(u,x),u(a)=0,u(b)=1
其双线性形式为a(u,v)=p*u'*v'的一重积分
"""
import numpy as np
import pandas as pd
import scipy as sp
import numpy.linalg as LA
import matplotlib.pyplot as plt
from scipy.special.orthogonal import p_roots
import time
start_time=time.time()
def gauss_xw(m=100):"""默认用100个点求高斯——勒让德节点point和权weight,并返回"""point,weight=p_roots(m+1)return point,weightdef function_u(x):"""精确解"""return x**(3/4)def f_function(x):"""定义方程右端函数"""return 3/16*x**(-5/4)def hat_function(x,t,i,h,n):#构造试探函数空间的一组基函数:线性插值公式构造如下#参数:h 步长列表#t 节点列表#i 第i个节点#n 节点个数#x 变量if i==n:if t[n]>=x and x>=t[n-1]:return 1+(x-t[n])/h[n-1]else:return 0else:if t[i]>=x and x>=t[i-1]:return 1+(x-t[i])/h[i-1]elif t[i+1]>=x and x>=t[i]:return 1-(x-t[i])/h[i]else:return 0
"""
def hat_function(x,t,j,h,n):#h=(b-a)/(n)#t=[a+i*h for i in range(0,n+1)]if j==0:if x>=t[0] and x<=t[1]:return 1-abs(x-t[0])/h[j]else:return 0if j==n:if t[n-1]<=x and x<=t[n]:return 1-abs(x-t[n])/h[j]else:return 0if j>0 and j<n:if t[j-1]<=x and x<=t[j+1]:return 1-abs(x-t[j])/h[j]else:return 0
"""
def compute_right_of_equation(point,weight,t,h,n,a,b):"""计算方程的右端项参数:h 步长列表a,b 对应区间[a,b]"""c=(b-a)/2s=(b-a)/2*point+(a+b)/2 #把区间[a,b]变化到[-1,1]f=[sum([c*weight[i]*f_function(s[i])*hat_function(s[i],t,j,h,n) for i in range(len(s))])for j in range(1,len(t)-1)]f[-1]=f[-1]+1/h[-1] #添加右端初值条件return f
def compute_stiffness_matrix(point,weight,t,h,n,a,b):"""计算系数矩阵,即总的刚度矩阵"""A=np.zeros((n-1,n-1))for j in range(1,len(t)-1):if j==1:A[j-1,j-1]=+sp.integrate.quad(lambda x: h[j]**(-2),t[j-1],t[j+1])[0]#sp.integrate.quad(lambda x: h[j]**(-2),t[j],t[j+1])[0]A[j-1,j]=sp.integrate.quad(lambda x: -h[j]**(-2),t[j],t[j+1])[0]#A[j,j-1]=sp.integrate.quad(lambda x: h[j]**(-2),t[j-1],t[j])[0]elif j==n-1:A[j-1,j-2]=sp.integrate.quad(lambda x: -h[j-1]**(-2),t[j-1],t[j])[0]A[j-1,j-1]=+sp.integrate.quad(lambda x: h[j-1]**(-2),t[j-1],t[j+1])[0]else:A[j-1,j-2]=sp.integrate.quad(lambda x: -h[j-1]**(-2),t[j-1],t[j])[0]A[j-1,j-1]=+sp.integrate.quad(lambda x: h[j-1]**(-2),t[j-1],t[j])[0]+\sp.integrate.quad(lambda x: h[j]**(-2),t[j],t[j+1])[0]A[j-1,j]=sp.integrate.quad(lambda x: -h[j]**(-2),t[j],t[j+1])[0]return np.array(A)def solve_ui(A,b):"""计算ui"""return np.linalg.solve(A,np.array(b))
def solve_uh(ui,t,h,n,a,b):"""计算uh和误差"""uh=[sum([ui[j]*hat_function(x,t,j+1,h,n) for j in range(n-1)]) for x in t if x!=1]error=function_u(np.array(t[:-1]))-uhreturn uh,LA.norm(error,np.inf)
if __name__ == '__main__':a=0 #区间a,bb=1n=20h=(b-a)/(n) #等距步长t=[a+i*h for i in range(0,n+1)] #节点h=[(b-a)/(n) for i in range(n)] #子区间步长集合point,weight=gauss_xw(m=100) #高斯勒让德节点point和权weightf=compute_right_of_equation(point,weight,t,h,n,a,b) #计算方程右端项A=compute_stiffness_matrix(point,weight,t,h,n,a,b) #计算系数矩阵,即总刚度矩阵ui=solve_ui(A,f) #计算uiuh,error=solve_uh(ui,t,h,n,a,b) #计算uh和误差print('the error is:',error)plt.plot(t,[0]+list(ui)+[1],'+b',label='approximate solution')plt.plot(t,function_u(np.array(t)),'r',label='exact solution')plt.title('Fig2.(where a=0,b=1,n=20,$u(x)=x^{3/4}$,$u(1)=0,u(1)=1$)')plt.xlabel('x')plt.ylabel('u')plt.legend(loc='upper center')print('all cost time:',time.time()-start_time)
运行结果
the error is: 4.12305198818e-05
all cost time: 0.23321843147277832
有限元求解两点边值问题之二相关推荐
- 有限元求解两点边值问题之一
问题一: $$\frac{d^2u}{dx^2}=f(x)$$ 边界条件:$$u(a)=u(b)=0$$, 精确解为:$$u=x^4-1$$ 运用Galerkin方法建立有限元方程(其过程可参考< ...
- python有限元传热求解_二维稳态热传导基本方程的有限元求解(2)
四节点矩形单元 在二维稳态热传导基本方程的有限元求解(1)这篇文章中,我们仅仅给出了有限元单元方程的一种比较标准的推导步骤,并未涉及某种具体的单元.且在式(20)中,单元 上温度 的近似函数表示成节点 ...
- 【微分方程数值解】有限差分法(二)两点边值问题数值算例(附python代码)
1.两点边值问题形式 一般的两点边值问题形式为: −uxx=f(x)-u_{xx}=f(x)−uxx=f(x) u(a)=α,u(b)=βu(a)=\alpha,u(b)=\betau(a)=α,u ...
- matlab 读取图片后分区域编号_你的第一个有限元求解器——仅十行MATLAB代码
有限元分析话题中有不少讨论有限元求解器的问题,但大都停留在概念层面,未见实际代码.望本文能略起抛砖引玉之作用. 以下代码是基于MATLAB编写. 问题描述 考虑一平面有界区域 ,设其边界为 .我们求解 ...
- 有限差分法matlab两点边值代码,两点边值问题的有限差分法.doc
两点边值问题的有限差分法.doc 学生实验报告实验课程名称偏微分方程数值解开课实验室数统学院学院数统年级2013专业班信计2班学生姓名学号开课时间2015至2016学年第2学期总成绩教师签名数学与统计 ...
- 洛谷题库P5735距离函数C语言,扩展有限元求解弱不连续问题..docx
扩展有限元求解弱不连续问题. 扩展有限元求解弱不连续问题数学计算机学院 信息与计算科学专业 2011届 段斯琦摘要:扩展有限元 ( XFEM ) 是在标准有限元的框架下,提出来的一种有利于解决裂纹.孔 ...
- 【CAE】优秀的开源有限元求解器
工业有限元求解器是工业仿真中最重要,最核心的,目前市场上有较多的开源求解器,下面列举几个比较有名的开源求解器,仅供各位对国产CAE感兴趣的朋友参考.借鉴,学习CAE软件开发的框架,思路等. 1. Op ...
- EIT正问题求解--利用有限元求解节点电势
为了理解EIT正问题利用有限元求解的方法,自己通过在网上查找程序,结合EIDORS软件包所建立的EIT模型,编写了求解节点电势的matlab程序,并且进行了验证. 1.利用EIDORE建立模型,建立模 ...
- matlab微分方程组边值,matlab求解常微分方程边值问题的方法
matlab求解常微分方程边值问题的方法 Matlab 求解常微分方程边值问题的方法:bvp4c 函数 常微分方程的边值问题,即 boundary value problems ,简称 BVP 问题, ...
最新文章
- 南邮 AAencode
- STL源代码分析(ch2 内存分配)uninitialized_fill_n
- MyBatis实现SaveOrUpdate
- 产品经理必须知道的一点知识:三种方法判断一个产品该不该做
- ENSP配置 实例五 RIP配置
- android全局计时_Android定时器AlarmManager
- ZZULIOJ 1158: 又是排序(指针专题)
- 高手如何应对复杂系统架构的演进
- 下一个技术之城:长沙
- .5-浅析express源码之Router模块(1)-默认中间件
- Epub,Mobi,Azw3电子书格式的区别,有什么好用的安卓epub阅读器
- 关于学历与面试的一些看法
- [MICCAI2019] A Partially Reversible U-Net for Memory-Efficient Volumetric Image Segmentation
- java contains 大小写_用.contains方法忽略大小写的选项?
- SQL基础语法_刘世民
- 一、系统间的通信技术
- gdiplus图像库的使用
- 荣耀20公测鸿蒙,荣耀 20、30 系列等机型,将开始逐步适配华为鸿蒙系统
- [附源码]java毕业设计房屋租赁管理系统
- pytorch tensor求向量的模长