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}
]
]
I think this should get you started. I leave to you the plot of the second data set.