mmcdara

4924 Reputation

17 Badges

7 years, 6 days

MaplePrimes Activity


These are replies submitted by mmcdara

@NanaP 

There are always questions to which several players contribute and others that remain more "confidential".

You'll quickly understand all this with a little practice.

The most upsetting thing (even for me) is that you'll never know why your question remains unanswered (is it that no one has an answer to offer, or is it that the problem is so uninteresting or so poorly presented that no one cares to answer it?)

So a good practice is probably to present your problem in the most consise, but nevertheless detailed, form and to deliver what you have been able to achieve with Maple (generally people don't like an image or a reference, or a piece of some article to figure out what the OP [the acronym used to designate the person who asked the question] wants).

Concerning successive questions, I agree that when no one have given you a reply, it's very tempting to bring up a question again. But you have to account for jet lag and half a day may separate a question and its answer.

Concerning the difficulty to find Maple solutions or algorithms in the web, I agree that it is extremely difficult.
Theitems you get refer generally to Mapleprimes. An advice: even if the Search/Advance Search feature could be improved, do not hesitate to browse amonq the Questions and Post to seek if similar problems have not been alerady answeerd.
Another interesting link is https://www.maplesoft.com/applications/ 

Hope to see you soon.

To conclude, I believe, this thread, here is a more complete version of BDS_2.mw where I added several estimators of the left and right endpoints of the underlying distribution, given a sample from it and the lnowledge that it is left and right bounded.
Tou will see that the estimators are quite good is the distribution is "smooth" and quite rough if it is sharp (here left peaked).
BSD_2b.mw 

If you'd like to dig up further on the subject, here are some companions articles (their understanding resquires some knowledge in Extreme Values Theory).​​
https://boris.unibe.ch/36886/1/10687_2007_Article_52.pdf
https://arxiv.org/pdf/1412.3972.pdf
https://projecteuclid.org/journals/annals-of-statistics/volume-12/issue-4/Estimating-an-Endpoint-of-a-Distribution-with-Resampling-Methods/10.1214/aos/1176346811.full

https://www.d.umn.edu/~yqi/downloads/19sinica.pdf
https://mendel-journal.org/index.php/mendel/article/view/28/30

@NanaP 

UPDATED
(simpler definition of the SubjectiveBeta Distribution)

COMMENT 1
You should stop posting replies while the previous one doesn't receive an answer.

COMMENT 2:
Note from @NanaP: The Use_formulas.mw file by @mmcdara stopped at the 25 percentile "theoretical" calculation and did not go beyond this point until a histogram of simulated values is created

You're partially wrong: the first file I sent you (BDS.mw) "goes as far as the histogram", maybe you didn't pay attention?

COMMENT 3: (in response to your next reply)
My file Use_formulas.mw contains indeed some errors: I used the moments to compute the Skewness and Kurtosis instead of the central moments.
Apologies.

Here is a file which "[achieve] the ultimate goal of being able to simulate from a custom defined distribution of SB(minimum, mode, mean, maximum) .",  although I didn't particularly like the way you worded it.

The procedure SubjectiveBeta builds a "Subjective Beta Distribution" given its mean, mode, min value and max value

SubjectiveBeta := proc(me, mo, P, Q)
   # me = mean
   # mo = mode
   # P  = xmin
   # Q  = xmax

restart:

with(Statistics):

# Basically you want to parameterize  a Betadistribution by its mean and its mode
# (the smin and xmax stuff is just a translation and a scaling)
#
# S let's derive the expression of she usual parameters alpha and beta from the
# mode and the mean.

U  := RandomVariable('Beta' (alpha, beta)):
ME := Mean(U);
MO := Mode(U);

solve({ME=MyMean, MO=MyMode}, {alpha, beta}); # no solution found
 

ME := alpha/(alpha+beta)

 

piecewise(alpha < 1 and beta < 1, {0, 1}, 1 < alpha and 1 < beta, (-1+alpha)/(alpha+beta-2), alpha < 1 and 1 <= beta, 0, 1 <= alpha and beta < 1, 1, FAIL)

(1)

# No solution being found, lets assume alpha > 1, beta > 1 and let's try to find one

MO := Mode(U) assuming alpha > 1, beta > 1;

reparam := solve({ME=MyMean, MO=MyMode}, {alpha, beta});

(-1+alpha)/(alpha+beta-2)

 

{alpha = -MyMean*(2*MyMode-1)/(MyMean-MyMode), beta = (MyMean-1)*(2*MyMode-1)/(MyMean-MyMode)}

(2)

# Here we got a solution provided MuMean <> MyMode.
#
# The developpements which follows maje sense only if alpha > 1, beta > 1 and
# MyMean <>MyMode.
#
# I already pointed out the absurdity of such a re-parameterization but it is
# your responsibility to do it and not of my concern.
 

# Preliminary work to get the desired formulas

U := RandomVariable('Beta'(alpha, beta)):
V := x__min + (x__max-x__min)*U:

eval(reparam, {MyMean=me, MyMode=mo}):
eval(reparam, {MyMean=(me-P)/(Q-P), MyMode=(mo-P)/(Q-P)}):
simplify(%);

{alpha = (-2*mo+P+Q)*(-me+P)/((me-mo)*(-Q+P)), beta = -(-2*mo+P+Q)*(-me+Q)/((me-mo)*(-Q+P))}

(3)

SubjectiveBeta := proc(me, mo, P, Q)
   # me = mean
   # mo = mode
   # P  = xmin
   # Q  = xmax

   local A  := (-2*mo+P+Q)*(-me+P)/((me-mo)*(-Q+P));
   local B  := -(-2*mo+P+Q)*(-me+Q)/((me-mo)*(-Q+P));
   local f0 := unapply(PDF('Beta'(A, B), t), t);
   local F0 := unapply(CDF('Beta'(A, B), t), t);

   printf("According to Maple help pages : %a\n", [nu=A, omega=B]);

   Distribution(
     PDF = unapply(f0((x-P)/(Q-P)) / (Q-P), x),

     CDF = unapply(F0((x-P)/(Q-P)), x),

     Mean = me,

     Variance = (Q-P)^2*Variance('Beta'(A, B)),

     Skewness = simplify(CentralMoment('Beta'(A, B), 3) / Variance('Beta'(A, B))^(3/2)),

     Kurtosis = simplify(CentralMoment('Beta'(A, B), 4) / Variance('Beta'(A, B))^2),

     Mode = mo,

     Conditions = [Q > P],

     RandomSample = proc(N::nonnegint)
                      P +~ (Q-P) *~ Sample('Beta'(A, B), N)
                    end proc
   )
end proc:

BSD := RandomVariable(SubjectiveBeta(1500, 1400, 1000, 2100))

According to Maple help pages : [nu = 15/11, omega = 18/11]

 

_R8

(4)

N := 10^4:
S := Sample(BSD, N):
f := PDF(BSD, x):

plots:-display(
  Histogram(S),
  plot(f, x=1000..2100, color=red, thickness=2)
);

 

plots:-display(
  ScatterPlot(sort(S), [$1..N]/~N, symbol=point, legend="Empirical CDF"),
  plot(CDF(BSD, x), x=1000..2100, color=red, thickness=2, legend="CDF")
);

 

'Mean'          = evalf(Mean(BSD));
'EmpiricalMean' = Mean(S);
 print():

'Variance'          = evalf(Variance(BSD));
'EmpiricalVariance' = Variance(S);
 print():
'Skewness'          = evalf(Skewness(BSD));
'EmpiricalSkewness' = Skewness(S);
 print():

'Kurtosis'          = evalf(Kurtosis(BSD));
'EmpiricalKurtosis' = Kurtosis(S);
 print():

'Mode'          = evalf(Mode(BSD));
'EmpiricalMode' = Mode(S);

Statistics:-Mean = 1500.

 

EmpiricalMean = HFloat(1496.475069708961)

 

 

Statistics:-Variance = 75000.

 

EmpiricalVariance = HFloat(75095.54978281037)

 

 

Statistics:-Skewness = .1460593486

 

EmpiricalSkewness = HFloat(0.16168690502581273)

 

 

Statistics:-Kurtosis = 2.026666667

 

EmpiricalKurtosis = HFloat(2.0347480957357007)

 

 

Statistics:-Mode = 1400.

 

EmpiricalMode = HFloat(1568.508189104639)

(5)

BSD := RandomVariable(SubjectiveBeta(1500, 2000, 1000, 2100))

According to Maple help pages : [nu = 9/11, omega = 54/55]

 

_R17

(6)

plot(PDF('Beta'(9/11, 54/55), t), t=0..1);
Mode('Beta'(9/11, 54/55));

 

{0, 1}

(7)

S := Sample(BSD, 10^4):
f := PDF(BSD, x):

plots:-display(
  Histogram(S),
  plot(f, x=1000..2100, color=red, thickness=2)
);


'Mean'          = evalf(Mean(BSD));
'EmpiricalMean' = Mean(S);
 print():

'Variance'          = evalf(Variance(BSD));
'EmpiricalVariance' = Variance(S);
 print():
'Skewness'          = evalf(Skewness(BSD));
'EmpiricalSkewness' = Skewness(S);
 print():

'Kurtosis'          = evalf(Kurtosis(BSD));
'EmpiricalKurtosis' = Kurtosis(S);
 print():

'Mode'          = evalf(Mode(BSD));
'EmpiricalMode' = Mode(S);

 

Statistics:-Mean = 1500.

 

EmpiricalMean = HFloat(1489.5647600531224)

 

 

Statistics:-Variance = 107142.8571

 

EmpiricalVariance = HFloat(107409.17519359972)

 

 

Statistics:-Skewness = .1607921296

 

EmpiricalSkewness = HFloat(0.1944815270127249)

 

 

Statistics:-Kurtosis = 1.780701754

 

EmpiricalKurtosis = HFloat(1.7984991021830183)

 

 

Statistics:-Mode = 2000.

 

EmpiricalMode = HFloat(1005.7178731055807)

(8)

BSD := RandomVariable(SubjectiveBeta(1500, 1500, 1000, 2100));

Error, (in SubjectiveBeta) numeric exception: division by zero

 

 

Download BSD_2b.mw

I hope that's all right with you, because I'm done with this subject.

@NanaP 

I'm not sure I understood correctly so correct me if I'm wrong:

  • You have a random variable X with a SB (short for "subjective Beta") distribution of parameters xmin, xmax, alpha, beta (I'm not sure of this point, maybe they are xmin, xmax, mode, mean ?).
  • You have a sample S of X, let's say of size N, from which you compute some empirical statistics, for instance
    the mean 
    m := (1/N)*add(s[n], n=1..N)
    
  • Some facts are:
    • the above empirical statistics is the best unbiased asymptotic estimator of the Mean(X) (of course you already knew that);
    • the classical estimation of the variance 
      v := (1/(N-1))*add((s[n]-m)^2, n=1..N)

      also has this same property.

    • The case of Skewness and Kurtosis is a little bit more complex. If I'm not mistaken there are no formulas for unbiased estimator (nearly unbiased estimators do exist).

    • The case of xmin and xmax is even more complex.
      I have some material concerning what is called "the endpoint distribution", that is the estimation of xmin or xmax from sample S.
      Here is a first reference, https://arxiv.org/pdf/1412.3972.pdf, do not hesitate to ask for some more.

    • Concerning the estimation of the mode there is no unbiased estimator, but a lot of estimators do exist.

@NanaP 

Thank you for introducing me to something I didn't know about, even though I've been working in a parallel domain for a long time.

Hope my contribution will be useful

@acer 
I had observed on my own that isolate didn't work for this example.
So you are right saying that the imaginary unit is not the cause of my difficulty in the sense that I could faced them without having I.

It's just that for the particular case b=1 isolate works if c <> I, and not instead, so my quick jump to an erroneous conclusion.

@acer 

Thank you acer for all these informations.

@NanaP 

The case 

Min := 1000; Mlikely := 1400; Avg := 1500; Max := 2100;

can be solved by using the definitions of Mean, Variance, Skewness, Kurtosis, Mode and Quantile:

The Beta Subjective Distribution

restart

 

with(Statistics):

Min := 1000; Mlikely := 1400; Avg := 1500; Max := 2100

2100

(1)

Mid := (Min+Max)*(1/2)

1550

(2)

alpha := 2*(Avg-Min)*(Mid-Mlikely)/((Avg-Mlikely)*(Max-Min))

15/11

(3)

beta := alpha*(Max-Avg)/(Avg-Min)

18/11

(4)

 

f := simplify(piecewise(Min <= x and x <= Max, (x-Min)^(alpha-1)*(Max-x)^(beta-1)/(Beta(alpha, beta)*(Max-Min)^(alpha+beta-1)), 0))

piecewise(x < 1000, 0, x <= 2100, (1/1210000)*(x-1000)^(4/11)*(2100-x)^(7/11)/Beta(15/11, 18/11), 2100 < x, 0)

(5)

MD := Distribution(PDF = unapply(f, x), Conditions = [`and`(Min < Mlikely and Mlikely < Max and Min < Avg and Avg < Max and Mlikely <> Avg, piecewise(Avg < Mlikely, Mlikely > Mid, Mid > Mlikely))])

(6)

X := RandomVariable(MD)

_R

(7)

plot(f, x=Min..Max, gridlines=true)

 

mom := r -> evalf(Int(x^r*f, x=Min..Max));

_Mean     := mom(1);
_Variance := mom(2)-_Mean^2;
_Skewness := mom(3) / mom(2)^(3/2);
_Kurtosis := mom(4) / mom(2)^2:

Optimization:-Maximize(f, x=Min..Max);
_Mode := eval(x, %[2]);

Q := z -> evalf(Int(f, x=Min..z)):
p := 0.25;
_Q__||p := fsolve('Q(z)'=p, z=`if`(p < 0.5, Min, Max));
 

proc (r) options operator, arrow; evalf(Int(x^r*f, x = Min .. Max)) end proc

 

1500.000000

 

75000.000

 

1.048051997

 

[HFloat(0.001181159064529679), [x = HFloat(1399.9999998692258)]]

 

HFloat(1399.9999998692258)

 

.25

 

1274.300220

(8)

 

Download Use_formulas.mw

Note that the mode is not necesarily  unique (Beta(a, a) has two modes at x=0 and x=1 when a < 1 and a continuum of modes whan a=1).
So, I don't know what you want to do, but this re-parameterization you want to use is extremely dangerous.

NULL

The Beta Subjective Distribution

restart

 with(Statistics):

Min := 0; 1; Mlikely := 1/2; 1; Avg := 1/2; 1; Max := 1

1

(1)

Mid := (Min+Max)*(1/2)

1/2

(2)

alpha := 2*(Avg-Min)*(Mid-Mlikely)/((Avg-Mlikely)*(Max-Min))

Error, numeric exception: division by zero

 

beta := alpha*(Max-Avg)/(Avg-Min)

alpha

(3)

 

Download Distribution-Beta-Subjective_Failure.mw

About the Beta distribution there are only, up to my knowledge, these parameterizations:

  • shape (a) and scale (b), the most common ant the one Maple adopts;
  • mean (M) and variance (V);
  • mean and precision (inverse of variance, noted tau below);

As you can see there is a 1-to-1 map between each of pair of parameterization

m := Mean(Beta(a, b)):
v := Variance(Beta(a, b)):
solve({m=M, v=V}, {a, b});

solve({m=M, 1/v=tau}, {a, b});
        /        / 2        \      / 2        \        \ 
        |      M \M  - M + V/      \M  - M + V/ (M - 1)| 
       < a = - --------------, b = -------------------- >
        |            V                      V          | 
        \                                              / 
  /      3        2              / 2                \        \ 
 { a = -M  tau + M  tau - M, b = \M  tau - M tau + 1/ (M - 1) }
  \                                                          / 

This VOSE article you refer about seems to be oriented to risk assesment, for which this Subjective Beta Distribition might be, according to the authors, and likely under some restrictions an interesting tool to achieve this goal.
I'm working on the domain of Reliability Analysis (assessment of probability of failures due to extreme events in engineering, something probably not that far from the "risk assessment" VOSE talk about?]) for more than 25 years and I had never heard of the approach developped by VOSE (but each domain develops it's own method with little inter-operability).

Some intersting links:

Probability ontology
http://www.math.wm.edu/~leemis/2008amstat.pdf
https://en.wikipedia.org/wiki/Relationships_among_probability_distributions


Different distributions
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5013898/bin/supp_btw170_Swat_et_al_Supplementary1_ProbOnto_KBprintout.pdf

Different parameterizations
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5013898/bin/supp_btw170_Swat_et_al_Supplementary2_ProbOnto_useInPharmML.pdf

@NanaP 

Does it mean you want to re-parameterize a Generalized Beta distribution by min, max, mode, min instead of alpha, beta, min and max ?

If it so you will have a problem because the relation (alpha, beta) --> (mode, mean) is not a 1-to-1 mapping.
Take the simplest case alpha > 1, beta=alpha, xmin=0, xmax=1 (thus the Generalized Beta is just Beta(alpha, alpha)).
As Mean(Beta(alpha, alpha)) = Mode(Beta(alpha, alpha)) = 1/2 whatever alpha > 1, a re-parametrization of the Beta distribution in term of mode and mean doesn't say anything about the pdf, cdf, variance, ... of this distribution.

@C_R 

but the expression 

I*Int(f(x), x) = something - 2*I*Int(f(x), x);

comes from some computation and I didn't want to extract the Int component from it's left hand side and isolate it from the equality.
In fact I wrote f(x) but I don't even know a priori what f(x) is.
So you can do something like this:

expr := I*Int(f(x), x) = something - 2*I*Int(f(x), x):
J := select(has, [op(lhs(expr))], Int)[]:
isolate(expr, J)
                    /                         
                   |             1            
                   |  f(x) dx = -- I something
                   |             3            
                  /                           


The problem is that you replace I by some constant c, then isolate works as expected:
 

c*Int(f(x), x) = something - 2*c*Int(f(x),x):
isolate(%, lhs(%))
                   /  /        \              
                   | |         |   1          
                 c | |  f(x) dx| = - something
                   | |         |   3          
                   \/          /              
I*Int(f(x), x) = something - 2*I*Int(f(x),x):
isolate(%, lhs(%))
          /  /        \                   /  /        \
          | |         |                   | |         |
        I | |  f(x) dx| = something - 2 I | |  f(x) dx|
          | |         |                   | |         |
          \/          /                   \/          /

 

@sursumCorda 

Thank you for your explanation.

https://www.mapleprimes.com/questions/236509-How-Do-I-Solve-This-One.

As I'm in a good mood today, here's a solution

# A maple expression

PutMeInTheBox = `#mrow(msubsup(mo("&Sigma;"),mrow(mi("k"),mo("=120")),mo("1300")),mo("i",mathcolor="white"),msup(mrow(mo("&lpar;"),mo("2"),mo("i",mathcolor="white"),mi("x-k"),mo("&rpar;")),mi("k")))`

PutMeInTheBox = `#mrow(msubsup(mo("&Sigma;"),mrow(mi("k"),mo("=120")),mo("1300")),mo("i",mathcolor="white"),msup(mrow(mo("&lpar;"),mo("2"),mo("i",mathcolor="white"),mi("x-k"),mo("&rpar;")),mi("k")))`

(1)

 

NULL

Download Maple_expression.mw

Can't you provide a Maple code for us to see what you are capable of and then guide you towards the solution?
The formulation of your question sounds like "do my homework, I don't want to bother with it".

@KIRAN SAJJAN 

Could you explain more precisely what you want because your last worksheet doesn't contain any error and your request is not clear (at least for me).
I imagine M is the Mach number, Cf is the friction coefficient and Nu is the Nusselt number, aren't they?
Do you mean that you want to plot , for instance, aNusselt number defined by 

(-a4-Rd)*diff(T(Y), Y)

for different values of the Mach number?

Could it be this:
 

NULL

restart; with(plots)

PDEtools[declare]((F, T, G, H)(Y), prime = Y):

F(Y)*`will now be displayed as`*F

 

T(Y)*`will now be displayed as`*T

 

G(Y)*`will now be displayed as`*G

 

H(Y)*`will now be displayed as`*H

 

`derivatives with respect to`*Y*`of functions of one variable will now be displayed with '`

(1)

p1 := 0.1e-1:

rf := 1050:

sigma1 := 25000:

sigma2 := 0.210e-5:

sigma3 := 6.30*10^7:

sigma4 := 10^(-10):

sigma5 := 1.69*10^7:

sigma6 := 4.10*10^7:

``

S1 := .5:

alp := 2:

``

``

B1 := 1+2.5*p+6.2*p^2:

a1 := B1*p1+B2*p2+B3*p3:

a2 := 1-p1-p2-p3+p1*rs1/rf+p2*rs2/rf+p3*rs3/rf:

a3 := 1-p1-p2-p3+p1*rs1*cps1/(rf*cpf)+p2*rs2*cps2/(rf*cpf)+p3*rs3*cps3/(rf*cpf):

a4 := B4*p1+B5*p2+B6*p3:

``

a5 := 1+3*((p1*sigma1+p2*sigma2+p3*sigma3)/sigmaf-p1-p2-p3)/(2+(p1*sigma1+p2*sigma2+p3*sigma3)/((p1+p2+p3)*sigmaf)-((p1*sigma1+p2*sigma2+p3*sigma3)/sigmaf-p1-p2-p3)):

``

a6 := B1*p1+B2*p2+B3*p3:

a7 := 1-p1-p2-p3+p1*rs4/rf+p2*rs5/rf+p3*rs6/rf:

a8 := 1-p1-p2-p3+p1*rs4*cps4/(rf*cpf)+p2*rs5*cps5/(rf*cpf)+p3*rs6*cps6/(rf*cpf):

a9 := B7*p1+B8*p2+B9*p3:

``

a10 := 1+3*((p1*sigma4+p2*sigma5+p3*sigma6)/sigmaf-p1-p2-p3)/(2+(p1*sigma4+p2*sigma5+p3*sigma6)/((p1+p2+p3)*sigmaf)-((p1*sigma4+p2*sigma5+p3*sigma6)/sigmaf-p1-p2-p3)):

W := sum(b[i]*Y^i, i = 0 .. 7);

Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0]

 

