Skip to main content
2 of 4
Change the plot range to be consistent with OP's question
Jack LaVigne
  • 14.5k
  • 2
  • 26
  • 37

Introduction

The first thing that I noticed is that the functional form for all of your CDD and LHS functions are identical.

I propose that rather than defining six functions for each, we define one function and change the input arguments to cover the various cases.

cdd[s_, a_, b_, c_] := 10^(a + b*c - b*s)*Exp[-10^(s - c)]

lhs[s_?NumericQ, a_?NumericQ, b_?NumericQ, c_?NumericQ] := 
 NIntegrate[Log[10]*10^u*cdd[u, a, b, c], {u, s, Infinity}]

Note that we use NumericQ for all of the input arguments to lhs. This is important when we plot the results.

As an example, here is the output for the average parameters you had set for the first case

With[
 {
  a1 = -23.29,
  b1 = 0.88,
  c1 = 20.82,
  s = 20.3
  },
 {cdd[s, a1, b1, c1], lhs[s, a1, b1, c1]}
 ]
(* {1.08754*10^-23, 0.00291862} *)

Interval

The problem you were experiencing is that you wanted to use Interval but it was not accepted by NIntegrate so one was stuck.

The strategy will be to create all permutations of the average inputs parameters +/- their deltas and compute the minimum and maximum.

So for the first step, imagine that we have three parameters a, b and c and when their respective delta is subtracted we have the minus version, am, bm and cm.

We can get a list of all permutation of minus

Flatten[Outer[List, {a, am}, {b, bm}, {c, cm}], 2]
(* {{a, b, c}, {a, b, cm}, {a, bm, c}, {a, bm, cm}, {am, b, 
  c}, {am, b, cm}, {am, bm, c}, {am, bm, cm}} *)

Clearly we can perform the same operation for adding the delta.

Next we will define two functions that use lhs but compute the minimum of subtracting the delta values and the maximum of adding the delta values.

The case is covered where delta is 0 for a particular parameter.

lhsMinus[s_?NumericQ, a_?NumericQ, da_?NumericQ, b_?NumericQ, 
  db_?NumericQ, c_?NumericQ, dc_?NumericQ] := Module[
  {
   aList,
   bList,
   cList
   },
  aList = If[da == 0, {a}, {a, a - da}];
  bList = If[db == 0, {b}, {b, b - db}];
  cList = If[dc == 0, {c}, {c, c - dc}];
  
  Min[Map[lhs[s, Sequence @@ #] &, 
    Flatten[Outer[List, aList, bList, cList], 2]]]
  ]

lhsPlus[s_?NumericQ, a_?NumericQ, da_?NumericQ, b_?NumericQ, 
  db_?NumericQ, c_?NumericQ, dc_?NumericQ] := Module[
  {
   aList,
   bList,
   cList
   },
  aList = If[da == 0, {a}, {a, a + da}];
  bList = If[db == 0, {b}, {b, b + db}];
  cList = If[dc == 0, {c}, {c, c + dc}];
  
  Max[Map[lhs[s, Sequence @@ #] &, 
    Flatten[Outer[List, aList, bList, cList], 2]]]
  ]

Plot

Now I'll plot the results for the first set of parameters. For the large spread in the s values, 20.3 to 23.1 I had to increase the delta's or the min and max were too close to be visible.

With[
 {
  a1 = -23.29,
  da = 0.08,
  b1 = 0.88,
  db = 0.08,
  c1 = 20.82,
  dc = 0.04
  },
 LogPlot[
  {
   lhsMinus[s, a1, da, b1, db, c1, dc],
   lhs[s, a1, b1, c1],
   lhsPlus[s, a1, da, b1, db, c1, dc]
   },
  {s, 20.3, 23.1},
  {PlotStyle -> {Blue, Black, Red}},
   PlotRange -> {10^-35, 1}
  ]
 ]

Mathematica graphics

I think this should get you started. I leave to you the plot of the second data set.

Jack LaVigne
  • 14.5k
  • 2
  • 26
  • 37