Код проблемы конечного значения: реверсирование векторов внутри функции - PullRequest
Новая
       95

Код проблемы конечного значения: реверсирование векторов внутри функции

0 голосов
/ 14 апреля 2020

В настоящее время я пишу кусок кода для решения системы дифференциальных уравнений, две проблемы начальных значений, а последняя проблема конечных значений. Естественно, для последнего я просто перевернул уравнение и решил его как 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)

Здесь я перевернул векторы sol, sol2 и столбец результирующей матрицы, по сути сделав Это должно быть на единицу больше, чем окончательное значение в моем наборе индексов x, а затем считывать векторы в обратном направлении, используя, например, sol [tt-times, 2].

Однако, если я сделаю следующее вместо использования функции 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)
...