Y^7*c[7]+Y^6*c[6]+Y^5*c[5]+Y^4*c[4]+Y^3*c[3]+Y^2*c[2]+Y*c[1]+c[0]

 

Y^4*d[4]+Y^3*d[3]+Y^2*d[2]+Y*d[1]+d[0]

 

Y^4*h[4]+Y^3*h[3]+Y^2*h[2]+Y*h[1]+h[0]

(2)

F := a1*(1+1/bet)*(diff(W, `$`(Y, 2)))+a2*Ra*(diff(W, Y))+A-a5*M*W-S2*W^2+a2*Gr*Theta-S*betu*(W-U) = 0;

1-.5*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^2+1.622380952*Y*b[2]-1.296703274*M*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])+64.1780496*Y^5*b[7]+45.8414640*Y^4*b[6]+30.5609760*Y^3*b[5]+18.3365856*Y^2*b[4]+9.1682928*Y*b[3]+5.678333332*Y^6*b[7]+4.867142856*Y^5*b[6]+4.055952380*Y^4*b[5]+3.244761904*Y^3*b[4]+2.433571428*Y^2*b[3]+0.5e-1*d[0]-0.5e-1*b[0]+.8111904760*b[1]+3.0560976*b[2]+.8111904760*c[0]-0.5e-1*b[5]*Y^5-0.5e-1*b[6]*Y^6-0.5e-1*b[7]*Y^7+.8111904760*c[1]*Y+.8111904760*c[2]*Y^2+.8111904760*c[3]*Y^3+.8111904760*c[4]*Y^4+.8111904760*c[5]*Y^5+.8111904760*c[6]*Y^6+.8111904760*c[7]*Y^7+0.5e-1*d[1]*Y+0.5e-1*d[2]*Y^2+0.5e-1*d[3]*Y^3+0.5e-1*d[4]*Y^4-0.5e-1*b[1]*Y-0.5e-1*b[2]*Y^2-0.5e-1*b[3]*Y^3-0.5e-1*b[4]*Y^4 = 0

