智障。 - Modelica
type: channel
link: http://www.dmcxe.com/index.php/tag/Modelica/
items: 数值计算三——微分方程组的求解与Modelica

数值计算三——微分方程组的求解与Modelica
type: item
date: 2023-04-22T10:40:03+00:00
description: 本次将通过改进RK4（四阶龙格-库塔法）法构建ODE求解器以获得对于微分方程组的求解，并通过相关大作业的实力内容进行计算验证。通过此方程组的建立，对于Modelica运行时的一些event耗时问题也产生了一些思考与理解。此外，代码全部上传至https://github.com/DMCXE/Numerical_ComputationODE.py Homework.moRK4方法的回顾常用的ODE求解方法为龙格-库塔方法。即增加欧拉法中的节点数提高精度。RK方法可以用公式表示为： y n + 1 = y n + h φ ( x n , y n , h ) , 其中 φ ( x n , y n , h ) = ∑ i = 1 r c i K i K 1 = f ( x n , y n ) K i = f ( x n + λ i h , y n + h ∑ j = 1 i − 1 μ i j K j ) , i = 2 , ⋯ , r 当对r取不同的取值，则对应了不同的方法。四阶RK4方法具备较高的精度，常用于数值计算中。通过数学计算，可以得到RK4的一般表达式为： { y n + 1 = y n + h 6 ( K 1 + 2 K 2 + 2 K 3 + K 4 ) , K 1 = f ( x n , y n ) , K 2 = f ( x n + h 2 , y n + h 2 K 1 ) , K 3 = f ( x n + h 2 , y n + h 2 K 2 ) , K 4 = f ( x n + h , y n + h K 3 ) . 但通常微分方程是通过微分方程组的形式出现，比如说传染病模型SEIR。而在《数值计算二》中构建的ODE求解器只能够求解单独的微分方程，局限性较高。开发通用的微分方程组求解器是必须需要解决的。RK4方法对方程组的适应改进在这里，RK4是一种前向方法。因此当函数化为方程组时，只需要将方程指定为向量格式、并在自变量统一时，将因变量写成向量格式。若f(x,y)表示的形式是f(t,x)，即写为f(t,x1,x2,x3,x4)。RK4的数学形式是支持向量化的，因为按照上述方法处理后，各个变量的相互关系能够通过向量组缩影体现出来。在求解器中，求解过程完全一致，但与求解单方程不同的是，我们需要指定方程组在的方程个数N，并通过该数构建用于储存结果的Nx1矩阵。我们把主要的处理放在了对于函数组表达的预处理上。import numpy as np class RK4_for_equations: def __init__(self, F,n,min,max,totstep,begin): self.F = F self.n = n #指定方程组的个数 self.min = min self.max = max self.step = (self.max-self.min)/totstep self.begin = begin self.totstep = totstep def slover(self): f = self.F x = 0 y = self.begin yn = np.zeros((self.n,1)) yn[:,0] = y step = self.step flag = 1 for i in range(1,self.totstep): K1 = f(x,y) #print("K1:",K1) K2 = f(x+0.5*step,y+0.5*step*K1) #print("K2:",K2) K3 = f(x+0.5*step,y+0.5*step*K2) #print("K3:",K3) K4 = f(x+step,y+step*K3) #print("K4:",K4) y = y + (step/6)*(K1+2*K2+2*K3+K4) yn = np.append(yn,y.reshape(self.n,1),axis=1) flag = flag + 1 x = i * step return yn对于传入函数的预处理简单的来说，我们需要把传入函数写为F(x,y)=[f1(x,y,...,fn(x,y)]格式。以下通过一个实际案例进行描述：已知蒸汽发生器水位控制系统的微分方程如下： x ̇ 1 = G 1 ( q e − q v ) x ̇ 2 = − 1 τ 0 x 2 + G 2 τ 0 q v x ̇ 3 = 1 τ 2 x 4 − G 3 β τ q e x ̇ 4 = − x 3 − 2 ξ τ x 4 + ( 2 ξ β − 1 ) G 3 q e y ( t ) = x 1 + x 2 + x 3 x 1 ( 0 ) = x 2 ( 0 ) = x 3 ( 0 ) = x …