## 792 Reputation

9 years, 234 days

## First step: define a correct number of I...

You have to fix this big problem before

```# Order of eq1
D_terms  := select(has, indets(eq1, function), D):
Eq1Order := max(degree~(eval(eval(%, f = (eta -> exp(k*eta))), eta=0)));
3

# Total number of IC/BC you impose
ficbc := select(has, {ics, bcs}, f):
NumberIcBc := numelems(%);
4
```

The number of IC/BC is not equal to the order of eq1.
Same thing happen to eq2: a 2nd order ode with 3 IC/BC.

Once you have fix these two problems you will be able (or maybe not...) to solve your differential system:

```# Example 1 : eq1 with 3 ICs and eq2 with 2 ICs
#
# As you then have "free" parameters [a, n, m] use option "parameters".

sol_1 := dsolve(
{eq1, eq2} union select(has, {ics}, {f, theta})
, numeric
, parameters=[a, b, m]
)
proc(x_rkf45)  ...  end;

# Example 2 : eq1 with 2 ICs + 1 BC and eq2 with 2 ICs.
#
# In this case you solve a coupled bv problem and you cannot use the option "parameters".
# You must instanciate [a, b, m] BEFORE calling dsolve, for instance:

data := [a=1, b=2, m=3]:

sys   := eval({eq1, eq2} union {f(0)=0, D(f)(0)=0, D(f)(10)=1, theta(0) = 1, (D(theta))(0) = b}, data):
sol_2 := dsolve(
sys
, numeric
)

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging
```

This error is quite common when you solve bv problems.
There exist several ways to get rid of this error but they are problem-dependent, in particular they may depend here on the values in data. So i will not waste tume to solve a problem wich could not be the one you are interested in.

Last point: writting things like Cp[f] := 4179 (for instance) is not safe (f is also the name of the function ). Look at this

```restart
Cp[f] := 4179
4179
f := x^2;
2
x
Cp[f]
2
Cp[x ]
```

Write Cp__f := 4179 instead (same thing for all the quantities whose names nontain [f])

## remove...

```r := [[5,7],[],[1,3],[],[5,4]];
remove(type, r, [])```

## Convergence!...

F(z) is not convergent if |z| > 1

```restart:
alias(po=pochhammer):
F := z -> sum(po(1,n)^2/po(1/2,n)^2*z^n/(n+1),n=1..infinity):
F(1/2)
/              [3  3   ]  1\
hypergeom|[1, 2, 2, 2], [-, -, 3], -|
\              [2  2   ]  2/
F(2)
infinity
Digits:=10:
eval(F(z), z=2);
evalf[5](%);
/              [3  3   ]   \
4 hypergeom|[1, 2, 2, 2], [-, -, 3], 2|
\              [2  2   ]   /
-3.8666 + 3.5627 I
```

As serie F is not convergent at point z=2 you cannot substutite the value z=2 in the general expression of G(z) because this expression is only valid only if G converges at this point.

A simpler example

```n := 'n':
G := z -> sum(z^n, n=0..infinity);
infinity
-----
\
)     n
z ->   /     z
-----
n = 0
'G(2)' = G(2), 'G(z)' = G(z), 'eval(G(z), z=2)' = eval(G(z), z=2);
1
G(2) = infinity,       G(z) = - ------,       eval(G(z), z = 2) = -1
-1 + z
```

## A "while free" solution...

I found it interesting to show one can find the result without using a while-loop, even though this code is pretty crude and can probably be improved.

```NT := proc(p, q)
local ip, iq, ip1, iq1, div, pf, qf, rf:
ip  := ifactors(p)[2];
iq  := ifactors(q)[2];
ip1 := map2(op, 1, ip):
iq1 := map2(op, 1, iq):
if `subset`({iq1[]}, {ip1[]}) then
div := `intersect`({ip1[]} , {iq1[]});
pf  := select((x -> op(1, x) in div), ip);
qf  := select((x -> op(1, x) in div), iq);
rf  := min(iquo~(map2(op, 2, pf), map2(op, 2, qf)));
printf("%d divides %d  %d times \n", q, p, rf)
else
printf("%d does not divide %d \n", q, p)
end if;
end proc:

NT(294912, 8);
8 divides 294912  5 times

NT(294912, 8*3);
24 divides 294912  2 times

NT(90,3);
3 divides 90  2 times
```

