rfunctionplotnumerical-methodsnon-linear

Plotting a nonlinear function in R

I am currently struggling to plot the next function in R.

(F1 / F2^m * (R1 - F1)^((1 - s) / s)) - (m * (((R1 - F1))^(1 / s) + ((R2 - F2))^(1 / s)) / (F1^m + F2^m)) = 0

Where I want to plot F2 as a function of F1. That is, F2 in y-axis and F1 in the x-axis. The values R1,R2 are fixed constants and s,m are also fixed.

As the function is non-linear, I was suggested to use a numeric approximation to find the values of F2 in a domain of F1 previously defined.

I developed the next code defining first the constants, then the function, and the root finding function as:

rm(list=ls())

library(ggplot2)

R1_values <- c(100)  # Example values for R1
R2_values <- c(100)  # Example values for R2
s_values <- c(1)  # Example values for s
m_values <- c(1)  # Example values for m

# Define the original function to solve for F2
F2_function <- function(F1, F2, R1, R2, s, m) {
  (F1 / F2^m * (R1 - F1)^((1 - s) / s)) - (m * (((R1 - F1))^(1 / s) + ((R2 - F2))^(1 / s)) / (F1^m + F2^m))
}


F1_values <- seq(0.01, 100, by = 0.1)
# Solve for F2
F2_values <- sapply(F1_values, function(F1) {
  tryCatch({
    # 
    uniroot(F2_function, c(0.01, 100), F1 = F1, R1 = R1, R2 = R2, s = s, m = m, tol = 1e-8)$root
  }, error = function(e) NA)
})

data <- data.frame(F1 = F1_values, F2 = F2_values)

ggplot(data, aes(x = F1, y = F2)) +
  geom_line() +
  labs(title = "F2 as a Function of F1",
       x = "F1",
       y = "F2") +
  theme_minimal()

Everything however, returns a NA or extremly large values (randomly knowing when or how). Even when the function is clearly defined in other software (Geogebra) picture attached.

enter image description here

Hence, I would like to ask help of how to do the plot (and recover the values of F2 for each part of the domain of F1).

The equation in mathematical terms looks like:

enter image description here

PD: I think F2 has two results for each F1, but I am also unable to sort this out.

Solution

  • We can solve this explicitly using Ryacas0. Note that there are two roots and that the argument of sqrt is negative if F1 > 50 so using 50 as the upper bound instead of 100 and plotting the two roots (positive root is black) we have:

    library(Ryacas0)

m <- s <- 1
R1 <- R2 <- 100
F1 <- Sym("F1")
F2 <- Sym("F2")
z <- (F1 / F2^m * (R1 - F1)^((1 - s) / s)) -
  (m * (((R1 - F1))^(1 / s) + ((R2 - F2))^(1 / s)) / (F1^m + F2^m))
zs <- Simplify(z)
F2solve <- Solve(zs, F2)
F2solve
## Yacas vector:
## [1] F2 == (200 - 2 * F1 + sqrt((2 * F1 - 200)^2 - 4 * F1^2))/2
## [2] F2 == (200 - 2 * F1 - sqrt((2 * F1 - 200)^2 - 4 * F1^2))/2

f1 <- seq(0.1, 50, .1)
f2p <- (200 - 2 * f1 + sqrt((2 * f1 - 200)^2 - 4 * f1^2))/2
f2m <- (200 - 2 * f1 - sqrt((2 * f1 - 200)^2 - 4 * f1^2))/2

matplot(f1, cbind(f2p, f2m))

    screenshot