数值解常微分方程(三)

数值解常微分方程(三)

使用R语言数值求解ODE方程的包为deSolve，基本语法格式为

ode(y, times, func, parms, ...)

y为初始条件，time表示时间(自变量)，func为函数，parms表示参数。绝对和相对误差默认为atol = 1e-6,rtol = 1e-6

---

1 一元微分方程

使用deSolve包求解如下微分方程
$$y^{\prime }=sin(2x)+cos(y),y(0)=1$$

library(deSolve) 
# 一元微分方程
# 参数和初值
m <- 2; n <- 1; yini <- 1
# 微分方程
derivs <- function(x, y, parms) list(sin(m*x)+cos(n*y))
# 自变量区间离散化
x <- seq(from = -10, to = 5, by = 0.1) 
# ode求解
out <- ode(y = yini, times = x, func = derivs, parms = NULL)

yini <- 2
out2 <- ode(y = yini, times = x, func = derivs, parms = NULL)

yini <- 3
out3 <- ode(y = yini, times = x, func = derivs, parms = NULL)

plot(out,out2,out3, main = " ", lwd = 5,xlab = "x",
     ylab = "y",cex.lab = 1.5,cex.axis = 1.5)
grid(lwd = 1,col = "black")

---

2 多元微分方程组

对于多元一阶微分方程组，也可以使用对应的数值方法求解。例如
$$\begin{array}{rl}{y}_i^{\prime }=f_i(x,y_1,y_2,\dots y_N)& \\ y_i(0)& =y_i^0,i=1,2\dots N\end{array}$$
将其写作向量格式
$$\begin{array}{rl}& {\mathbf{y}}^{\prime }=\mathbf{f}(x,\mathbf{y})\\ & \mathbf{y}(x_0)={\mathbf{y}}_0\end{array}$$
对应的R-K四阶公式为
$$\{\begin{array}{lc}{\mathbf{y}}_{k+1}& ={\mathbf{y}}_k+h(\frac{1}{6}{\mathbf{K}}_1+\frac{2}{6}{\mathbf{K}}_2+\frac{2}{6}{\mathbf{K}}_3+\frac{1}{6}{\mathbf{K}}_4)\\ {\mathbf{K}}_1& =f(x_k,{\mathbf{y}}_k)\\ {\mathbf{K}}_2& =f(x_k+h/2,{\mathbf{y}}_k+h{\mathbf{K}}_1/2)\\ {\mathbf{K}}_3& =f(x_k+h/2,{\mathbf{y}}_k+h{\mathbf{K}}_2/2)\\ {\mathbf{K}}_4& =f(x_k+h,{\mathbf{y}}_k+h{\mathbf{K}}_3)\end{array}$$
使用deSlove包采用rk4方法即可求解。例如求解如下三元一阶常微分方程组
$$\{\begin{array}{l}{y}_1^{\prime }=sin(y_2{y}_3)\\ y_2^{\prime }=cos(y_1{y}_3)\\ y_3^{\prime }=sin(y_1-y_2)\\ (y_1,y_2,y_3)=(0,0,0),x\in (0,20)\end{array}$$

# 多元微分方程组
yini <- c(0,0,0)
multieq <- function(x,y,parms){
  dy1 <- sin(y[2] * y[3]+y[2])
  dy2 <- cos(y[1] * y[3])
  dy3 <- y[2]*y[1]
  list(c(dy1, dy2, dy3))
}
t <- seq(from = 0, to = 10, by = 0.05) 
out <- ode (times = t, y = yini, func = multieq, 
            parms = NULL, method = rkMethod("rk45ck"))
matplot(x = out[,1], y = out[,-1], type = "l", lwd = 2, lty = "solid", 
        col = c("red", "blue", "black"), xlab = "time", ylab = "y") 
legend("bottomleft", col = c("red", "blue", "black"), 
       legend = c("y1", "y2", "y3"), lwd = 2)
grid(lwd = 1,col = "black")

参数method可以替换为其他数值方法

# R-K方法
rkMethod()
# [1] "euler"   "rk2"     "rk4"     "rk23"    "rk23bs"  "rk34f"   "rk45f"  
# [8] "rk45ck"  "rk45e"   "rk45dp6" "rk45dp7" "rk78dp"  "rk78f"   "irk3r"  
# [15] "irk5r"   "irk4hh"  "irk6kb"  "irk4l"   "irk6l"   "ode23"   "ode45"

# 经典4阶RK方法
rkMethod("rk4")

将$y_1,y_2,y_3$绘制为3D坐标曲线

library(scatterplot3d) 
library(FactoClass)
scatterplot3d(out[,-1], type = "l", lwd = 2, xlab = "x", ylab = "y",
              zlab = "z",angle = 50,color = "red", grid=FALSE,box = FALSE,
              cex.lab = 1.3,cex.axis = 1.3)

addgrids3d(out[,-1], grid = c("xy", "xz", "yz"),angle = 50,col.grid = "blue",
           lty.grid = 5)

---

3 高阶微分方程

高阶微分方程可以通过换元法使高阶微分方程转化为一阶微分方程组。例如对于一般形式的高阶微分方程
$$y^{(m)}=f(x,y,y^{\prime },y^{\prime \prime }\dots y^{(m-1)})$$
其中初始条件
$$y(x_0)=y_0,y^{\prime }(x_0)=y_0^{\prime },\dots y^{(m-1)}(x_0)=y_0^{(m-1)}$$
引进新变量
$$y_1=y,y_2=y_1^{\prime },\dots ,y_m=y^{(m-1)}$$
得到如下一阶$m$元微分方程组
$$\{\begin{array}{l}{y}_1^{\prime }=y_2\\ y_2^{\prime }=y_3\\ ⋮\\ y_{m-1}^{\prime }=y_m\\ y_m^{\prime }=f(x,y_1,y_2,\dots ,y_m)\end{array}$$
对应初始条件为
$$y_1(x_0)=y_0,y_2(x_0)=y_0^{\prime },\dots ,y_m(x_0)=y_0^{(m-1)}$$
此时使用多元一阶微分方程组的数值解法即可。

---

参考书籍：

李庆扬等，数值分析[M]. 北京: 清华大学出版社, 2008.