CODEWITHSUNDEEP
rgraphparametersdifferential-equations
P3peM4th.
P3peM4th.

Reputation: 1

Graphing a solution of a system of differential equation varying a parameter respect of time in R

i have been doing this code to plot a solution of Kuramoto's model, my objective was to plot a solution varying the coupling coefficient K and then put all graphs together to see how the coupling evolutions and my attempt was here:

x11()
plot(times1, s1,type="l",xlim=c(0,180),col="white")
value <- seq(0, 3, by=0.1)
pe=1
for(numero in value){
  parameters <- c(N = 2, w1=0.4838171, w2=0.8842176, K=numero)
  state <- c(theta1=0, theta2=0)

  FN <- function (t, state, parameters) {
    with(as.list(c(state,parameters)),{
      dtheta1 <- 2*pi*w1 + (K/N)*(sin(theta1-theta1)+sin(theta2-theta1))
      dtheta2 <- 2*pi*w2 + (K/N)*(sin(theta1-theta2)+sin(theta2-theta2))


      list(c(dtheta1,dtheta2))
    })
  }
  times1 <- seq(1800*(value[pe]*10),1800*(value[1+pe]*10), by=1)*T
  if (pe<=29){
    pe=pe+1
  }
  library(deSolve)
  out <- ode(y = state,times = times1,func = FN,parms=parameters)
  head(out)


  s1 = (1/length(state))*(cos(out[,2])+cos(out[,3]))
  lines(times1, s1,type="l")
}

But now i want to have a K that depends on time and it changes, i.e., take K like some function, example K=K(t), and graphs a solution that at the same time that the time vector varies the K parameter changes and graph that changes with the solution... but i dont know how to solve a system with a parameter varying depending on time, i was reading some questions here, but none was helpful, can you give me some help or how to make a parameter depend on this?, thank you very much!

Upvotes: 0

Views: 90

Answers (1)

tpetzoldt
tpetzoldt

Reputation: 5813

The general approach is to define a "forcing function" that is interpolated with R's approxfun(). As your code contains a few issues that hinder reproducibility, I show this with another example. It is essentially a logistic growth model with a time dependent carrying capacity parameter K:

library(deSolve)

FN <- function (t, state, parms, K_fun) {
  with(as.list(c(state, parms)), {
    K <- K_fun(t)
    dN <- r * N * (1 - N / K)
    list(c(dN))
  })
}

parms <- c(r = 1)
state <- c(N = 1)

times <- seq(0, 50)
value <- seq(10, 20, length.out = length(times))
K_fun <- approxfun(times, value, rule = 2)
out <- ode(y = state, times = times, func = FN, parms = parms, K_fun = K_fun)
plot(out)

More about this can be found in the ?eventshelp page of package deSolve and in the following tutorial: https://tpetzoldt.github.io/deSolve-forcing/deSolve-forcing.html

Upvotes: 1

Related Questions