plotting two functions in one 2D graph in R

Question

I have the following function (for example):

f<-function(x,a,b) a*sinh(x/b)

where x is real from -L to L. Also, a and b, discretely increments from a1 to a2 and from b1 to b2, respectively. (e.g. for loop with some steps)

I'm wondering how I can plot f_max (maximum of f for -L < x < L ) and f_avg (average of f over -L < x < L ) versus a both in the same plot for when b is also changing. I'm looking for something like this:

If anything is missing in my question, please make any assumption as you wish (without asking). It should help anyways as my question is general.

As you may notice I am trying to visualize 2 different 3D graphs (f_max vs a vs b) and (f_avg vs a vs b) in one 2D graph. ggplot2 plus coloring and linear curve fitting is highly appreciated.

Answer 1

As the comments indicated, you need to create a grid of values first, but the optim and integrate suggestions aren't needed for this sort of task.

grid.xab <- expand.grid(a = c(4,5,6), 
                        b = c(1,1.5,2), 
                        x = seq(-2, 3, length.out = 50))

    a   b  x
  1 4 1.0 -2
  2 5 1.0 -2
  3 6 1.0 -2
  4 4 1.5 -2
  5 5 1.5 -2
  6 6 1.5 -2
         ...
445 4 1.5  3
446 5 1.5  3
447 6 1.5  3
448 4 2.0  3
449 5 2.0  3
450 6 2.0  3

You see this makes every possible combination of the a and b values with the range of x values. Depending on the precision you need, you can change length.out =. Now you evaluate f(x) at every point.

---

library(dplyr)
library(ggplot2)

grid.xabf <- grid.xab %>% mutate(fx = a*sinh(x/b)) 

  a   b  x         fx
  1 4 1.0 -2 -14.507442
  2 5 1.0 -2 -18.134302
  3 6 1.0 -2 -21.761162
  4 4 1.5 -2  -7.060142
  5 5 1.5 -2  -8.825177
  6 6 1.5 -2 -10.590212
                    ...
445 4 1.5  3  14.507442
446 5 1.5  3  18.134302
447 6 1.5  3  21.761162
448 4 2.0  3   8.517118
449 5 2.0  3  10.646397
450 6 2.0  3  12.775677

---

You can visualize this:

ggplot(grid.xabf, aes(x, fx)) + 
  geom_line(aes(color = factor(a), linetype = factor(b)))

---

Now we just need to collapse all the x values to get the maximum and average values of f(x) for each a and b.

f.hat <- grid.xabf %>% 
  group_by(a, b) %>% 
  summarise(fmax = max(fx), favg = mean(fx))

# A tibble: 9 x 4
# Groups:   a [?]
      a     b      fmax     favg
  <dbl> <dbl>     <dbl>    <dbl>
1     4   1.0 40.071500 5.203415
2     4   1.5 14.507442 2.113929
3     4   2.0  8.517118 1.307466
4     5   1.0 50.089375 6.504269
5     5   1.5 18.134302 2.642412
6     5   2.0 10.646397 1.634333
7     6   1.0 60.107250 7.805122
8     6   1.5 21.761162 3.170894
9     6   2.0 12.775677 1.961199

---

And now we can finally make the figure you want:

p1 <- ggplot(f.hat) + 
  geom_line(aes(a, fmax, color = factor(b), linetype = "f(x) max")) +
  geom_line(aes(a, favg, color = factor(b), linetype = "f(x) avg")) +
  geom_point(aes(a, fmax, color = factor(b)), shape = 16) +
  geom_point(aes(a, favg, color = factor(b)), shape = 1)

p1

But we can also make it a bit nicer:

p1 + 
scale_linetype_manual(values = c("f(x) max" = "solid",
                                   "f(x) avg" = "dashed"),
                        labels = c("f(x) max" = expression(bold(max)~italic(f)(x)), 
                                   "f(x) avg" = expression(bar(italic(f)(x)))),
                        name = expression(widehat(italic(f)(x)))) +
  labs(y = expression(widehat(italic(f)[-2]^3*(x))),
       color = "b",
       title = expression(italic(f)(x)==a%.%bold(sinh)*bgroup("(",frac(x, b),")")))