R语言常微分方程数值解海强作业
**
**
一階微分方程:一元以人口增長為例
我們使用R package deSolve ODEs 函數。
A simplified form of the syntax for solving ODEs is: ode(y, times, func, parms, …)
where times holds the times at which output is wanted, y holds the initial conditions, func is the name of the R function that describes the differential equations,andparmscontainsthe parametervalues(oris NULL)
dy/dx.
Ode()的y相當于y(0) times相當與dx,function 相當與ry(1-y/k)
r <- 1
K <- 10
yini <- 2
derivs <- function(t, y, parms) list(r * y * (1-y/K))
library(deSolve)
times <- seq(from = 0, to = 20, by = 0.2)
out <- ode(y = yini, times = times, func = derivs, parms = NULL)
我們換個y(0)=12,輸出out2.
r <- 1
K <- 10
yini <- 12
derivs <- function(t, y, parms) list(r * y * (1-y/K))
library(deSolve)
times <- seq(from = 0, to = 20, by = 0.2)
out2 <- ode(y = yini, times = times, func = derivs, parms = NULL)
畫出out1,out2.
plot(out, out2, main = “logistic growth”, lwd = 2)
紅線和黑線是y=f(t)的圖像
二.一階多元線性微分方程組,以The Lorenz Model 為例
a <- -8/3
b <- -10
c <- 28
yini <- c(X = 1, Y = 1, Z = 1)
Lorenz <- function (t, y, parms)
{with(as.list(y),{
dX <- aX+YZ
dY <- b*(Y-Z)
dZ<-XY+cY-Z
list(c(dX, dY, dZ))})}
times <- seq(from = 0, to = 100, by = 0.01)
out <- ode(y = yini,times = times, func = Lorenz, parms = NULL)
~~
- 自編函數
~~
歐拉方法 及其改進
Euler<-function
(x,h=0.1,func=f(),inni=c(0,0))
{ y<-inni[2]
for (i in 0:((x-inni[1])/h)) {
x=x+h
y=y+h*f(x,y)}
return (y)
}
adEuler<-function(x,h=0.1,func=f(),inni=c(0,0))
{ y<-inni[2]
for (i in 0:((x-inni[1])/h)) {
p=y+hf(x,y)
x=x+h
c=y+hf(x,p)
y=0.5*(p+c)}
return (y)
}
Runge<-function(x,h=0.1,func=f(),inni=c(0,0))
{ y<-inni[2]
for (i in 0:((x-inni[1])/h)) {
q=hf(x,y)
w=hf(x+0.5h,y+0.5q)
e=hf(x+0.5h,y+0.5w)
r=hf(x+0.5h,y+e)
x=x+h
y=y+(1/6)(q+2w+2e+r)}
return (y)
}
f<-function(x,y){return(y*(1-y**2))}
總結
以上是生活随笔為你收集整理的R语言常微分方程数值解海强作业的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 解释相机中的弥散现象
- 下一篇: canvas实现图片压缩和缩放