（知识点）一阶方程的初值问题进行欧拉方法数值求解

2 局部截断误差和顺序

3 代码实现

常微分方程数值求解（ODE）是利用计算机求解微分方程的近似解。将自变量离散化，即

近似解。具体来说，将离散化区间等分，则

其中 是步长， 是节点。 为了使用数值方法求解 ODE，需要首先对 ODE 进行离散化。以下是以下形式的一阶方程：

使用欧拉方法对初值问题进行数值求解。

1 欧拉法

数值解的离散化方法有欧拉法、后向欧拉法、梯形法和修正欧拉法等，下面对上述方法进行说明。

1.1 欧拉法

欧拉法可以通过多种推导方法得到

然后

省略二阶项，我们得到

对微分方程两边在区间上求积分

这里将小矩形左端的函数值作为矩形的高度。因为，所以

1.2 后向欧拉法

使用数值积分法计算矩形的高度也可以使用矩形的右端值，即

然后

这种方法称为欧拉回归法。 上式的两端都包含在内。 右端可以使用欧拉公式进行初步预测，并使用迭代方法求解。

1.3 梯形法

给定一个区间和一个由函数曲线界定的近似梯形，使用梯形面积公式

然后

同样，上式的两端也都包含在内。 首先用欧拉法初步求解并代入上式右端，进行迭代收敛。

1.4 修正的欧拉方法

梯形法计算复杂，每次迭代都要计算函数值，计算量较大。可以使用下面修正的欧拉法

2 局部截断误差和顺序

上述方法可以统一写成如下形式

其中有一个多元函数，称为增量函数（）。 当包含时，该方法是隐式的，否则是显式的。 其中，欧拉法为显式方法，欧拉退却法、梯形法和修正欧拉法为隐式方法。

定义：如果精确解是，则称为

为局部截断误差，该方法具有一阶精度。 这里假设该步骤之前没有错误，从下一步开始就有错误。

这里使用泰勒公式，并进行泰勒展开。 其中，欧拉公式法的精度为1。

其中，后向欧拉公式法的精度为1。

3 代码实现

##---------------Euler-----------------
Euler <- function(fun,a,b,n,y0){
  # 欧拉公式法
  # fun:一阶微分方程函数
  # a,b:求解区间
  # n:区间分段数
  # y0初始值
  x <- numeric()
  y <- numeric()
  h <- (b-a)/n
  x[1] <- a
  y[1] <- y0
  for(k in 1:n){
    x[k+1] = x[k]+h
    y[k+1] = y[k]+h*fun(x[k],y[k])
  }
  result <- data.frame(x = x,y = y)
}

##---------------Euler_back-----------------
Euler_back <- function(fun,a,b,n,y0){
  # 后退欧拉公式法
  # fun:一阶微分方程函数
  # a,b:求解区间
  # n:区间分段数
  # y0初始值
  x <- numeric()
  y <- numeric()
  h <- (b-a)/n
  x[1] <- a
  y[1] <- y0
  for(k in 1:n){
    x[k+1] = x[k]+h
    # 迭代法求y[k+1]
    y0 <- y[k]+h*fun(x[k],y[k])
    for(i in 1:5){
      y1 <- y[k]+h*fun(x[k+1],y0)
      if (abs(y1-y0) <= 0.0001){
        break
      }
      y0 <- y1
    }
    y[k+1] <- y1
  }
  result <- data.frame(x = x,y = y)
  return(result)
}

##---------------Euler_trap-----------------
Euler_trap <- function(fun,a,b,n,y0){
  x <- numeric()
  y <- numeric()
  h <- (b-a)/n
  x[1] <- a
  y[1] <- y0
  for(k in 1:n){
    x[k+1] = x[k]+h
    # 迭代法求y[k+1]
    y0 <- y[k]+h*fun(x[k],y[k])
    for(i in 1:5){
      y1 <- y[k]+h/2*(fun(x[k+1],y0)+fun(x[k],y[k]))
      if (abs(y1-y0) <= 0.0001){
        break
      }
      y0 <- y1
    }
    y[k+1] <- y1
  }
  result <- data.frame(x = x,y = y)
  return(result)
}

##---------------Euler_correct-----------------
Euler_correct <- function(fun,a,b,n,y0){
  x <- numeric()
  y <- numeric()
  h <- (b-a)/n
  x[1] <- a
  y[1] <- y0
  for(k in 1:n){
    x[k+1] = x[k]+h
    y0 <- y[k]+h*fun(x[k],y[k])
    y[k+1] <- y[k]+h/2*(fun(x[k+1],y0)+fun(x[k],y[k]))
  }
  result <- data.frame(x = x,y = y)
  return(result)
}

myfun <- function(x,y){
  sin(5*y)+cos(2*x)
}
a = -10;b = 10;n = 80;y0 = 1
res1 <- Euler(myfun,a,b,n,y0)
plot(res1$x,res1$y,type = "l",lwd = 2,xlab = "x",ylab = "y",col = "black",
     cex.lab = 1.5,cex.axis = 1.5)

res2 <- Euler_back(myfun,a,b,n,y0)
lines(res2$x,res2$y,type = "l",lwd = 2,col = "red")

res3 <- Euler_trap(myfun,a,b,n,y0)
lines(res3$x,res3$y,type = "l",lwd = 2,col = "blue")

res4 <- Euler_correct(myfun,a,b,n,y0)
lines(res4$x,res4$y,type = "l",lwd = 2,col = "green")

grid(lwd = 1,col = "black")

legend("top",legend = c("Euler","Euler_back","Euler_trap","Euler_correct"),
       col = c("black","red","blue","green"),lty = rep(1,4),lwd = rep(2,4),
       bty  = "n", ncol = 2,text.width = 0.8,cex = 1,y.intersp = 0.8)