(3)

T := (a4+Rd)*(diff(Theta, `$`(Y, 2)))+a3*Pr*Ra*(diff(Theta, Y))+Q*Theta+Pr*alp*S*bett*(Theta-Phi)+Pr*Ec*((1+1/bet)*a1*(diff(W, Y))^2+a5*M*W^2+(1+1/bet)*a1*S1*W^2+S2*W^3+S*betu*(W-U)) = 0;

8.022256200*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^2-.525*d[0]-2.10*h[0]+.525*b[0]+2.20*c[0]+10.23967791*c[1]+2.266560888*c[2]+.525*b[5]*Y^5+.525*b[6]*Y^6+.525*b[7]*Y^7+2.20*c[1]*Y+2.20*c[2]*Y^2+2.20*c[3]*Y^3+2.20*c[4]*Y^4+2.20*c[5]*Y^5+2.20*c[6]*Y^6+2.20*c[7]*Y^7-.525*d[1]*Y-.525*d[2]*Y^2-.525*d[3]*Y^3-.525*d[4]*Y^4+.525*b[1]*Y+.525*b[2]*Y^2+.525*b[3]*Y^3+.525*b[4]*Y^4+16.04451240*(7*Y^6*b[7]+6*Y^5*b[6]+5*Y^4*b[5]+4*Y^3*b[4]+3*Y^2*b[3]+2*Y*b[2]+b[1])^2+5.25*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^3-2.10*Y*h[1]+13.61538438*M*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^2+47.59777865*Y^5*c[7]+33.99841332*Y^4*c[6]+22.66560888*Y^3*c[5]+13.59936533*Y^2*c[4]+6.799682664*Y*c[3]+71.67774537*Y^6*c[7]+61.43806746*Y^5*c[6]+51.19838955*Y^4*c[5]+40.95871164*Y^3*c[4]+30.71903373*Y^2*c[3]+20.47935582*Y*c[2]-2.10*Y^4*h[4]-2.10*Y^3*h[3]-2.10*Y^2*h[2] = 0

