r - 终值问题代码：反转函数内部的向量

Question

我目前正在编写一大段代码来解决微分方程组，其中两个是初值问题，最后一个是终值问题。当然，对于后者，我只是颠倒了方程式并将其作为 IVP 求解。（顺便说一句，我正在使用 deSolve）。然而，我发现有多种方法可以为我的问题反转向量，只有一种方法似乎给出了一个合理的答案，其他方法都疯狂地爆炸了，我不知道为什么。

这是有效代码的一个片段（下面的完整代码用于上下文）：

linearzreverse<-function(times,y,parms){
    tt=tail(x,n=1)+1
    dw.dx<-(-(sol[tt-times,2]/(a*dt)+(p[1]-p[2]+a)/a))*y[1]-(sol2[tt-times,2]-resultsmatrix[tt-times,val])/(a*dt)
    return(list(dw.dx))
  }
  y2<-(-tail(sol2[,2],n=1)/tail(sol[,2],n=1))
  sol3<-ode(y=y2,times=x,func=linearzreverse,parms=p)

在这里，我通过本质上使 tt 比我的索引集 x 中的最终值多一个，然后使用 sol[tt-times,2] 向后读取向量来反转向量 sol、sol2 和一列 resultsmatrix .

但是，如果我使用 rev() 函数执行以下操作

linearzreverse<-function(times,y,parms){
    dw.dx<-(-(rev(sol[,2])[times]/(a*dt)+(p[1]-p[2]+a)/a))*y[1]-(rev(sol2[,2])[times]-rev(resultsmatrix[,val])[times])/(a*dt)
    return(list(dw.dx))
  }
  y2<-(-tail(sol2[,2],n=1)/tail(sol[,2],n=1))
  sol3<-ode(y=y2,times=x,func=linearzreverse,parms=p)

我得到了一个爆炸的解决方案。即使只是在将反向向量放入函数之前定义它们也会爆炸。

有谁知道为什么我所做的只是以不同的方式反转向量时结果会发生变化？

非常感谢，对于完整的上下文，我将在此处包含完整的代码：

library(deSolve)
tstar=1
n<-100 #determines time intervals
dt<-tstar/n
t=seq(from=0,to=tstar,by=dt)
L=100
M=100
x=seq(from=1,to=L+M,by=1)
v0<-ifelse(1-exp(x-L)<0,0,1-exp(x-L)) #initial condition
#chosen constants for the problem
a=0.25
r=0.12
D=0.1
#blank matricies to store each successive output 
resultsmatrix=matrix(data=0,nrow=length(x),ncol=length(t))
solmatrix=matrix(data=0,nrow=length(x),ncol=length(t)-1)
sol2matrix=matrix(data=0,nrow=length(x),ncol=length(t)-1)
sol3matrix=matrix(data=0,nrow=length(x),ncol=length(t)-1)
#first column is initial condition
resultsmatrix[,1]=v0

riccati<-function(times,y,parms){
  dU.dx<-1-y[1]*((p[1]-p[2]+a)/a)-(y[1]^2)/(a*dt)
  return(list(dU.dx))
}

#loop for solving the system of ODEs at each time step
for(val in (1:length(t[-1]))){
  p<-c(r,D)
  y0<-0
  sol<-ode(y=y0,times=x,func=riccati,parms=p)
  solmatrix[,val]=sol[,2]
  linear<-function(times,y,parms){
    du.dx<-(sol[times,2]/(a*dt))*(resultsmatrix[times,val]-y[1])#if this break, try adding [] to make sol[[times,2]] and make change at * too
    return(list(du.dx))
  }
  y1<-1
  sol2<-ode(y=y1,times=x,func=linear,parms=NULL)
  sol2matrix[,val]=sol2[,2]
  linearzreverse<-function(times,y,parms){
    tt=tail(x,n=1)+1
    dw.dx<-(-(sol[tt-times,2]/(a*dt)+(p[1]-p[2]+a)/a))*y[1]-(sol2[tt-times,2]-resultsmatrix[tt-times,val])/(a*dt) #*
    return(list(dw.dx))
  }
  y2<-(-tail(sol2[,2],n=1)/tail(sol[,2],n=1))
  sol3<-ode(y=y2,times=x,func=linearzreverse,parms=p)
  sol3matrix[,val]=rev(sol3[,2])
  resultsmatrix[,val+1]=sol[,2]*rev(sol3[,2])+sol2[,2]
}
#Transforming the problem into the form needed for my work
S=50
yvals=x-L
K=S*exp(yvals)
DiscountResults=S*resultsmatrix%*%diag(exp(-D*t))
truncno=101 #101 for L=M=100
persp(x=head(K,n=truncno),y=t,z=head(DiscountResults,n=truncno),
      xlab="K",ylab="t",zlab="V(K,t)",main="Stock Price S=50", theta = 50, phi = 25,
      ticktype="detailed",axes=T)
#plot(head(K,n=truncno),head(S*resultsmatrix[,101],n=truncno),type="l")
#Plots to track the three individual differential equations
persp(x=x,y=t[-1],z=solmatrix,
      xlab="x",ylab="t",zlab="U_n(x,tau)",main="Stock Price S=50", theta = 50, phi = 25,
      ticktype="detailed",axes=T)
persp(x=x,y=t[-1],z=sol2matrix,
      xlab="x",ylab="t",zlab="u_n(x,tau)",main="Stock Price S=50", theta = 50, phi = 25,
      ticktype="detailed",axes=T)
persp(x=x,y=t[-1],z=sol3matrix,
      xlab="x",ylab="t",zlab="z_n(x,tau)",main="Stock Price S=50", theta = 50, phi = 25,
      ticktype="detailed",axes=T)

标签： rreversedifferential-equations