## A simple comment to add to the previous ...

 > restart
 > ineq := log[2](7/10*x) - log[3](3*x-1) <= 0;
 (1)
 > lineq := lhs(ineq); plot(   lineq   , x=1/3..20   , color=blue   , legend=typeset(lineq) )
 > # find values of x such that lineq = 0 s  := solve(lineq); r  := allvalues(s); rf := evalf(r);
 (2)
 > # lineq being continuous for x > 1/3, there is a point x=M between r[1] and r[2] # where lineq reaches its extremal value. diff(lhs(ineq), x); M := solve(%=0); verify(M, r[1]..r[2], interval);
 (3)
 > # At x = M the second ordre derivative of lineq being positive, lineq reaches its minimum: diff(lineq, x\$2); simplify(eval(%, x=M)); is(%, positive)
 (4)
 > # Then lineq < 0 for each x in (r[1], r[2]). # Or, if you prefer: x in RealRange(Open(rf[1]), Open(rf[2])) implies lineq < 0; x in RealRange(rf[1], rf[2]) implies lineq <= 0
 (5)
 >

## Begin to look to functions like simplif...

Begin to look to functions like simplify, normal, factor, collect to cite a few and try to do some things by yourself.

The one to your second question is not and it is impossible that Maple, or any other CAS, can provide the "simplified" expression you want, just because this simplified form is subjective (meaning it is who claims it is a simple form but others could say another form is simpler).
If guide Maple to the result you want you will get it, but Maple cannot finf it by itself

```expr1 := 4*(theta - 1/4)*gamma1^2:
simplify(expr1)
2
(4 theta - 1) gamma1
expr2 := (alpha + gamma1)^2*epsilon^2 + 2*(alpha + gamma1)*(beta + gamma1)*epsilon - 4*alpha*gamma1*theta + beta^2 + 2*beta*gamma1 - 4*(theta - 1/4)*gamma1^2:

Student:-Precalculus:-CompleteSquare(expr2, beta):
op(1, %) + factor(simplify(%-op(1, %)));
2
(beta + (alpha + gamma1) epsilon + gamma1)   - 4 gamma1 theta (alpha + gamma1)
```

solutions.mw

## Wouldn't that be simpler?...

```parms:=[colour=red,thickness=`if`(leader=0, 3, 0),linestyle =dash]:
```

## My opinion, for what it worth...

Personnaly I prefer writting my frequency table (which is not that difficult) instead of using the build-in procedure which I feel quite limited (IMO).