(4)

G := Ra*(diff(U, Y))+betu*(W-U) = 0;

2.0*Y^3*d[4]+1.5*Y^2*d[3]+1.0*Y*d[2]+.5*d[1]+.5*b[7]*Y^7+.5*b[6]*Y^6+.5*b[5]*Y^5+.5*b[4]*Y^4-.5*d[4]*Y^4+.5*b[3]*Y^3-.5*d[3]*Y^3+.5*b[2]*Y^2-.5*d[2]*Y^2+.5*b[1]*Y-.5*d[1]*Y+.5*b[0]-.5*d[0] = 0

(5)

H := Ra*(diff(Phi, Y))+bett*(Theta-Phi) = 0;

2.0*Y^3*h[4]+1.5*Y^2*h[3]+1.0*Y*h[2]+.5*h[1]+.5*c[7]*Y^7+.5*c[6]*Y^6+.5*c[5]*Y^5+.5*c[4]*Y^4-.5*Y^4*h[4]+.5*c[3]*Y^3-.5*Y^3*h[3]+.5*c[2]*Y^2-.5*Y^2*h[2]+.5*c[1]*Y-.5*Y*h[1]+.5*c[0]-.5*h[0] = 0

(6)

BCS := (D(W))(0) = 0, (D(Theta))(0) = 0, W(1) = -delta*(1+1/bet)*(D(W))(1), Theta(1) = 1+g*(D(Theta))(1), U(1) = -delta*(1+1/bet)*(D(W))(1), Phi(1) = 1+g*(D(Theta))(1):

