常微分方程的解析解(方法归纳)以及基于Python的微分方程数值解算例实现 |
您所在的位置:网站首页 › 差分法求解常微分方程 › 常微分方程的解析解(方法归纳)以及基于Python的微分方程数值解算例实现 |
常微分方程的解析解(方法归纳)以及基于Python的微分方程数值解算例实现
本文归纳常见常微分方程的解析解解法以及基于Python的微分方程数值解算例实现。 常微分方程的解析解考虑常微分方程的解析解法,我们一般可以将其归纳为如下几类: 可分离变量的微分方程(一阶) 一阶齐次(非齐次)线性微分方程(一阶) 二阶常系数微分方程(二阶) 高阶常系数微分方程(这类微分方程可以变形成如下形式:
某些方程看似不可分离变量,但是经过换元之后,其实还是可分离变量的,不要被这种方程迷惑。 一阶齐次(非齐次)线性微分方程(一阶)形如
解法: (直接套公式)
伯努利方程
形如
形如
求解此类方程分两步: 求出齐次通解 求出非齐次特解原方程的解=齐次通解+非齐次特解 齐次通解的求法首先假设 情况一:方程有两个不同的实数解 假设两个实数解分别是 对于
若0不是方程特征解 则方程有特解 若0是方程的单特征解 则方程有特解 若0是方程的二重特征解 则方程有特解 其中 若 则方程有特解 若 则方程有特解 若 则方程有特解 若 则方程有特解 若 则方程有特解 其中 形如
求解此类方程分两步: 求出齐次通解 求出非齐次特解原方程的解=齐次通解+非齐次特解 齐次通解的求法(参考二阶常系数微分方程解法) 非齐次特解的求法(参考二阶常系数微分方程解法) 算例考虑带有第三类边界条件的二阶常系数微分方程边值问题
问题一:两点边值问题的解析解 由于此方程是非齐次的,故求解此类方程分两步: 求出齐次通解 求出非齐次特解原方程的解=齐次通解+非齐次特解 齐次通解的求法首先假设 此时方程的齐次通解是
由于 原方程的解=齐次通解+非齐次特解 于是,原方程的全解为
综上所述,原两点边值问题的解为
对一般的二阶微分方程边值问题
把区间 对式中的二阶导数仍用数值微分公式
讨论差分方程组的解是否收敛到原微分方程的解,估计误差. 这里就不再详细介绍. 数值算例考虑带有第三类边界条件的二阶常系数微分方程边值问题
问题二:有限差分方法算出其数值解及误差 对于原问题,取步长h=0.2,用有限差分求其近似解,并将结果与精确解y(x)=-x-1进行比较. 解: 因为
代入上述差分方程组公式得到差分格式
二阶线性常系数微分方程边值问题的差分法Python代码 先以将区间划分为5份为例,求出数值解 import numpy as np import matplotlib.pyplot as plt # 准备工作 def initial(L, R, NS): x = np.linspace(L, R, NS + 1) h = (R - L) / NS return x, h # 右端函数f def f(x): return 3 * x + 1 # 解方程 def solve(NS): # 矩阵A A1 = np.array([-2 * h - 3] + [b] * (NS - 1) + [3 + 2 * h]) # print("A1",A1) A2 = np.array([a] * (NS - 1) + [-4]) # print("A2", A2) A3 = np.array([4] + [c] * (NS - 1)) # print("A3", A3) A4 = np.array([-1] + [0] * (NS - 2)) # print("A4", A4) A5 = np.array([0] * (NS - 2) + [1]) # print("A5", A5) A = np.diag(A1) + np.diag(A2, -1) + np.diag(A3, 1) + np.diag(A4, 2) + np.diag(A5, -2) # print("A", A) # 矩阵b br = 2 * h ** 2 * f(x) + np.array( [2 * h - 2 * h ** 2 * f(x[0])] + [0] * (NS - 1) + [-8 * h + 2 * h ** 2 * f(x[NS])]) uh = np.linalg.solve(A, br) return uh L = 1.0 R = 2.0 NS = 5 x, h = initial(L, R, NS) a = 2 - 2 * h b = -4 - 6 * h ** 2 c = 2 + 2 * h uh = solve(NS) print("h=%f时数值解\n" % h, uh)结果: h=0.200000时数值解 [-1.48151709 -1.56543883 -1.6245972 -1.6531625 -1.64418896 -1.58929215]是不是解出数值解就完事了呢?当然不是。我们可以将问题的差分格式解与问题的真解进行比较,以得到解的可靠性。通过数学计算我们得到问题的真解为 结果: h=0.200000时真解 [-2. -2.2 -2.4 -2.6 -2.8 -3. ] L_∞(最大模)范数下的误差 1.4107078471946164 L_2(平方和)范数下的误差 1.0483548020308784接下来绘图比较 结果: 哎呀呀! 根据输出显示的误差,发现此时数值解的求解效果令人难以接受此处首先认为解析解或者查分格式计算错了,休息了一晚上起来重新推导了一遍后并未发现明显错误. 后来在玩GTA的时候忽然想起导师之前提到过步长会影响数值解稳定性···,于是重新编写程序,修改步长后发现确实如此. 将区间划分为 结果: h=0.000781时数值解 [-1.99798816 -1.99876784 -1.99954751 ... -2.99297729 -2.99375427 -2.99453125] h=0.000781时真解 [-2. -2.00078125 -2.0015625 ... -2.9984375 -2.99921875 -3. ] L_∞(最大模)范数下的误差 0.005468750663166322 L_2(平方和)范数下的误差 0.0035976563635910903绘图比较 最后,我们还可以从数学的角度分析所采用的差分格式的一些性质。因为差分格式的误差为 结果: 网格密度加倍时误差变化情况 步长h=0.200000 最大模误差=1.410708 步长h=0.100000 最大模误差=0.701417 步长h=0.050000 最大模误差=0.350182 步长h=0.025000 最大模误差=0.175023 步长h=0.012500 最大模误差=0.087503 步长h=0.006250 最大模误差=0.043750 步长h=0.003125 最大模误差=0.021875 步长h=0.001563 最大模误差=0.010938 步长h=0.000781 最大模误差=0.005469 步长h=0.000391 最大模误差=0.002734 完整代码 # 开发者: Leo 刘 # 开发环境: macOs Big Sur # 开发时间: 2021/10/5 6:51 下午 # 邮箱 : [email protected] # @Software: PyCharm # ---------------------------------------------------------------------------------------------------------- import numpy as np import matplotlib.pyplot as plt # 准备工作 def initial(L, R, NS): x = np.linspace(L, R, NS + 1) h = (R - L) / NS return x, h # 右端函数f def f(x): return 3 * x + 1 # 解方程 def solve(NS): # 矩阵A A1 = np.array([-2 * h - 3] + [b] * (NS - 1) + [3 + 2 * h]) # print("A1",A1) A2 = np.array([a] * (NS - 1) + [-4]) # print("A2", A2) A3 = np.array([4] + [c] * (NS - 1)) # print("A3", A3) A4 = np.array([-1] + [0] * (NS - 2)) # print("A4", A4) A5 = np.array([0] * (NS - 2) + [1]) # print("A5", A5) A = np.diag(A1) + np.diag(A2, -1) + np.diag(A3, 1) + np.diag(A4, 2) + np.diag(A5, -2) # print("A", A) # 矩阵b br = 2 * h ** 2 * f(x) + np.array( [2 * h - 2 * h ** 2 * f(x[0])] + [0] * (NS - 1) + [-8 * h + 2 * h ** 2 * f(x[NS])]) uh = np.linalg.solve(A, br) return uh L = 1.0 R = 2.0 NS = 1280 x, h = initial(L, R, NS) a = 2 - 2 * h b = -4 - 6 * h ** 2 c = 2 + 2 * h uh = solve(NS) print("h=%f时数值解\n" % h, uh) # 真解函数 def ture_u(x): ture_u = -x - 1 return ture_u t_u = ture_u(x) print("h=%f时真解\n" % h, t_u) # 误差范数 def err(ture_u, uh): ee = max(abs(ture_u - uh)) e0 = np.sqrt(sum((ture_u - uh) ** 2) * h) return ee, e0 ee, e0 = err(t_u, uh) print('L_∞(最大模)范数下的误差', ee) print('L_2(平方和)范数下的误差', e0) # 绘图比较 plt.figure() plt.plot(x, uh, label='Numerical solution') plt.plot(x, t_u, label='Exact solution') plt.title("solution") plt.legend() plt.show() # 误差与网格关系讨论 N = [5, 10, 20, 40, 80, 160, 320, 640, 1280, 2560] print("网格密度加倍时误差变化情况") for i in range(len(N)): NS = N[i] x, h = initial(L, R, NS) a = 2 - 2 * h b = -4 - 6 * h ** 2 c = 2 + 2 * h uh = solve(NS) t_u = ture_u(x) ee, e0 = err(t_u, uh) print('步长h=%f' % h, end='\t') print('最大模误差=%f' % ee)数值算例来源: 《微分方程数值解》-M.Ran |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |