0

I want to use deSolve to solve coupled partial differential equations, and am trying to get used to the package by solving easy differential equations with ode. I know the solution to y' = sqrt(1-y^2) is sin(t + c), so I tested it with the starting conditions y(0) = 0 and y(0) = 1, expecting to get sin(t) and sin(t + pi/2), respectively. But the plots of the solutions are completely different, and very strange.

yini <- c(y = 0)
derivs <- function(t, y)
  list(sqrt(1 - y^2))
times <- seq(from = 0, to = 4, by = 0.2)
out <- ode(y = yini, times = times, func = derivs)
head(out, n =3)

yini <- c(y = 1)
out2 <- ode(y = yini, times = times, func = derivs)
plot(out, out2, main = "", lwd = 2)

enter image description here

Can anyone tell me why?

1 Answers1

0

This is as expected, the right side sqrt(1-y^2) is always positive on the domain with value zero at y=1, so solutions are increasing to the constant solution y=1. As the right side is not Lipschitz there, numerical solvers might get a problem when reaching close to that value and interrupt the integration prematurely.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51