W     := unapply(W, Y);Theta := unapply(Theta, Y);
U     := unapply(U, Y);Phi   := unapply(Phi, Y);
F     := unapply(F, Y);T     := unapply(T, Y);
G     := unapply(G, Y);H     := unapply(H, Y);



 

proc (Y) options operator, arrow; Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0] end proc

 

proc (Y) options operator, arrow; Y^7*c[7]+Y^6*c[6]+Y^5*c[5]+Y^4*c[4]+Y^3*c[3]+Y^2*c[2]+Y*c[1]+c[0] end proc

 

proc (Y) options operator, arrow; Y^4*d[4]+Y^3*d[3]+Y^2*d[2]+Y*d[1]+d[0] end proc

 

proc (Y) options operator, arrow; Y^4*h[4]+Y^3*h[3]+Y^2*h[2]+Y*h[1]+h[0] end proc

 

proc (Y) options operator, arrow; 1+0.5e-1*d[0]-0.5e-1*b[0]+.8111904760*b[1]+3.0560976*b[2]+.8111904760*c[0]-.5*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^2-0.5e-1*Y^7*b[7]-0.5e-1*Y^6*b[6]-0.5e-1*Y^5*b[5]-0.5e-1*Y^4*b[4]-0.5e-1*Y^3*b[3]-0.5e-1*Y^2*b[2]-0.5e-1*Y*b[1]+.8111904760*Y^7*c[7]+.8111904760*Y^6*c[6]+.8111904760*Y^5*c[5]+.8111904760*Y^4*c[4]+.8111904760*Y^3*c[3]+.8111904760*Y^2*c[2]+.8111904760*Y*c[1]+0.5e-1*Y^4*d[4]+0.5e-1*Y^3*d[3]+0.5e-1*Y^2*d[2]+0.5e-1*Y*d[1]+1.622380952*Y*b[2]-1.296703274*M*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])+64.1780496*Y^5*b[7]+45.8414640*Y^4*b[6]+30.5609760*Y^3*b[5]+18.3365856*Y^2*b[4]+9.1682928*Y*b[3]+5.678333332*Y^6*b[7]+4.867142856*Y^5*b[6]+4.055952380*Y^4*b[5]+3.244761904*Y^3*b[4]+2.433571428*Y^2*b[3] = 0 end proc

 

