Solve a Non-Linear ODE numerically with no initial condition

224 Views Asked by At

I need to solve the following differential equation numerically: $$ C \frac{dT}{dt}=(1-\alpha (T))Q-\epsilon \sigma T^4 $$

Where $Q,\sigma,\epsilon$ are constants and $\alpha(t)=0.5-(0.2)\text{tanh}(\frac{T-265}{10})$

So far I tried a few built in MATLAB like dsolve functions but none of them seem to work. This is the governing equation for the global mean surface temperature (Equation 2.10) from Kaeper and Engels "Mathematics and Climate".

1

There are 1 best solutions below

1
On

Here is an implementation in R programming language of the Euler-scheme to numerically solve the given ODE with some arbitrarily chosen values for $Q,\sigma, \epsilon$ (and presumably $C$ is also a constant?). An initial condition is required. (Are you sure the authors do not use one?)

plot of solution

# Parameters:
C <- 5
Q <- 0.4
epsilon <- 0.02
sigma <- 0.01

# alpha function
alpha <- function(y)
{
  0.5-0.2*tanh((y-265)/10)
}
# the function in the ODE describing the rate of change
f <- function(t, y)
{
  (1/C)*(1-alpha(y))*Q-epsilon*sigma*y^4
}
# Euler scheme implementation
ode_solve <- function(f, y0, t0 = 0, tn = 1, n = 1000)
{
  h <- (tn-t0)/n
  y <- matrix(0, n+1)
  # n+1 points, n sub-intervals
  tt <- seq(t0, tn, length.out = n+1)
  y[1] <- y0
  for(i in 2:(n+1))
  {
    y[i] <- y[i-1] + h*f(tt[i], y[i-1])
  }
  return(data.frame(time = tt, state = y))
}
# Solving with initial condition T(0)=1, for 500 years (assuming time is in years)
z <- ode_solve(f = f, y0 = 1, tn = 500, n = 5000)
plot(z, type = "l")

Comment for clarifications, questions or corrections please.