Here is an example of a procedure (it's name F reflects my lack of imagination):
MyFrequencyTable.mw

It acts differently depending wether you want to consider the data are observations of a discrete process or of a continuous one.
You can also define your own header names (for instance using french words)  manage some colors and chose the display precision for percentages.
It should be easy for you to customize it as you want.

## undesired '==' ...

You wrote

`p2b==g2b(x-c*t)`

`p2b=g2b(x-c*t)`

Once corrected you will get what you want.

In case you find that the expressions you get are a little bit obscure or not synthetic enough, this file can give you some ideas for alternative displays forms.mw

## For what it worth...

As a table is a list of equalities indices=entries, it all depends on what you want these indices and entries are.

For instance it could be natural (@Rouben Rostamian  answer) to think that A is the index of  `A=70` and 70 its entry. But what about "DOT=106"?
Is DOT the corresponding index or is it "DOT"?
Is "106" the corresponding entry or is it 106"?

The attached file presents five different ways (among many others while you do not define indices and entries)  to construct a table from your list st.

 > restart:
 > with(StringTools):
 > st := [`A=70`, `B=17`, `C=27`, `D=37`, `E=74`, `F=57`, `G=67`, `H=08`, `I=81`, `J=28`, `K=38`, `L=48`, `M=58`, `N=68`, `O=90`, `P=19`, `Q=29`, `R=39`, `S=49`, `T=59`, `U=96`, `V=010`, `W=110`, `X=210`, `Y=310`, `Z=410`, "SPACE=105", "DOT=106"]
 (1)
 > sst := convert~(st, string):

Option 1

 > map(x -> Split(x,"="), sst): map(x -> x[1]=x[2], %): convert(%, table)
 (2)

Option 2

 > map(x -> Split(x,"="), sst): map(x -> x[1]=parse(x[2]), %): convert(%, table)
 (3)

Option 3

 > convert(parse~(st), table)
 (4)

Option 4

 > map(x -> if x::string then `=`(op(Split(x, "="))) else parse(x) end if, st): convert(%, table)
 (5)

Option 5

 > map(x -> if x::string then op(1, Split(x, "="))=parse(op(2, Split(x, "="))) else parse(x) end if, st): convert(%, table)
 (6)

....

## Exact value...

With Maple 2015.2

 > restart
 > f := evalf(Pi*int(2*299792458^2*662607015*10^(-8)*10^(-34)/((exp(299792458*662607015*10^(-8)*10^(-34)/(1380649*10^(-6)*10^(-23)*lambda*5772.0)) - 1)*lambda^5), lambda = 0 .. infinity));
 (1)
 > g :=Pi*Int(A/((exp(B/lambda) - 1)*lambda^5), lambda = 0 .. infinity);
 (2)
 > value(%) assuming A > 0, B > 0
 (3)
 > evalf(eval(%, [A=2*299792458^2*662607015*10^(-8)*10^(-34), B=299792458*662607015*10^(-8)*10^(-34)/(1380649*10^(-6)*10^(-23)*5772.0)]));
 (4)
 >

In the attached file you will also find an estimation of the error associated to the value (4) given the uncertainties on the Planck and Boltzmann contants and on the solar effective temperature.

## No problem with Maple 2015...

Download and run the attached file to check it works correctly for you.
If it's not the case this is likely a version issue

## A simple typo...

Replace

`SymbolicMatrix := Matrix(6\$3, symbol=m):`

by

`SymbolicMatrix := Matrix(6\$2, symbol=m):`

Do not try to simplify the symbolic solution, it is useless.
Trying to simplify the C_solution is likely a dead end for the expressions are far too complex.

What I suggest you is to do a preliminary work on the equations fn1..fn6 before doing the steps

```MatrixCoeffs := Matrix(6\$2, (i, j) -> coeff(fn||i, C||j)):
RhsCoeffs    := Vector(6, i -> -eval(fn||i, [C1,C2,C3,C4,C5,C6]=~0)):
```

This is not a simple task but I doubt hopping tha build-in Maple functions will prove usefull to get "simple" expressions for C_solution.

As an example a way to simplify the expression of equatiois provided at the end of the attached file
6eqns_mmcdara.mw

M
y last advice is to transform the C_solution into a function of tha many parameters in order to be capable to investigate it as these parameters take different values.

## Any function of a random variable (RV) i...

Any function of a random variable (RV) is named "a statistic".

For instance the PDF is a statistic, and so is the Expectation (commonly named Mean too) or the Variance.
Unfortunately the term average is not a term used to refer to a common statistic. So I have to suppose that speaking of average means speaking of Expectation (~Mean).

To illustrate this, the Expectation/Mean of the face value of a fair dice is 7/2 ( = (1+2+3+4+5+6)/6 )
To obtain the Mean of a RV named sigma with the Statistics package just do

```with(Statistics):
sigma := RandomVariable(...some distribution ...): # see help pages
Mean(sigma)```

You used the function Median which computes, as its name indicates, a statistic named "Median" of a random varriable .
The Median and the Mean have the same value if and only if the probability density function of your (continuous) RV is symetric wrt its mean. This is the case of a Gaussian RV, for instance, but not of an Exponential one.

Then you are wrong writing A0:=Median(sigma) : the correct expression is A0:=Mean(sigma) .

I take it that the notation < Q > means "Expectation of RV Q" isn'it?

So < sigma^2 > means "Expectation of the RV equal to sigma^2".
To get its value just type

`Mean(sigma^2)`

It is true that  < RV >^2 is never equal to < RV^2 > whichever the distribution of the random variable.
The demonstration is extremely simple. Given that the variance of RV is defind by  Var(RV) = < RV^2 > - < RV >^2 ,
then < RV^2 > = < RV >^2 means that Var(RV) = 0 and thus that RV is not a random variable but a "certain variable".

 > restart
 > with(Statistics):
 > sigma := RandomVariable(EmpiricalDistribution([0, 2, 4, 6]));
 (1)
 > A0 := Mean(sigma);
 (2)
 > A := A0^2
 (3)
 > B := Mean(sigma^2)
 (4)
 > Variance(sigma);
 (5)
 > B - A;
 (6)
 >

Final remark:
The quantity < Q^2 >, where Q is some RV, is a statistic name "2nd raw moment of Q" (or "moment of order 2 of Q".

`B := Mean(sigma^2)`

you can simply write

`B := Moment(sigma, 2)`

## A few remarks about integration in a dis...

Point 1:
About your last sentence: "Final goal: select an arbitrary located circular domain for integration and a bivariate distribution with different variance in x and y"
There will be no close form expression of this integral.
To getone you must integrate over ellipses whose half axes are equal to the standard deviations of x and y.
In the more general case of non centered correlated components x and y, the integration has to be done done within elliptic domains centered at the mean x-y point, and whose axes are aligned along the eigenvectors of the variance matrix.

See wiki  item Cumulative distribution function

Point 2:
A main result in Probability Theory (see reference above) is that (I only consider here your simple case of centered uncorrelated components with equal standard deviations sigma) is that the probability that the random vector < X, Y > (whose PDF is your formula (1)) has a square norm less than r__mx2 is related to the ChiSquare distribution with 2 degrees of freedom:

```restart:
with(Statistics):
Z := RandomVariable(ChiSquare(2)):
Probability(Z <= (r__mx/sigma)^2) assuming r__mx > 0, sigma > 0
/        2 \
|   r__mx  |
1 - exp|- --------|
|         2|
\  2 sigma /
```

Point 3:
In case you woulf be interested in integrating a mulnormal PDF over ,o, canonical domains (ellipsoids), it's good to know there exist very specific fast algorithms to do this.

 1 2 3 4 5 Page 1 of 5
﻿