proc (Y) options operator, arrow; -.525*d[0]-2.10*h[0]+.525*b[0]+2.20*c[0]+10.23967791*c[1]+2.266560888*c[2]+8.022256200*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^2+16.04451240*(7*Y^6*b[7]+6*Y^5*b[6]+5*Y^4*b[5]+4*Y^3*b[4]+3*Y^2*b[3]+2*Y*b[2]+b[1])^2+5.25*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^3+13.61538438*M*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^2+47.59777865*Y^5*c[7]+33.99841332*Y^4*c[6]+22.66560888*Y^3*c[5]+13.59936533*Y^2*c[4]+6.799682664*Y*c[3]+71.67774537*Y^6*c[7]+61.43806746*Y^5*c[6]+51.19838955*Y^4*c[5]+40.95871164*Y^3*c[4]+30.71903373*Y^2*c[3]+20.47935582*Y*c[2]+.525*Y^7*b[7]+.525*Y^6*b[6]+.525*Y^5*b[5]+.525*Y^4*b[4]+.525*Y^3*b[3]+.525*Y^2*b[2]+.525*Y*b[1]+2.20*Y^7*c[7]+2.20*Y^6*c[6]+2.20*Y^5*c[5]+2.20*Y^4*c[4]+2.20*Y^3*c[3]+2.20*Y^2*c[2]+2.20*Y*c[1]-.525*Y^4*d[4]-.525*Y^3*d[3]-.525*Y^2*d[2]-.525*Y*d[1]-2.10*Y^4*h[4]-2.10*Y^3*h[3]-2.10*Y^2*h[2]-2.10*Y*h[1] = 0 end proc

 

proc (Y) options operator, arrow; 2.0*Y^3*d[4]+1.5*Y^2*d[3]+1.0*Y*d[2]+.5*d[1]+.5*Y^7*b[7]+.5*Y^6*b[6]+.5*Y^5*b[5]+.5*Y^4*b[4]-.5*Y^4*d[4]+.5*Y^3*b[3]-.5*Y^3*d[3]+.5*Y^2*b[2]-.5*Y^2*d[2]+.5*Y*b[1]-.5*Y*d[1]+.5*b[0]-.5*d[0] = 0 end proc

 

proc (Y) options operator, arrow; 2.0*Y^3*h[4]+1.5*Y^2*h[3]+1.0*Y*h[2]+.5*h[1]+.5*Y^7*c[7]+.5*Y^6*c[6]+.5*Y^5*c[5]+.5*Y^4*c[4]-.5*Y^4*h[4]+.5*Y^3*c[3]-.5*Y^3*h[3]+.5*Y^2*c[2]-.5*Y^2*h[2]+.5*Y*c[1]-.5*Y*h[1]+.5*c[0]-.5*h[0] = 0 end proc

(7)

 

z1 := (D(W))(0) = 0:

z2 := W(1) = -delta*(1+1/bet)*(D(W))(1):

z3 := (D(Theta))(0) = 0:

``

z4 := Theta(1) = 1+g*(D(Theta))(1):

 

z5 := U(1) = -delta*(1+1/bet)*(D(W))(1):

``

z6 := Phi(1) = 1+g*(D(Theta))(1):

z7 := F(0):

z8 := F(1):

``

z9 := (D(F))(0):

z10 := (D(F))(1):

z11 := (D(D(F)))(0):

z12 := (D(D(F)))(1):

``

z13 := T(0):

z14 := T(1):

z15 := (D(T))(0):

z16 := (D(T))(1):

z17 := (D(D(T)))(0):

z18 := (D(D(T)))(1):

z19 := G(0):

z20 := G(1):

z21 := H(0):

z22 := H(1):

z23 := (D(G))(0):

z24 := (D(G))(1):

z25 := (D(H))(0):

z26 := (D(H))(1):

``

NULL

NULL

``

Zs :=  [z1, z2, z3, z4, z5, z6, z7, z8, z9, z10, z11, z12, z13, z14,z15, z16,z17,z18,z19,z20,z21, z22, z23, z24,z25, z26]:
print~(Zs):

b[1] = 0

 

b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0] = -0.14e-1*b[7]-0.12e-1*b[6]-0.10e-1*b[5]-0.8e-2*b[4]-0.6e-2*b[3]-0.4e-2*b[2]-0.2e-2*b[1]

 

c[1] = 0

 

c[7]+c[6]+c[5]+c[4]+c[3]+c[2]+c[1]+c[0] = 1+.7*c[7]+.6*c[6]+.5*c[5]+.4*c[4]+.3*c[3]+.2*c[2]+.1*c[1]

 

d[4]+d[3]+d[2]+d[1]+d[0] = -0.14e-1*b[7]-0.12e-1*b[6]-0.10e-1*b[5]-0.8e-2*b[4]-0.6e-2*b[3]-0.4e-2*b[2]-0.2e-2*b[1]

 

h[4]+h[3]+h[2]+h[1]+h[0] = 1+.7*c[7]+.6*c[6]+.5*c[5]+.4*c[4]+.3*c[3]+.2*c[2]+.1*c[1]

 

1.+0.5e-1*d[0]-0.5e-1*b[0]+.8111904760*b[1]+3.0560976*b[2]+.8111904760*c[0]-.5*b[0]^2-1.296703274*M*b[0] = 0

 

1+0.5e-1*d[0]+0.5e-1*d[1]+0.5e-1*d[2]+0.5e-1*d[3]+0.5e-1*d[4]-0.5e-1*b[0]+.7611904760*b[1]+4.628478552*b[2]+11.55186423*b[3]+21.53134750*b[4]+34.56692838*b[5]+50.65860686*b[6]+69.80638293*b[7]+.8111904760*c[0]+.8111904760*c[1]+.8111904760*c[2]+.8111904760*c[3]+.8111904760*c[4]+.8111904760*c[5]+.8111904760*c[6]+.8111904760*c[7]-1.296703274*M*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])-.5*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])^2 = 0

 

0.5e-1*d[1]-0.5e-1*b[1]+1.622380952*b[2]+9.1682928*b[3]+.8111904760*c[1]-1.296703274*M*b[1]-1.0*b[0]*b[1] = 0

 

0.5e-1*d[1]+.10*d[2]+.15*d[3]+.20*d[4]-0.5e-1*b[1]+1.522380952*b[2]+13.88543566*b[3]+46.20745691*b[4]+107.6567375*b[5]+207.4015703*b[6]+354.6102480*b[7]+.8111904760*c[1]+1.622380952*c[2]+2.433571428*c[3]+3.244761904*c[4]+4.055952380*c[5]+4.867142856*c[6]+5.678333332*c[7]-1.296703274*M*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])-1.0*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1]) = 0

 

.10*d[2]-.10*b[2]+4.867142856*b[3]+36.6731712*b[4]+1.622380952*c[2]-2.0*b[0]*b[2]-2.593406548*M*b[2]-1.0*b[1]^2 = 0

 

.10*d[2]+.30*d[3]+.60*d[4]-.10*b[2]+4.567142856*b[3]+55.54174262*b[4]+231.0372846*b[5]+645.9404251*b[6]+1451.810992*b[7]+1.622380952*c[2]+4.867142856*c[3]+9.734285712*c[4]+16.22380952*c[5]+24.33571428*c[6]+34.06999999*c[7]-1.0*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])^2-1.0*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])*(42*b[7]+30*b[6]+20*b[5]+12*b[4]+6*b[3]+2*b[2])-1.296703274*M*(42*b[7]+30*b[6]+20*b[5]+12*b[4]+6*b[3]+2*b[2]) = 0

 

-.525*d[0]-2.10*h[0]+.525*b[0]+2.20*c[0]+10.23967791*c[1]+2.266560888*c[2]+8.022256200*b[0]^2+16.04451240*b[1]^2+5.25*b[0]^3+13.61538438*M*b[0]^2 = 0

 

-2.10*h[4]-.525*d[0]-.525*d[1]-.525*d[2]-.525*d[3]-.525*d[4]-2.10*h[0]-2.10*h[1]-2.10*h[2]-2.10*h[3]+.525*b[0]+.525*b[1]+.525*b[2]+.525*b[3]+.525*b[4]+.525*b[5]+.525*b[6]+.525*b[7]+2.20*c[0]+12.43967791*c[1]+24.94591671*c[2]+39.71871639*c[3]+56.75807697*c[4]+76.06399843*c[5]+97.63648078*c[6]+121.4755240*c[7]+5.25*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])^3+13.61538438*M*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])^2+16.04451240*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])^2+8.022256200*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])^2 = 0

 

-.525*d[1]-2.10*h[1]+.525*b[1]+2.20*c[1]+20.47935582*c[2]+6.799682664*c[3]+15.75*b[0]^2*b[1]+64.17804960*b[1]*b[2]+27.23076876*M*b[0]*b[1]+16.04451240*b[0]*b[1] = 0

 

-8.40*h[4]-.525*d[1]-1.050*d[2]-1.575*d[3]-2.100*d[4]-2.10*h[1]-4.20*h[2]-6.30*h[3]+.525*b[1]+1.050*b[2]+1.575*b[3]+2.100*b[4]+2.625*b[5]+3.150*b[6]+3.675*b[7]+2.20*c[1]+24.87935582*c[2]+74.83775012*c[3]+158.8748656*c[4]+283.7903848*c[5]+456.3839906*c[6]+683.4553654*c[7]+32.08902480*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])*(42*b[7]+30*b[6]+20*b[5]+12*b[4]+6*b[3]+2*b[2])+15.75*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])^2*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])+27.23076876*M*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])+16.04451240*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1]) = 0

 

-1.050*d[2]-4.20*h[2]+1.050*b[2]+4.40*c[2]+61.43806746*c[3]+27.19873066*c[4]+128.3560992*b[2]^2+192.5341488*b[1]*b[3]+31.50*b[0]^2*b[2]+27.23076876*M*b[1]^2+31.50*b[0]*b[1]^2+54.46153752*M*b[0]*b[2]+32.08902480*b[0]*b[2]+16.04451240*b[1]^2 = 0

 

32.08902480*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])*(210*b[7]+120*b[6]+60*b[5]+24*b[4]+6*b[3])+15.75*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])^2*(42*b[7]+30*b[6]+20*b[5]+12*b[4]+6*b[3]+2*b[2])+27.23076876*M*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])^2+31.50*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])^2-25.20*h[4]-1.050*d[2]-3.150*d[3]-6.300*d[4]-4.20*h[2]-12.60*h[3]+1.050*b[2]+3.150*b[3]+6.300*b[4]+10.500*b[5]+15.750*b[6]+22.050*b[7]+4.40*c[2]+74.63806746*c[3]+299.3510005*c[4]+794.3743279*c[5]+1702.742309*c[6]+3194.687934*c[7]+32.08902480*(42*b[7]+30*b[6]+20*b[5]+12*b[4]+6*b[3]+2*b[2])^2+27.23076876*M*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])*(42*b[7]+30*b[6]+20*b[5]+12*b[4]+6*b[3]+2*b[2])+16.04451240*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])^2+16.04451240*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])*(42*b[7]+30*b[6]+20*b[5]+12*b[4]+6*b[3]+2*b[2]) = 0

 

.5*d[1]+.5*b[0]-.5*d[0] = 0

 

1.5*d[4]+1.0*d[3]+.5*d[2]+.5*b[7]+.5*b[6]+.5*b[5]+.5*b[4]+.5*b[3]+.5*b[2]+.5*b[1]+.5*b[0]-.5*d[0] = 0

 

.5*h[1]+.5*c[0]-.5*h[0] = 0

 

1.5*h[4]+1.0*h[3]+.5*h[2]+.5*c[7]+.5*c[6]+.5*c[5]+.5*c[4]+.5*c[3]+.5*c[2]+.5*c[1]+.5*c[0]-.5*h[0] = 0

 

1.0*d[2]+.5*b[1]-.5*d[1] = 0

 

4.0*d[4]+1.5*d[3]+3.5*b[7]+3.0*b[6]+2.5*b[5]+2.0*b[4]+1.5*b[3]+1.0*b[2]+.5*b[1]-.5*d[1] = 0

 

1.0*h[2]+.5*c[1]-.5*h[1] = 0

 

4.0*h[4]+1.5*h[3]+3.5*c[7]+3.0*c[6]+2.5*c[5]+2.0*c[4]+1.5*c[3]+1.0*c[2]+.5*c[1]-.5*h[1] = 0

(8)

IZs := indets(Zs) minus {M};

numelems(Zs), numelems(IZs);

{b[0], b[1], b[2], b[3], b[4], b[5], b[6], b[7], c[0], c[1], c[2], c[3], c[4], c[5], c[6], c[7], d[0], d[1], d[2], d[3], d[4], h[0], h[1], h[2], h[3], h[4]}

 

26, 26

(9)

Zsol := proc(Ma) fsolve(eval(Zs, M=Ma)) end proc:

# example

Zsol(1)

{b[0] = .4528477725, b[1] = 0., b[2] = -.5456992938, b[3] = 0.9787892139e-1, b[4] = 0.5392589969e-1, b[5] = -.1004469452, b[6] = 0.5280749340e-1, b[7] = -0.9643936946e-2, c[0] = 1.688266264, c[1] = 0., c[2] = -2.715762558, c[3] = 8.023205965, c[4] = -18.19126343, c[5] = 22.70336190, c[6] = -13.86047069, c[7] = 3.251216625, d[0] = .2118383626, d[1] = -.2410094099, d[2] = -.1205047050, d[3] = .1233638336, d[4] = 0.2798182974e-1, h[0] = 1.242893264, h[1] = -.4453730000, h[2] = -.2226865000, h[3] = .4041352773, h[4] = -0.8041495798e-1}

(10)

F := proc(Y, Ma) eval(sum(b[i]*Y^i,i=0..7), Zsol(Ma)) end proc:
T := proc(Y, Ma) eval(sum(c[i]*Y^i,i=0..7), Zsol(Ma)) end proc:
G := proc(Y, Ma) eval(sum(d[i]*Y^i,i=0..4), Zsol(Ma)) end proc:
H := proc(Y, Ma) eval(sum(h[i]*Y^i,i=0..4), Zsol(Ma)) end proc:

plot([seq(F(Y, Ma), Ma = 1 .. 3)], Y = 0 .. 1, axes = boxed, color = [red, green, blue], legend = `~`[cat]("M=", [`$`(1 .. 3)]), thickness = 2, labels = [Y, W])

 

plot([seq(T(Y, Ma), Ma = 1 .. 3)], Y = 0 .. 1, axes = boxed, color = [red, green, blue], legend = `~`[cat]("M=", [`$`(1 .. 3)]), thickness = 2, labels = [Y, T])

 

NULL

-1.275852823

(11)

NULL

1.149666788

(12)

``

Cf := 'Cf':
Nu := 'Nu':

plot(
  [seq(a1*(1+1/bet)*diff(F(Y, Ma), Y), Ma = 1 .. 3)]
  , Y=0..1
  , color=[red, green, blue]
  , legend = `~`[cat]("M=", [`$`(1 .. 3)])
  , labels=[Y, "Cf"]
);

print():
plot(
  [seq((-a4-Rd)*diff(T(Y, Ma), Y), Ma = 1 .. 3)]
  , Y=0..1
  , color=[red, green, blue]
  , legend = `~`[cat]("M=", [`$`(1 .. 3)])
  , labels=[Y, "Nu"]
)

 

 

 

NULL

 

 

 

 

``


 

Download AGM-2_mmcdara_1.mw

@ogunmiloro 

Ok.

Your model can be summarized as 

diff(G[dp](T), T) = psi[1]*F[d](T)+psi[2]*L[d](T)+psi[3]*E[c](T)-delta[3]*G[dp](T)
=
diff(G[dp](T), T) = F(T)-delta[3]*G[dp](T);

where
F(T) = psi[1]*F[d](T)+psi[2]*L[d](T)+psi[3]*E[c](T)

You say that F[d], L[d] and E[c] are indeed functions of time.
Let's assume that a correct fitting is such that psi[1], psi[2], psi[3], F[d](T), L[d](T) and E[c](T) get these values

psi[1] = p[1], psi[2] = p[2], psi[3] = p[3], F[d](T)=A(T), L[d](T)=B(T), E[c](T)=C(T)

Thus this solution is such that 

p[1]*A(T) + p[2]*B(T) + p[3]*C(T)

and any 6-uple 

(psi[1], psi[2], psi[3], F[d](T), L[d](T), E[c](T) )

which verifies 

psi[1]*F[d](T)+psi[2]*L[d](T)+psi[3]*E[c](T) = p[1]*A(T) + p[2]*B(T) + p[3]*C(T)

gives an exacyly as good fit than 

psi[1] = p[1], psi[2] = p[2], psi[3] = p[3], F[d](T)=A(T), L[d](T)=B(T), E[c](T)=C(T)

Thus your parameters are not identifiable.

If F[d], L[d] and E[c] arr known functions of time the problem remains the same: if F[d], L[d] and E[c] are not linearly independent there is an infinity of triples (psi[1], psi[2], psi[3]) which give exactly fit, meaning the psi are not identifiable.
Can you give the expressions of F[d], L[d] and E[c]? Or at least there values at the different years?

More of this you said nothing about the initial condition on G[dp](T)


Please enlighten us

  • What does (for instance) year 2022 = 3, 023, 744. 14 mean?
    Are 3 quantities collected each year... or only one (GDP) and do you use commas to separate thousands?
  • So, what are the data?
  • What do you mean by "Fit the following model to the data" ?
    Do you mean "solve the ode
    diff(G[dp](T), T) = psi[1]*F[d](T)+psi[2]*L[d](T)+psi[3]*E[c](T)-delta[3]*G[dp](T);
    for some initial condition
    • Which one?  
      G[dp](some_unknown_date) = 1023744 ?
    • What about the parameters 
      delta[3] = .19, psi[1] = .63, psi[2] = .84, psi[3] = .72
      
      Are these the values to use in the previous ode?
      Or are there only indicative values (initial values) of these four parameters you aimed to assess through some minimization procedure?
    • Are F[d], L[d] andE[c] functions of time or constants?
       

Perhaps you think that the people you're asking for help just have to muddle through your draft, and that any effort on your part is superfluous?

1 2 3 4 5 6 7 Last Page 1 of 107