Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 174 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

Function undefined at required values...

tomleslie 13876
February 06 2019
2 0

so required derivatives are undefined.

One may get (mathematically invalid) "approximations" of the required derivatives by considering limits as x->0.01 from abov and below, together with limits as x->0.0418 from above and below.

See the attached


 

> 

restart;
> 

y := (1-Heaviside(x-0.10e-1))*(cos(9.0218219*x)-.99533285*sin(9.0218219*x)-.99999991*cosh(9.0218219*x)+.99533285*sinh(9.0218219*x))+(Heaviside(x-0.10e-1)-Heaviside(x-0.418e-1))*(.24369100*cos(7.7520047*x)-.36109859*sin(7.7520047*x)-.23739778*cosh(7.7520047*x)+.19615343*sinh(7.7520047*x))+Heaviside(x-0.418e-1)*(.95680995*cos(9.0218219*x)-.80884870*sin(9.0218219*x)-1.0381918*cosh(9.0218219*x)+1.1704059*sinh(9.0218219*x));

(1-Heaviside(x-0.10e-1))*(cos(9.0218219*x)-.99533285*sin(9.0218219*x)-.99999991*cosh(9.0218219*x)+.99533285*sinh(9.0218219*x))+(Heaviside(x-0.10e-1)-Heaviside(x-0.418e-1))*(.24369100*cos(7.7520047*x)-.36109859*sin(7.7520047*x)-.23739778*cosh(7.7520047*x)+.19615343*sinh(7.7520047*x))+Heaviside(x-0.418e-1)*(.95680995*cos(9.0218219*x)-.80884870*sin(9.0218219*x)-1.0381918*cosh(9.0218219*x)+1.1704059*sinh(9.0218219*x))

 (1)
> 

#
# Eliminate the Heaviside functions by converting
# the expression to piecewise. Note that the
# resulting expression is UNDEFINED at x=0.01 and
# x=0.0418
#
  yy:=convert(y, piecewise);;

yy := piecewise(x < 0.10e-1, cos(9.0218219*x)-.99533285*sin(9.0218219*x)-.99999991*cosh(9.0218219*x)+.99533285*sinh(9.0218219*x), x = 0.10e-1, Float(undefined), x < 0.418e-1, .24369100*cos(7.7520047*x)-.36109859*sin(7.7520047*x)-.23739778*cosh(7.7520047*x)+.19615343*sinh(7.7520047*x), x = 0.418e-1, Float(undefined), 0.418e-1 < x, .95680995*cos(9.0218219*x)-.80884870*sin(9.0218219*x)-1.0381918*cosh(9.0218219*x)+1.1704059*sinh(9.0218219*x))

 (2)
> 

#
# One can still generate the derivative of the
# piecewsie function - although it will still be
# UNDEFINED at x=0.01 and x=0.0418
#
  diffyy:=diff(yy,x);

diffyy := piecewise(x < 0.1000000000e-1, -9.021821896*sin(9.021821896*x)-8.979715709*cos(9.021821896*x)-9.021821084*sinh(9.021821896*x)+8.979715709*cosh(9.021821896*x), x = 0.1000000000e-1, Float(undefined), x < 0.4180000000e-1, -1.889093776*sin(7.752004694*x)-2.799237963*cos(7.752004694*x)-1.840308705*sinh(7.752004694*x)+1.520582310*cosh(7.752004694*x), x = 0.4180000000e-1, Float(undefined), 0.4180000000e-1 < x, -8.632168958*sin(9.021821896*x)-7.297288914*cos(9.021821896*x)-9.366381519*sinh(9.021821896*x)+10.55919358*cosh(9.021821896*x))

 (3)
> 

#
# Out of idle curiosity, avoid the undefined value
# of the derivative at 0.01 by examining its limit
# as one approaches from the left and the right
#
  limit( op(2, diffyy), x=0.01, left);
  limit( op(6, diffyy), x=0.01, right);

-1.554777372

  

-1.554777657

 (4)
> 

#
# Perform the same limiting operation at x=0.0418
#
  limit(op(6,  diffyy), x=0.0418, left);
  limit(op(10, diffyy), x=0.0418, right);

-2.260736868

  

-2.26073702

 (5)
> 

 

> 

 


 

Download undef.mw

Pragmatic alternative to VV's clever sol...

tomleslie 13876
February 03 2019
0 1

sometimes looking at the problem numerically can guide you towards a viable solution.

See the attached (which for some reason won't display here!!)


numInt.mw

First of all...

tomleslie 13876
February 03 2019
1 2

you should avoid using eq(1), eq(2) etc to indicate distinct equations - eq(1) will probably be  interpreted as the function eq() evaluated with the argument '1' and so on. Similarly for the constants c(1)..c(12). Just don't use round brackets ie '()' for purposes for whihc they were never intended. Only confusion will result.

Editing your worksheet so that eq(1)..eq(12) are replaced by eq1..eq12 and c(2)..c(12) are replaced by c2..c12, then using the (free) add-on DirectSearch(), available from the MAple Application Centre, I can obtain a numerical solution as shown in the following. There may be other solutions, although locating any/all of these would be made much easier if you were able to place some restrictions on the expected ranges for these values.

Check the attached: NB I'm guessing that you do not have the DirectSearch() package installed, so you will not be able to re-execute

> 

restart:
> 

eq1 := 1+c3+c5+c7+c9+c11 = 0:
> 

eq2 := 1.0379456*c2*(omega^2)^(1/4)+1.0379456*c4*(omega^2)^(1/4)+.89185524*c6*(omega^2)^(1/4)+.89185524*c8*(omega^2)^(1/4)+1.0379456*c10*(omega^2)^(1/4)+1.0379456*c12*(omega^2)^(1/4) = 0:
> 

eq3 := 1.1182111*(omega^2)^(3/4)*sin(.27713148*(omega^2)^(1/4))-1.1182111*c2*(omega^2)^(3/4)*cos(.27713148*(omega^2)^(1/4))+1.1182111*c3*(omega^2)^(3/4)*sinh(.27713148*(omega^2)^(1/4))+1.1182111*c4*(omega^2)^(3/4)*cosh(.27713148*(omega^2)^(1/4))+.70938680*c5*(omega^2)^(3/4)*sin(.23812535*(omega^2)^(1/4))-.70938680*c6*(omega^2)^(3/4)*cos(.23812535*(omega^2)^(1/4))+.70938680*c7*(omega^2)^(3/4)*sinh(.23812535*(omega^2)^(1/4))+.70938680*c8*(omega^2)^(3/4)*cosh(.23812535*(omega^2)^(1/4))+1.1182111*c9*(omega^2)^(3/4)*sin(.27713148*(omega^2)^(1/4))-1.1182111*c10*(omega^2)^(3/4)*cos(.27713148*(omega^2)^(1/4))+1.1182111*c11*(omega^2)^(3/4)*sinh(.27713148*(omega^2)^(1/4))+1.1182111*c12*(omega^2)^(3/4)*cosh(.27713148*(omega^2)^(1/4))+.20600541*omega^2*(cos(.27713148*(omega^2)^(1/4))+c2*sin(.27713148*(omega^2)^(1/4))+c3*cosh(.27713148*(omega^2)^(1/4))+c4*sinh(.27713148*(omega^2)^(1/4))+c5*cos(.23812535*(omega^2)^(1/4))+c6*sin(.23812535*(omega^2)^(1/4))+c7*cosh(.23812535*(omega^2)^(1/4))+c8*sinh(.23812535*(omega^2)^(1/4))+c9*cos(.27713148*(omega^2)^(1/4))+c10*sin(.27713148*(omega^2)^(1/4))+c11*cosh(.27713148*(omega^2)^(1/4))+c12*sinh(.27713148*(omega^2)^(1/4))-0.10275662e-1*(omega^2)^(1/4)*sin(.27713148*(omega^2)^(1/4))+0.10275662e-1*c2*(omega^2)^(1/4)*cos(.27713148*(omega^2)^(1/4))+0.10275662e-1*c3*(omega^2)^(1/4)*sinh(.27713148*(omega^2)^(1/4))+0.10275662e-1*c4*(omega^2)^(1/4)*cosh(.27713148*(omega^2)^(1/4))-0.88293670e-2*c5*(omega^2)^(1/4)*sin(.23812535*(omega^2)^(1/4))+0.88293670e-2*c6*(omega^2)^(1/4)*cos(.23812535*(omega^2)^(1/4))+0.88293670e-2*c7*(omega^2)^(1/4)*sinh(.23812535*(omega^2)^(1/4))+0.88293670e-2*c8*(omega^2)^(1/4)*cosh(.23812535*(omega^2)^(1/4))-0.10275662e-1*c9*(omega^2)^(1/4)*sin(.27713148*(omega^2)^(1/4))+0.10275662e-1*c10*(omega^2)^(1/4)*cos(.27713148*(omega^2)^(1/4))+0.10275662e-1*c11*(omega^2)^(1/4)*sinh(.27713148*(omega^2)^(1/4))+0.10275662e-1*c12*(omega^2)^(1/4)*cosh(.27713148*(omega^2)^(1/4))) = 0:
> 

eq4 := 1.1182111*(omega^2)^(3/4)*sin(.27713148*(omega^2)^(1/4))-1.1182111*c2*(omega^2)^(3/4)*cos(.27713148*(omega^2)^(1/4))+1.1182111*c3*(omega^2)^(3/4)*sinh(.27713148*(omega^2)^(1/4))+1.1182111*c4*(omega^2)^(3/4)*cosh(.27713148*(omega^2)^(1/4))+.70938680*c5*(omega^2)^(3/4)*sin(.23812535*(omega^2)^(1/4))-.70938680*c6*(omega^2)^(3/4)*cos(.23812535*(omega^2)^(1/4))+.70938680*c7*(omega^2)^(3/4)*sinh(.23812535*(omega^2)^(1/4))+.70938680*c8*(omega^2)^(3/4)*cosh(.23812535*(omega^2)^(1/4))+1.1182111*c9*(omega^2)^(3/4)*sin(.27713148*(omega^2)^(1/4))-1.1182111*c10*(omega^2)^(3/4)*cos(.27713148*(omega^2)^(1/4))+1.1182111*c11*(omega^2)^(3/4)*sinh(.27713148*(omega^2)^(1/4))+1.1182111*c12*(omega^2)^(3/4)*cosh(.27713148*(omega^2)^(1/4))-.20600541*omega^2*(-0.79735630e-1*(omega^2)^(1/4)*sin(.27713148*(omega^2)^(1/4))+0.79735630e-1*c2*(omega^2)^(1/4)*cos(.27713148*(omega^2)^(1/4))+0.79735630e-1*c3*(omega^2)^(1/4)*sinh(.27713148*(omega^2)^(1/4))+0.79735630e-1*c4*(omega^2)^(1/4)*cosh(.27713148*(omega^2)^(1/4))-0.68512877e-1*c5*(omega^2)^(1/4)*sin(.23812535*(omega^2)^(1/4))+0.68512877e-1*c6*(omega^2)^(1/4)*cos(.23812535*(omega^2)^(1/4))+0.68512877e-1*c7*(omega^2)^(1/4)*sinh(.23812535*(omega^2)^(1/4))+0.68512877e-1*c8*(omega^2)^(1/4)*cosh(.23812535*(omega^2)^(1/4))-0.79735630e-1*c9*(omega^2)^(1/4)*sin(.27713148*(omega^2)^(1/4))+0.79735630e-1*c10*(omega^2)^(1/4)*cos(.27713148*(omega^2)^(1/4))+0.79735630e-1*c11*(omega^2)^(1/4)*sinh(.27713148*(omega^2)^(1/4))+0.79735630e-1*c12*(omega^2)^(1/4)*cosh(.27713148*(omega^2)^(1/4))+0.99000000e-2*cos(.27713148*(omega^2)^(1/4))+0.99000000e-2*c2*sin(.27713148*(omega^2)^(1/4))+0.99000000e-2*c3*cosh(.27713148*(omega^2)^(1/4))+0.99000000e-2*c4*sinh(.27713148*(omega^2)^(1/4))+0.99000000e-2*c5*cos(.23812535*(omega^2)^(1/4))+0.99000000e-2*c6*sin(.23812535*(omega^2)^(1/4))+0.99000000e-2*c7*cosh(.23812535*(omega^2)^(1/4))+0.99000000e-2*c8*sinh(.23812535*(omega^2)^(1/4))+0.99000000e-2*c9*cos(.27713148*(omega^2)^(1/4))+0.99000000e-2*c10*sin(.27713148*(omega^2)^(1/4))+0.99000000e-2*c11*cosh(.27713148*(omega^2)^(1/4))+0.99000000e-2*c12*sinh(.27713148*(omega^2)^(1/4))) = 0:
> 

eq5 := cos(0.51897280e-3*(omega^2)^(1/4))+c2*sin(0.51897280e-3*(omega^2)^(1/4))+c3*cosh(0.51897280e-3*(omega^2)^(1/4))+c4*sinh(0.51897280e-3*(omega^2)^(1/4)) = c5*cos(0.44592762e-3*(omega^2)^(1/4))+c6*sin(0.44592762e-3*(omega^2)^(1/4))+c7*cosh(0.44592762e-3*(omega^2)^(1/4))+c8*sinh(0.44592762e-3*(omega^2)^(1/4)):
> 

eq6 := -1.0379456*(omega^2)^(1/4)*sin(0.51897280e-3*(omega^2)^(1/4))+1.0379456*c2*(omega^2)^(1/4)*cos(0.51897280e-3*(omega^2)^(1/4))+1.0379456*c3*(omega^2)^(1/4)*sinh(0.51897280e-3*(omega^2)^(1/4))+1.0379456*c4*(omega^2)^(1/4)*cosh(0.51897280e-3*(omega^2)^(1/4)) = -.89185524*c5*(omega^2)^(1/4)*sin(0.44592762e-3*(omega^2)^(1/4))+.89185524*c6*(omega^2)^(1/4)*cos(0.44592762e-3*(omega^2)^(1/4))+.89185524*c7*(omega^2)^(1/4)*sinh(0.44592762e-3*(omega^2)^(1/4))+.89185524*c8*(omega^2)^(1/4)*cosh(0.44592762e-3*(omega^2)^(1/4)):
> 

eq7 := -1.0773311*sqrt(omega^2)*cos(0.51897280e-3*(omega^2)^(1/4))-1.0773311*c2*sqrt(omega^2)*sin(0.51897280e-3*(omega^2)^(1/4))+1.0773311*c3*sqrt(omega^2)*cosh(0.51897280e-3*(omega^2)^(1/4))+1.0773311*c4*sqrt(omega^2)*sinh(0.51897280e-3*(omega^2)^(1/4)) = -4.4787661*c5*sqrt(omega^2)*cos(0.44592762e-3*(omega^2)^(1/4))-4.4787661*c6*sqrt(omega^2)*sin(0.44592762e-3*(omega^2)^(1/4))+4.4787661*c7*sqrt(omega^2)*cosh(0.44592762e-3*(omega^2)^(1/4))+4.4787661*c8*sqrt(omega^2)*sinh(0.44592762e-3*(omega^2)^(1/4)):
> 

eq8 := 1.1182110*(omega^2)^(3/4)*sin(0.51897280e-3*(omega^2)^(1/4))-1.1182110*c2*(omega^2)^(3/4)*cos(0.51897280e-3*(omega^2)^(1/4))+1.1182110*c3*(omega^2)^(3/4)*sinh(0.51897280e-3*(omega^2)^(1/4))+1.1182110*c4*(omega^2)^(3/4)*cosh(0.51897280e-3*(omega^2)^(1/4)) = 3.9944110*c5*(omega^2)^(3/4)*sin(0.44592762e-3*(omega^2)^(1/4))-3.9944110*c6*(omega^2)^(3/4)*cos(0.44592762e-3*(omega^2)^(1/4))+3.9944110*c7*(omega^2)^(3/4)*sinh(0.44592762e-3*(omega^2)^(1/4))+3.9944110*c8*(omega^2)^(3/4)*cosh(0.44592762e-3*(omega^2)^(1/4)):
> 

eq9 := c5*cos(0.28806924e-1*(omega^2)^(1/4))+c6*sin(0.28806924e-1*(omega^2)^(1/4))+c7*cosh(0.28806924e-1*(omega^2)^(1/4))+c8*sinh(0.28806924e-1*(omega^2)^(1/4)) = c9*cos(0.33525643e-1*(omega^2)^(1/4))+c10*sin(0.33525643e-1*(omega^2)^(1/4))+c11*cosh(0.33525643e-1*(omega^2)^(1/4))+c12*sinh(0.33525643e-1*(omega^2)^(1/4)):
> 

eq10 := -.89185524*c5*(omega^2)^(1/4)*sin(0.28806924e-1*(omega^2)^(1/4))+.89185524*c6*(omega^2)^(1/4)*cos(0.28806924e-1*(omega^2)^(1/4))+.89185524*c7*(omega^2)^(1/4)*sinh(0.28806924e-1*(omega^2)^(1/4))+.89185524*c8*(omega^2)^(1/4)*cosh(0.28806924e-1*(omega^2)^(1/4)) = -1.0379456*c9*(omega^2)^(1/4)*sin(0.33525643e-1*(omega^2)^(1/4))+1.0379456*c10*(omega^2)^(1/4)*cos(0.33525643e-1*(omega^2)^(1/4))+1.0379456*c11*(omega^2)^(1/4)*sinh(0.33525643e-1*(omega^2)^(1/4))+1.0379456*c12*(omega^2)^(1/4)*cosh(0.33525643e-1*(omega^2)^(1/4)):
> 

eq11 := -1.0773311*c9*sqrt(omega^2)*cos(0.33525643e-1*(omega^2)^(1/4))-1.0773311*c10*sqrt(omega^2)*sin(0.33525643e-1*(omega^2)^(1/4))+1.0773311*c11*sqrt(omega^2)*cosh(0.33525643e-1*(omega^2)^(1/4))+1.0773311*c12*sqrt(omega^2)*sinh(0.33525643e-1*(omega^2)^(1/4)) = -4.4787661*c5*sqrt(omega^2)*cos(0.28806924e-1*(omega^2)^(1/4))-4.4787661*c6*sqrt(omega^2)*sin(0.28806924e-1*(omega^2)^(1/4))+4.4787661*c7*sqrt(omega^2)*cosh(0.28806924e-1*(omega^2)^(1/4))+4.4787661*c8*sqrt(omega^2)*sinh(0.28806924e-1*(omega^2)^(1/4)):
> 

eq12 := 1.1182110*c9*(omega^2)^(3/4)*sin(0.33525643e-1*(omega^2)^(1/4))-1.1182110*c10*(omega^2)^(3/4)*cos(0.33525643e-1*(omega^2)^(1/4))+1.1182110*c11*(omega^2)^(3/4)*sinh(0.33525643e-1*(omega^2)^(1/4))+1.1182110*c12*(omega^2)^(3/4)*cosh(0.33525643e-1*(omega^2)^(1/4)) = 3.9944110*c5*(omega^2)^(3/4)*sin(0.28806924e-1*(omega^2)^(1/4))-3.9944110*c6*(omega^2)^(3/4)*cos(0.28806924e-1*(omega^2)^(1/4))+3.9944110*c7*(omega^2)^(3/4)*sinh(0.28806924e-1*(omega^2)^(1/4))+3.9944110*c8*(omega^2)^(3/4)*cosh(0.28806924e-1*(omega^2)^(1/4)):
> 

#fsolve( [eq1,eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12]):
> 

sol:=DirectSearch:-SolveEquations( [eq1,eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12]);

sol := [4.11865776433866*10^(-10), Vector(12, {(1) = HFloat(1.0762416638554484e-5), (2) = HFloat(1.90806372051692e-6), (3) = HFloat(1.2218119666807776e-20), (4) = HFloat(1.2218119666807311e-20), (5) = HFloat(1.5363809819788995e-5), (6) = HFloat(3.994719428393346e-7), (7) = HFloat(-3.9877655566680003e-13), (8) = HFloat(-5.577535711243513e-21), (9) = HFloat(-7.494238328693537e-6), (10) = HFloat(-1.6006008702439414e-7), (11) = HFloat(-4.037170031660053e-13), (12) = HFloat(9.665940715840068e-21)}), [c10 = -1.14166608668719, c11 = -1.16402644730934, c12 = 6.10899166984361, c2 = 3.95810011086111, c3 = -.999988666477194, c4 = 2.96016686865916, c5 = -2.94430113081394, c6 = 1.84874549443748, c7 = 2.94429710072666, c8 = 2.41428757555186, c9 = 1.16402990629045, omega = -1.39782776397407*10^(-14)], 6689]

 (1)
> 

 

NULL

Download DSprob.mw

As far as I can tell...

tomleslie 13876
February 03 2019
1 0

there is no way to set an alternative currency to be the default.

It is relatively simple to use AddUnit() to define a "new" currency in terms of the default (ie USD).

Entering data with this new unit is then pretty trivial - no simplification or conversion to the default (USD) takes place.

However subsequent 'calculations' appear to always(?) convert to USD, so one has to convert back to the 'new' currency - which is a bit of a pain :-(

> 

restart;
with(Units):
AddUnit(dollar, context=AUD, abbreviations=AUD, conversion=0.75*USD);
BC:=6000*Unit(AUD); #Battery Cost
BS:=3*Unit(kWh/d); #Battery Storage
EC:=.264*Unit(AUD/kWh); #Electricity Cost
BC/EC/(1+.035)^t/BS=t*Unit(yr);
pt:=fsolve(%,t);
FIT:=0.09*Unit(AUD/kWh); #Feed In Tarrif rate
convert( BS*FIT*pt*Unit(yr),'units', 'AUD');
convert( subs(t=pt,convert(EC*BS, units,AUD/yr)*(1+.035)^t*t*Unit(year)),'units','AUD');

Automatically loading the Units[Simple] subpackage

  

6000*Units:-Unit(AUD)

  

3*Units:-Unit(kWh/d)

  

.264*Units:-Unit(AUD/kWh)

  

654545454.7/1.035^t = 31556925.22*t

  

13.18033222

  

0.9e-1*Units:-Unit(AUD/kWh)

  

1299.783619*Units:-Unit(AUD)

  

6000.000001*Units:-Unit(AUD)

 (1)
> 

 

Download unitsAgain.mw

How I'd do it...

tomleslie 13876
January 31 2019
0 0

is shown in the attached - you will have to change the file path/name to something appropriate for your machine

> 

restart;
S:= ExcelTools:-Import("C:/Users/TomLeslie/Desktop/data.xlsx");
seq
( [ seq
    ( S[j, k],
      k = 1..2
    )
  ],
  j = 1..op([1, 1], S)
);

S := Matrix(25, 2, {(1, 1) = 41115.0, (1, 2) = 6.97, (2, 1) = 41257.0, (2, 2) = 18.9, (3, 1) = 41270.0, (3, 2) = 15.67, (4, 1) = 41316.0, (4, 2) = 1.92, (5, 1) = 41344.0, (5, 2) = .7, (6, 1) = 41376.0, (6, 2) = .3, (7, 1) = 41498.0, (7, 2) = 0., (8, 1) = 41620.0, (8, 2) = 0.2e-1, (9, 1) = 41621.0, (9, 2) = 0., (10, 1) = 41751.0, (10, 2) = 0., (11, 1) = 41872.0, (11, 2) = 0., (12, 1) = 41977.0, (12, 2) = 0., (13, 1) = 42118.0, (13, 2) = 0., (14, 1) = 42241.0, (14, 2) = 0., (15, 1) = 42366.0, (15, 2) = 0.4e-1, (16, 1) = 42376.0, (16, 2) = 0.4e-1, (17, 1) = 42557.0, (17, 2) = 2.0, (18, 1) = 42685.0, (18, 2) = 1.2, (19, 1) = 42887.0, (19, 2) = 2.7, (20, 1) = 42961.0, (20, 2) = 1.03, (21, 1) = 43045.0, (21, 2) = .32, (22, 1) = 43153.0, (22, 2) = .12, (23, 1) = 43272.0, (23, 2) = 0.9e-1, (24, 1) = 43396.0, (24, 2) = .64, (25, 1) = 43495.0, (25, 2) = 4.06})

  

[41115.0, 6.97], [41257.0, 18.9], [41270.0, 15.67], [41316.0, 1.92], [41344.0, .7], [41376.0, .3], [41498.0, 0.], [41620.0, 0.2e-1], [41621.0, 0.], [41751.0, 0.], [41872.0, 0.], [41977.0, 0.], [42118.0, 0.], [42241.0, 0.], [42366.0, 0.4e-1], [42376.0, 0.4e-1], [42557.0, 2.0], [42685.0, 1.2], [42887.0, 2.7], [42961.0, 1.03], [43045.0, .32], [43153.0, .12], [43272.0, 0.9e-1], [43396.0, .64], [43495.0, 4.06]

 (1)
> 

``

Download getData.mw

How I'd (probably) do it...

tomleslie 13876
January 30 2019
3 2

would be something like

  restart;
  randomize():
  numThrows:= 10:
  r:= rand(1..2):
  res:=Vector[row](numThrows, i->["H","T"][r()]);

 

My $0.02...

tomleslie 13876
January 30 2019
0 4
  1. Your method of generating samples using the RandomTools() package is very inefficient. It can be made much faster
  2. In order to obtain a valid comparison, you have to check whether "hardware" or "software" floats are being generated and whether there is any hidden conversion between the two types
  3. You should ensure that the resulting datatype (vector, list, whatever) is the same for all methods

I have attempted to address all of the above issues in the attached.

rand() calls routines from the RandomTools() package: it returns individual software floats which have to be assemebled into a vetor

Statistics:-Sample() command appears to always return a vector of hardware floats by default.

RandomTools:-Generate() returns individual software floats which have to be assembled into a vector

RandomTools:-GenerateFloat64() returns individual hardware floats which have to be assembled into a vector

Timings on my machine are

Method                                                  time                 float  datatype
                                                            (secs)            hardware/software

rand()                                                     4.2                     software
Statistics:-Sample()                              0.015                  hardware
RandomTools:-Generate()                    0.218                  software
RandomTools:-GenerateFloat64()        0.172                  hardware

Check the attached

> 

###########################################
# First three execution groups are OP's
# original code
###########################################
> 

restart; with(RandomTools); CodeTools:-Usage(seq((rand(0. .. 1.0))(), i = 1 .. 100000))

memory used=0.50GiB, alloc change=37.00MiB, cpu time=4.70s, real time=4.59s, gc time=296.40ms

  
> 

restart; X := Statistics:-RandomVariable(Uniform(0, 1)); SX := Statistics:-Sample(X); CodeTools:-Usage(SX(100000))

memory used=0.78MiB, alloc change=0 bytes, cpu time=16.00ms, real time=7.00ms, gc time=0ns

  
> 

restart; with(RandomTools); CodeTools:-Usage(seq(Generate(distribution(Uniform(0, 1))), i = 1 .. 100000))

memory used=1.20GiB, alloc change=253.00MiB, cpu time=43.02s, real time=42.60s, gc time=858.01ms

  
> 

##############################################
#
# First method - simple use of rand()
#
# This is "simple" but (relatively) slow. Takes
# about the same time as OP's equivalent
#
##############################################
  restart;
  with(Statistics):
  with(CodeTools):
#
# Comment out the randomize() command if you want
# the same answer each time the code executes
# (which can be useful for debug)
#
  randomize():
  a:= rand(0.0..1.0);
#
# Notice that use of rand() actually returns a call
# to RandomTools:-Generate() !!
#
# Produces a vector of software floats
#
  R:= Usage
      ( Vector
        ( 100000,
          i-> a()
        )
      ):
#
# Histogram the output data - a crude check to
# see how "uniform" it is. Should be 10 bins
# with ~10000 entries in each bin
#
  Histogram
  ( R,
    bincount = 10,
    frequencyscale = absolute
  );
#
# Check the type of data returned - it is a
# "software float"
#
  type(R[1], hfloat);
  type(R[1], sfloat);

proc () options operator, arrow; RandomTools:-Generate(float('range' = 0. .. 1.0, 'method' = 'uniform')) end proc

  

memory used=487.19MiB, alloc change=37.00MiB, cpu time=4.26s, real time=4.12s, gc time=327.60ms

  

 

false

  

true

 (1)
> 

####################################################
#
# Second method - using Statistics:-RandomVariable()
#
# This is
#
# 1. about the same time as OP's equivalent,
# 2. ~250X faster than using simple calls to 'rand()'
# 3. produces a vector of hardware floats by default
#
####################################################
  restart;
  with(Statistics):
  with(CodeTools):
#
# Comment out the randomize() command if you want
# the same answer each time the code executes
# (which can be useful for debug)
#
  randomize():
#
# Define the distribution
#
  X:= RandomVariable(Uniform(0, 1)):
#
# Generate the samples (produces a vector
# of hardware floats by default)
#
  R:= Usage(Sample(X, 100000)):
#
# Histogram the output data - a crude check to
# see how "uniform" it is. Should be 10 bins
# with ~10000 entries in each bin
#
  Histogram( R,
             bincount = 10,
             frequencyscale = absolute
           );
#
# Check the type of data returned - it is a
# "hardware float"
#
  type(R[1], hfloat);
  type(R[1], sfloat);

memory used=0.86MiB, alloc change=0 bytes, cpu time=15.00ms, real time=7.00ms, gc time=0ns

  

 

Vector[row]

  

true

  

false

 (2)
> 

################################################
#
# Third method - using RandomTools:-Generate()
#
# This is
#
# 1. about 200X faster than OP's equivalent!!,
# 2. about 20X faster than simple calls to rand,
# 3, about 15X slower than using the
#    Statistics:RandomVariable approach
# 4. produces software floats by default
#
################################################
  restart;
  with(Statistics):
  with(CodeTools):
  with(RandomTools[MersenneTwister]):
#
# Comment out the randomize() command if you
# the same answer each time the code executes
# (which can be useful for debug)
#
  randomize():
  R:= Usage
      ( Vector
        ( 100000,
          i->GenerateFloat()
        )
      ):
#
# Histogram the output data - a crude check to
# see how "uniform" it is. Should be 10 bins
# with ~10000 entries in each bin
#
  Histogram
  ( R,
    bincount = 10,
    frequencyscale = absolute
  );
#
# Check the type of data returned - it is a
# "software float"
#
  type(R[1], hfloat);
  type(R[1], sfloat);

memory used=28.12MiB, alloc change=39.01MiB, cpu time=218.00ms, real time=215.00ms, gc time=15.60ms

  

 

false

  

true

 (3)
> 

#######################################################
#
# Fourth method - using RandomTools:-GenerateFloat64()
#
# which ought to return hardware floats by default
#
# This is
#
# 1. about 30X faster than simple calls to rand,
# 2. about 10X slower than using the
#    Statistics:RandomVariable approach
# 3. about 1.25X than using RandomTools:-Generate()
# 4. returns hardware floats
#
#######################################################
  restart;
  with(Statistics):
  with(CodeTools):
  with(RandomTools[MersenneTwister]):
#
# Comment out the randomize() command if you
# the same answer each time the code executes
# (which can be useful for debug)
#
  randomize():
#
# Generate the samples and ensure that these are stored
# as hardware floats
#
  R:= Usage
      ( Vector
        ( 100000,
          i-> GenerateFloat64(),
          datatype=float[8]
        )
      ):
#
# Histogram the output data - a crude check to
# see how "uniform" it is. Should be 10 bins
# with ~10000 entries in each bin
#
  Histogram
  ( R,
    bincount = 10,
    frequencyscale = absolute
  );
#
# Check the type of data returned - it is a
# "software float"
#
  type(R[1], hfloat);
  type(R[1], sfloat);

memory used=15.28MiB, alloc change=2.00MiB, cpu time=172.00ms, real time=174.00ms, gc time=0ns

  

 

true

  

false

 (4)
> 

 


 

Download randoms.mw

 

evalf()...

tomleslie 13876
January 28 2019
1 1

See the following

expression := 3*sin(z)+9*sin(x)*(1/4):
values := [x = 1, z = 1.570796327]:
#
# Use
#
evalf(eval(expression,values));
#
# rather than
#
eval( expression,values);
                         

 

Careful with the comparison!!...

tomleslie 13876
January 28 2019
0 0

When comparing the nested loop approach with multiple add() commands in your original, there is a difference between which parts of the code are being executed using exact arithmetic and which are using approximate (aka floating-point) arithmetic.

This completely invalidates the comparison!!

The attached shows the same two approaches where the separation between exact and floating point calculation is set at the same point. The timing difference between the two methods pretty much disappears

As a further comparison, rather than using a mixture of exact and floating-point arithmetic, I have included another couple of execution groups where the nested-loop and multiple add appraoches uses float[8] (aka" double precision") throughout. These both take about the same time: slightly longer than the earlier two execution groups, but somewhat more accurate (because more digits)

See the attached

> 

#
# OP's original nested for loop
#
  restart;
  FF := proc( xx::Array, nn::integer)
              local gg, i, j;
              gg:= 0;
              for i by 3 to nn do
                  for j from i+3 by 3 to nn do
                      gg:= gg+evalf
                              ( 1/((xx[i]-xx[j])^2+(xx[i+1]-xx[j+1])^2+(xx[i+2]-xx[j+2])^2)
                              )
                  end do
              end do;
              return gg
        end proc:
  XX:=Array( 1..9999,
             [ seq
               ( i^2,
                 i = 1..9999
               )
             ]
           ):
  CodeTools:-Usage(FF(XX, 9999));

memory used=1.74GiB, alloc change=1.00MiB, cpu time=12.50s, real time=12.50s, gc time=499.20ms

  

0.1661350362e-2

 (1)
> 

#
# Using nested add() - but with the split between
# exact and floating-point calculations adjusted
# so that it corresponds to the above
#
  restart;
  FF := proc( xx::Array, nn::integer)
              local gg, i, j;
              gg:= add
                   ( add
                     ( evalf
                       ( 1/((xx[i]-xx[j])^2+(xx[i+1]-xx[j+1])^2+(xx[i+2]-xx[j+2])^2)
                       ),
                       j = i+3..nn, 3
                     ),
                     i = 1..nn, 3
                   );
              return gg
        end proc:
  XX := Array( 1..9999,
               [ seq
                 ( i^2,
                   i = 1..9999
                 )
               ]
             ):
  CodeTools:-Usage(evalf(FF(XX, 9999)));

memory used=1.81GiB, alloc change=1.00MiB, cpu time=12.86s, real time=12.90s, gc time=577.20ms

  

0.1661398654e-2

 (2)
> 

#
# Use consistent datatypes thoughout the calculation
# by defining the datatype of the input Array as float[8]
#
# Note that this is equivalent to "double precision", so
# about 13 digits, rather than the (default) 10 digits
# used in the floating point calculations in the above
# execution groups
#
  restart;
  FF := proc( xx::Array, nn::integer)
              local gg, i, j;
              gg:= 0;
              for i by 3 to nn do
                  for j from i+3 by 3 to nn do
                      gg:= gg+1/((xx[i]-xx[j])^2+(xx[i+1]-xx[j+1])^2+(xx[i+2]-xx[j+2])^2)
                  end do
              end do;
              return gg
        end proc:
  XX:=Array( 1..9999,
             [ seq
               ( i^2,
                 i = 1..9999
               )
             ],
             datatype=float[8]
           ):
  CodeTools:-Usage(FF(XX, 9999));``

memory used=1.91GiB, alloc change=2.00MiB, cpu time=16.38s, real time=16.39s, gc time=1.42s

  

HFloat(0.0016613986903789885)

 (3)
> 

  restart;
  FF := proc( xx::Array, nn::integer)
              local gg, i, j;
              gg:= add
                   ( add
                     ( 1/((xx[i]-xx[j])^2+(xx[i+1]-xx[j+1])^2+(xx[i+2]-xx[j+2])^2),
                       j = i+3 .. nn, 3
                     ),
                     i = 1..nn, 3
                   );
              return gg
        end proc:
  XX := Array( 1..9999,
               [ seq
                 ( i^2,
                   i = 1..9999
                 )
               ],
               datatype=float[8]
             ):
  CodeTools:-Usage(evalf(FF(XX, 9999)));

memory used=1.91GiB, alloc change=2.00MiB, cpu time=15.91s, real time=15.93s, gc time=1.34s

  

HFloat(0.0016613986903788584)

 (4)
> 

 

Download loopVsAdd.mw

In the first case...

tomleslie 13876
January 24 2019
1 1

change

end proc():

to

end proc:

in the first procedure definition - I have no idea what the inclusion of '()'  with 'end proc' is intended to do, but Maple interprets this as a syntax error

Question is unclear...

tomleslie 13876
January 24 2019
1 3

so you are probably going to have to clarify.

The attached contains two execution groups.

The first solves (what seems to be??) your ODE system as a straightforward BVP

The second converts this BVP system to an IVP by finding an appropriate value for alpha, so that dA(0)=alpha, forces the original boundary condition uB*dA(2)+E*D(dA)(2)=0 to be satisfied. The latter is then replaced by the former, and the resulting IVP system is solved

> 

  restart;
#
# Define the parameters
#
  kL:=0.005: uB:=0.003: uC:=0.06: E:=0.0005: psi:=100:
#
# Define the ODE system and the BOUNDARY conditions
#
  odeSys:=  diff( dA(t),t,t)=-(uB/E)*diff(dA(t),t)-kL*(cA(t)-dA(t)/psi),
            diff(cA(t),t)=-(kL/uC)*(cA(t)-dA(t)/psi);
  bcs:= cA(0)=100,
        uB*dA(2)+E*D(dA)(2)=0,
        D(dA)(0)=0;
#
# Solve the system
#
  sol:=dsolve( [odeSys, bcs ], numeric);
#
# Plot the dependent variables
#
  plots:-odeplot( sol, [t, dA(t)], t=0..2);
  plots:-odeplot( sol, [t, cA(t)], t=0..2);
#
# Note that one cannot obtain results for t=3 or
# t=-1 since these values are not within the
# defined boundaries
#
# These commands will return errors
#
  sol(-1); sol(3);

diff(diff(dA(t), t), t) = -6.000000000*(diff(dA(t), t))-0.5e-2*cA(t)+0.5000000000e-4*dA(t), diff(cA(t), t) = -0.8333333333e-1*cA(t)+0.8333333333e-3*dA(t)

  

cA(0) = 100, 0.3e-2*dA(2)+0.5e-3*(D(dA))(2) = 0, (D(dA))(0) = 0

  

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(13, {(1) = .0, (2) = .16033428256059937, (3) = .32125650750484075, (4) = .4833165965147155, (5) = .6467671231035825, (6) = .8118027038832505, (7) = .978488846773951, (8) = 1.1468752993061262, (9) = 1.316911901278128, (10) = 1.4886498962831525, (11) = 1.6625924094841709, (12) = 1.835726343725359, (13) = 2.0}, datatype = float[8], order = C_order); Y := Matrix(13, 3, {(1, 1) = 100.0, (1, 2) = .15351698021756233, (1, 3) = .0, (2, 1) = 98.67278757513375, (2, 2) = .14876024609399527, (2, 3) = -0.51092340753756675e-1, (3, 1) = 97.35841916856353, (3, 2) = .1387640717335658, (3, 3) = -0.6997739898854996e-1, (4, 1) = 96.0524493524535, (4, 2) = .12679475330799297, (4, 3) = -0.7651967740484958e-1, (5, 1) = 94.75301693578267, (5, 2) = .1141017902834766, (5, 3) = -0.7832784426111718e-1, (6, 1) = 93.45881674860645, (6, 2) = .10116021525779254, (6, 3) = -0.7833040534069385e-1, (7, 1) = 92.16961366101715, (7, 2) = 0.8815488915403295e-1, (7, 3) = -0.7765063397658895e-1, (8, 1) = 90.88531410920558, (8, 2) = 0.7515628322793266e-1, (8, 3) = -0.767169715317628e-1, (9, 1) = 89.6065856765496, (9, 2) = 0.6219816329811077e-1, (9, 3) = -0.7569195791525911e-1, (10, 1) = 88.33332192553891, (10, 2) = 0.4928956461534264e-1, (10, 3) = -0.7463619088234699e-1, (11, 1) = 87.06215331503489, (11, 2) = 0.36400064455055975e-1, (11, 3) = -0.7356930655391997e-1, (12, 1) = 85.81505796122336, (12, 2) = 0.23753905111906254e-1, (12, 3) = -0.7251807172670598e-1, (13, 1) = 84.64830197885787, (13, 2) = 0.11922175697836164e-1, (13, 3) = -0.7153305418701704e-1}, datatype = float[8], order = C_order); YP := Matrix(13, 3, {(1, 1) = -8.333205402183156, (1, 2) = .0, (1, 3) = -.4999923241509891, (2, 1) = -8.222608330727162, (2, 2) = -0.51092340753756675e-1, (2, 3) = -.18680245534082404, (3, 1) = -8.113085960329325, (3, 2) = -0.6997739898854996e-1, (3, 3) = -0.669207637079312e-1, (4, 1) = -8.00426511675653, (4, 2) = -0.7651967740484958e-1, (4, 3) = -0.21137842595504604e-1, (5, 1) = -7.89598965950748, (5, 2) = -0.7832784426111718e-1, (5, 3) = -0.37923140226961086e-2, (6, 1) = -7.7881504285596295, (6, 2) = -0.7833040534069385e-1, (6, 3) = 0.26934063118937393e-2, (7, 1) = -7.680727675703238, (7, 2) = -0.7765063397658895e-1, (7, 3) = 0.5060143298905584e-2, (8, 1) = -7.57371354522816, (8, 2) = -0.767169715317628e-1, (8, 3) = 0.5879016458710249e-2, (9, 1) = -7.4671636409443645, (9, 2) = -0.7569195791525911e-1, (9, 3) = 0.61219290169715335e-2, (10, 1) = -7.361069085529953, (10, 2) = -0.7463619088234699e-1, (10, 3) = 0.61530001446181695e-2, (11, 1) = -7.255149109242322, (11, 2) = -0.7356930655391997e-1, (11, 3) = 0.6106892751568072e-2, (12, 1) = -7.15123503489497, (12, 2) = -0.7251807172670598e-1, (12, 3) = 0.6034328249374616e-2, (13, 1) = -7.054015229476247, (13, 2) = -0.7153305418701704e-1, (13, 3) = 0.595741133659772e-2}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(13, {(1) = .0, (2) = .16033428256059937, (3) = .32125650750484075, (4) = .4833165965147155, (5) = .6467671231035825, (6) = .8118027038832505, (7) = .978488846773951, (8) = 1.1468752993061262, (9) = 1.316911901278128, (10) = 1.4886498962831525, (11) = 1.6625924094841709, (12) = 1.835726343725359, (13) = 2.0}, datatype = float[8], order = C_order); Y := Matrix(13, 3, {(1, 1) = .0, (1, 2) = -0.4873773704418254e-16, (1, 3) = .0, (2, 1) = -0.1739862157916885e-12, (2, 2) = 0.1874276363286383e-8, (2, 3) = -0.11245673555383011e-7, (3, 1) = -0.5302758854198837e-13, (3, 2) = 0.4496738548389487e-9, (3, 3) = -0.26980460523357126e-8, (4, 1) = -0.47507499663740725e-13, (4, 2) = -0.6834029522918639e-10, (4, 3) = 0.41004336503729474e-9, (5, 1) = 0.4910073964130581e-13, (5, 2) = -0.10398335422273562e-9, (5, 3) = 0.6239007821355648e-9, (6, 1) = -0.2432521807917761e-13, (6, 2) = -0.4697696331034258e-10, (6, 3) = 0.28186235179123627e-9, (7, 1) = -0.16521753618128666e-13, (7, 2) = -0.9159087974075395e-11, (7, 3) = 0.5495494253883809e-10, (8, 1) = 0.3527858334917285e-15, (8, 2) = 0.3885707218657995e-11, (8, 3) = -0.23314126232636988e-10, (9, 1) = 0.4322138925869052e-13, (9, 2) = 0.5126044288679351e-11, (9, 3) = -0.30756221302272794e-10, (10, 1) = 0.4020135487531714e-13, (10, 2) = 0.32077841880083654e-11, (10, 3) = -0.1924655919393402e-10, (11, 1) = -0.6501309876320215e-13, (11, 2) = 0.14288722937576077e-11, (11, 3) = -0.8573006789082293e-11, (12, 1) = -0.1600071522969324e-13, (12, 2) = 0.4327058487863985e-12, (12, 3) = -0.25962099554841803e-11, (13, 1) = -0.8023258475130859e-13, (13, 2) = 0.3160457009305187e-13, (13, 3) = -0.1894656752702962e-12}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[13] elif outpoint = "order" then return 8 elif outpoint = "error" then return HFloat(1.1245673555383011e-8) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [3, 13, [cA(t), dA(t), diff(dA(t), t)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[13] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[13] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(3, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(13, 3, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(3, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(13, 3, X, Y, outpoint, yout, L, V) end if; [t = outpoint, seq('[cA(t), dA(t), diff(dA(t), t)]'[i] = yout[i], i = 1 .. 3)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[13] elif outpoint = "order" then return 8 elif outpoint = "error" then return HFloat(1.1245673555383011e-8) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [3, 13, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[13] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[13] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(13, 3, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(3, {(1) = 0., (2) = 0., (3) = 0.}); `dsolve/numeric/hermite`(13, 3, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 3)] end proc, (2) = Array(0..0, {}), (3) = [t, cA(t), dA(t), diff(dA(t), t)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [t = res[1], seq('[cA(t), dA(t), diff(dA(t), t)]'[i] = res[i+1], i = 1 .. 3)] catch: error end try end proc

  

 

 

Error, (in sol) solution is only defined in the range HFloat(0.0)..HFloat(2.0)

  

Error, (in sol) solution is only defined in the range HFloat(0.0)..HFloat(2.0)

  
> 

  restart;
#
# Set up a procedure which is useful for deriving
# the shooting method
#
  myFun:= proc( alpha, doSol:=0 )
                local kL:=0.005, uB:=0.003, uC:=0.06,
                      E:=0.0005, psi:=100,
                      odeSys, bcs, sol:
                odeSys:=  diff( dA(t),t,t)=-(uB/E)*diff(dA(t),t)-kL*(cA(t)-dA(t)/psi),
                          diff(cA(t),t)=-(kL/uC)*(cA(t)-dA(t)/psi);
                bcs:= cA(0)=100,
                      dA(0)=alpha,
                      D(dA)(0)=0;
                sol:=dsolve( [odeSys, bcs ], numeric);
                if   doSol=0
                then return uB*rhs(sol(2)[3])+E*rhs(sol(2)[4])
                else return sol
                fi;
         end proc:
> 

#
# Obtain an initial-value condition which guarantees
# the OP's boundary condition
#
# uB*dA(2)+E*D(dA)(2)=0
#
# and solve the resulting IVP problem
#
  ans:= myFun(fsolve(myFun), 1):
#
# Check the original boundary condition
#
#  uB*dA(2)+E*D(dA)(2)=0,
#
# This should return 0 (within numerical limits)
#
  eval( uB*rhs(ans(2)[3])+E*rhs(ans(2)[4]),
        [ uB=0.003, E=0.0005]
      );
#
# Note that one can now obtain results for t=3
# and t=-1 if desired, since the problem has been
# comverted from a boundary-value to an initial-value
# one
#
  ans(-1); ans(3);
#
# Plot the dependent variables - note the plot ranges
#
  plots:-odeplot( ans, [t, dA(t)], t=-1..5);
  plots:-odeplot( ans, [t, cA(t)], t=-1..5);

HFloat(-5.423253128045494e-14)

  

[t = -1., cA(t) = HFloat(108.69102019860324), dA(t) = HFloat(-5.426299673464453), diff(dA(t), t) = HFloat(34.00036113392377)]

  

[t = 3., cA(t) = HFloat(77.88017877150217), dA(t) = HFloat(-0.05671212844674486), diff(dA(t), t) = HFloat(-0.06581462188249138)]

  

 

 
> 

 


 

Download shoot.mw

 

Something like...

tomleslie 13876
January 24 2019
0 1

the printout in the attached - whihc seems to be correct except for a numerical problem with Quantile(.., 0)


 

> 

  with(Statistics):
  X := RandomVariable(Normal(0, 1)):
  printf("        "):
  seq( printf( "%8.3f",j), j=0..0.009,0.001);
  printf("\n");
  for k from 0 to 0.9 by 0.1 do
      for i from 0 to 0.09 by 0.01 do     
          printf( "%8.2f", k+i);
          seq( printf( "%8.3f", Quantile(X, k+i+j)), j=0..0.009,0.001);
          printf("\n");
      od;
 od;

           0.000   0.001   0.002   0.003   0.004   0.005   0.006   0.007   0.008   0.009
    0.00 -15.681  -3.090  -2.878  -2.748  -2.652  -2.576  -2.512  -2.457  -2.409  -2.366
    0.01  -2.326  -2.290  -2.257  -2.226  -2.197  -2.170  -2.144  -2.120  -2.097  -2.075
    0.02  -2.054  -2.034  -2.014  -1.995  -1.977  -1.960  -1.943  -1.927  -1.911  -1.896
    0.03  -1.881  -1.866  -1.852  -1.838  -1.825  -1.812  -1.799  -1.787  -1.774  -1.762
    0.04  -1.751  -1.739  -1.728  -1.717  -1.706  -1.695  -1.685  -1.675  -1.665  -1.655
    0.05  -1.645  -1.635  -1.626  -1.616  -1.607  -1.598  -1.589  -1.580  -1.572  -1.563
    0.06  -1.555  -1.546  -1.538  -1.530  -1.522  -1.514  -1.506  -1.499  -1.491  -1.483
    0.07  -1.476  -1.468  -1.461  -1.454  -1.447  -1.440  -1.433  -1.426  -1.419  -1.412
    0.08  -1.405  -1.398  -1.392  -1.385  -1.379  -1.372  -1.366  -1.359  -1.353  -1.347
    0.09  -1.341  -1.335  -1.329  -1.323  -1.317  -1.311  -1.305  -1.299  -1.293  -1.287
    0.10  -1.282  -1.276  -1.270  -1.265  -1.259  -1.254  -1.248  -1.243  -1.237  -1.232
    0.11  -1.227  -1.221  -1.216  -1.211  -1.206  -1.200  -1.195  -1.190  -1.185  -1.180
    0.12  -1.175  -1.170  -1.165  -1.160  -1.155  -1.150  -1.146  -1.141  -1.136  -1.131
    0.13  -1.126  -1.122  -1.117  -1.112  -1.108  -1.103  -1.098  -1.094  -1.089  -1.085
    0.14  -1.080  -1.076  -1.071  -1.067  -1.063  -1.058  -1.054  -1.049  -1.045  -1.041
    0.15  -1.036  -1.032  -1.028  -1.024  -1.019  -1.015  -1.011  -1.007  -1.003  -0.999
    0.16  -0.994  -0.990  -0.986  -0.982  -0.978  -0.974  -0.970  -0.966  -0.962  -0.958
    0.17  -0.954  -0.950  -0.946  -0.942  -0.938  -0.935  -0.931  -0.927  -0.923  -0.919
    0.18  -0.915  -0.912  -0.908  -0.904  -0.900  -0.896  -0.893  -0.889  -0.885  -0.882
    0.19  -0.878  -0.874  -0.871  -0.867  -0.863  -0.860  -0.856  -0.852  -0.849  -0.845
    0.20  -0.842  -0.838  -0.834  -0.831  -0.827  -0.824  -0.820  -0.817  -0.813  -0.810
    0.21  -0.806  -0.803  -0.800  -0.796  -0.793  -0.789  -0.786  -0.782  -0.779  -0.776
    0.22  -0.772  -0.769  -0.765  -0.762  -0.759  -0.755  -0.752  -0.749  -0.745  -0.742
    0.23  -0.739  -0.736  -0.732  -0.729  -0.726  -0.722  -0.719  -0.716  -0.713  -0.710
    0.24  -0.706  -0.703  -0.700  -0.697  -0.693  -0.690  -0.687  -0.684  -0.681  -0.678
    0.25  -0.674  -0.671  -0.668  -0.665  -0.662  -0.659  -0.656  -0.653  -0.650  -0.646
    0.26  -0.643  -0.640  -0.637  -0.634  -0.631  -0.628  -0.625  -0.622  -0.619  -0.616
    0.27  -0.613  -0.610  -0.607  -0.604  -0.601  -0.598  -0.595  -0.592  -0.589  -0.586
    0.28  -0.583  -0.580  -0.577  -0.574  -0.571  -0.568  -0.565  -0.562  -0.559  -0.556
    0.29  -0.553  -0.550  -0.548  -0.545  -0.542  -0.539  -0.536  -0.533  -0.530  -0.527
    0.30  -0.524  -0.522  -0.519  -0.516  -0.513  -0.510  -0.507  -0.504  -0.502  -0.499
    0.31  -0.496  -0.493  -0.490  -0.487  -0.485  -0.482  -0.479  -0.476  -0.473  -0.470
    0.32  -0.468  -0.465  -0.462  -0.459  -0.457  -0.454  -0.451  -0.448  -0.445  -0.443
    0.33  -0.440  -0.437  -0.434  -0.432  -0.429  -0.426  -0.423  -0.421  -0.418  -0.415
    0.34  -0.412  -0.410  -0.407  -0.404  -0.402  -0.399  -0.396  -0.393  -0.391  -0.388
    0.35  -0.385  -0.383  -0.380  -0.377  -0.375  -0.372  -0.369  -0.366  -0.364  -0.361
    0.36  -0.358  -0.356  -0.353  -0.350  -0.348  -0.345  -0.342  -0.340  -0.337  -0.335
    0.37  -0.332  -0.329  -0.327  -0.324  -0.321  -0.319  -0.316  -0.313  -0.311  -0.308
    0.38  -0.305  -0.303  -0.300  -0.298  -0.295  -0.292  -0.290  -0.287  -0.285  -0.282
    0.39  -0.279  -0.277  -0.274  -0.272  -0.269  -0.266  -0.264  -0.261  -0.259  -0.256
    0.40  -0.253  -0.251  -0.248  -0.246  -0.243  -0.240  -0.238  -0.235  -0.233  -0.230
    0.41  -0.228  -0.225  -0.222  -0.220  -0.217  -0.215  -0.212  -0.210  -0.207  -0.204
    0.42  -0.202  -0.199  -0.197  -0.194  -0.192  -0.189  -0.187  -0.184  -0.181  -0.179
    0.43  -0.176  -0.174  -0.171  -0.169  -0.166  -0.164  -0.161  -0.159  -0.156  -0.154
    0.44  -0.151  -0.148  -0.146  -0.143  -0.141  -0.138  -0.136  -0.133  -0.131  -0.128
    0.45  -0.126  -0.123  -0.121  -0.118  -0.116  -0.113  -0.111  -0.108  -0.105  -0.103
    0.46  -0.100  -0.098  -0.095  -0.093  -0.090  -0.088  -0.085  -0.083  -0.080  -0.078
    0.47  -0.075  -0.073  -0.070  -0.068  -0.065  -0.063  -0.060  -0.058  -0.055  -0.053
    0.48  -0.050  -0.048  -0.045  -0.043  -0.040  -0.038  -0.035  -0.033  -0.030  -0.028
    0.49  -0.025  -0.023  -0.020  -0.018  -0.015  -0.013  -0.010  -0.008  -0.005  -0.003
    0.50   0.000   0.003   0.005   0.008   0.010   0.013   0.015   0.018   0.020   0.023
    0.51   0.025   0.028   0.030   0.033   0.035   0.038   0.040   0.043   0.045   0.048
    0.52   0.050   0.053   0.055   0.058   0.060   0.063   0.065   0.068   0.070   0.073
    0.53   0.075   0.078   0.080   0.083   0.085   0.088   0.090   0.093   0.095   0.098
    0.54   0.100   0.103   0.105   0.108   0.111   0.113   0.116   0.118   0.121   0.123
    0.55   0.126   0.128   0.131   0.133   0.136   0.138   0.141   0.143   0.146   0.148
    0.56   0.151   0.154   0.156   0.159   0.161   0.164   0.166   0.169   0.171   0.174
    0.57   0.176   0.179   0.181   0.184   0.187   0.189   0.192   0.194   0.197   0.199
    0.58   0.202   0.204   0.207   0.210   0.212   0.215   0.217   0.220   0.222   0.225
    0.59   0.228   0.230   0.233   0.235   0.238   0.240   0.243   0.246   0.248   0.251
    0.60   0.253   0.256   0.259   0.261   0.264   0.266   0.269   0.272   0.274   0.277
    0.61   0.279   0.282   0.285   0.287   0.290   0.292   0.295   0.298   0.300   0.303
    0.62   0.305   0.308   0.311   0.313   0.316   0.319   0.321   0.324   0.327   0.329
    0.63   0.332   0.335   0.337   0.340   0.342   0.345   0.348   0.350   0.353   0.356
    0.64   0.358   0.361   0.364   0.366   0.369   0.372   0.375   0.377   0.380   0.383
    0.65   0.385   0.388   0.391   0.393   0.396   0.399   0.402   0.404   0.407   0.410
    0.66   0.412   0.415   0.418   0.421   0.423   0.426   0.429   0.432   0.434   0.437
    0.67   0.440   0.443   0.445   0.448   0.451   0.454   0.457   0.459   0.462   0.465
    0.68   0.468   0.470   0.473   0.476   0.479   0.482   0.485   0.487   0.490   0.493
    0.69   0.496   0.499   0.502   0.504   0.507   0.510   0.513   0.516   0.519   0.522
    0.70   0.524   0.527   0.530   0.533   0.536   0.539   0.542   0.545   0.548   0.550
    0.71   0.553   0.556   0.559   0.562   0.565   0.568   0.571   0.574   0.577   0.580
    0.72   0.583   0.586   0.589   0.592   0.595   0.598   0.601   0.604   0.607   0.610
    0.73   0.613   0.616   0.619   0.622   0.625   0.628   0.631   0.634   0.637   0.640
    0.74   0.643   0.646   0.650   0.653   0.656   0.659   0.662   0.665   0.668   0.671
    0.75   0.674   0.678   0.681   0.684   0.687   0.690   0.693   0.697   0.700   0.703
    0.76   0.706   0.710   0.713   0.716   0.719   0.722   0.726   0.729   0.732   0.736
    0.77   0.739   0.742   0.745   0.749   0.752   0.755   0.759   0.762   0.765   0.769
    0.78   0.772   0.776   0.779   0.782   0.786   0.789   0.793   0.796   0.800   0.803
    0.79   0.806   0.810   0.813   0.817   0.820   0.824   0.827   0.831   0.834   0.838
    0.80   0.842   0.845   0.849   0.852   0.856   0.860   0.863   0.867   0.871   0.874
    0.81   0.878   0.882   0.885   0.889   0.893   0.896   0.900   0.904   0.908   0.912
    0.82   0.915   0.919   0.923   0.927   0.931   0.935   0.938   0.942   0.946   0.950
    0.83   0.954   0.958   0.962   0.966   0.970   0.974   0.978   0.982   0.986   0.990
    0.84   0.994   0.999   1.003   1.007   1.011   1.015   1.019   1.024   1.028   1.032
    0.85   1.036   1.041   1.045   1.049   1.054   1.058   1.063   1.067   1.071   1.076
    0.86   1.080   1.085   1.089   1.094   1.098   1.103   1.108   1.112   1.117   1.122
    0.87   1.126   1.131   1.136   1.141   1.146   1.150   1.155   1.160   1.165   1.170
    0.88   1.175   1.180   1.185   1.190   1.195   1.200   1.206   1.211   1.216   1.221
    0.89   1.227   1.232   1.237   1.243   1.248   1.254   1.259   1.265   1.270   1.276
    0.90   1.282   1.287   1.293   1.299   1.305   1.311   1.317   1.323   1.329   1.335
    0.91   1.341   1.347   1.353   1.359   1.366   1.372   1.379   1.385   1.392   1.398
    0.92   1.405   1.412   1.419   1.426   1.433   1.440   1.447   1.454   1.461   1.468
    0.93   1.476   1.483   1.491   1.499   1.506   1.514   1.522   1.530   1.538   1.546
    0.94   1.555   1.563   1.572   1.580   1.589   1.598   1.607   1.616   1.626   1.635
    0.95   1.645   1.655   1.665   1.675   1.685   1.695   1.706   1.717   1.728   1.739
    0.96   1.751   1.762   1.774   1.787   1.799   1.812   1.825   1.838   1.852   1.866
    0.97   1.881   1.896   1.911   1.927   1.943   1.960   1.977   1.995   2.014   2.034
    0.98   2.054   2.075   2.097   2.120   2.144   2.170   2.197   2.226   2.257   2.290
    0.99   2.326   2.366   2.409   2.457   2.512   2.576   2.652   2.748   2.878   3.090

  
> 

 


 

Download quantiles.mw

Nothing much wrong with what you have...

tomleslie 13876
January 22 2019
1 3

but awkward if you want to change either the number of dice or the number of tries.

Consider the attached, which is more "flexible"

> 

  restart;
> 

#
# Define a procedure which does the work
#
  rolls:= proc( nD::posint, nT::posint )
                uses Statistics:
                local X, vals, probs;
              #
              # randomize the seed for the random number
              # generators. Comment this out if the same
              # answer is required each time the worksheet
              # is executed - which Can be useful for debug!
              #
                randomize():
              #
              # Define appropriate number of random
              # variables
              #
                X:= seq
                    ( RandomVariable(DiscreteUniform(1, 6)),
                      i = 1 .. nD
                    );
              #
              # Carry out 'nT' trials, compute the
              # frequencies, and sort by the die-sum
              # values
              #
                vals:= [ seq
                         ( [lhs(j), rhs(j)],
                           j in sort
                                ( Tally
                                  ( Sample
                                    ( add(X),
                                      nT
                                    )
                                  ),
                                  (i,j)->lhs(i)<lhs(j)
                                )
                         )
                       ];
              #
              # Compute probabilities for each die-sum value
              #
                probs:= [ seq
                          ( [ j[1],100.0*j[2]/nT ],
                            j in vals
                          )
                        ];
                return vals, probs
       end proc:
#
# R|un this procedure for any number of dice
# and trials, using
#
# rolls( numberOfDie, numberOfTrials);
#
# for example
#
   nV, pV:=rolls( 5, 100000);
#
# Plot the probabilities for each
# sum-value
#
  plot(pV, style=point);

[[HFloat(5.0), 16], [HFloat(6.0), 59], [HFloat(7.0), 186], [HFloat(8.0), 438], [HFloat(9.0), 916], [HFloat(10.0), 1648], [HFloat(11.0), 2624], [HFloat(12.0), 3883], [HFloat(13.0), 5410], [HFloat(14.0), 6834], [HFloat(15.0), 8305], [HFloat(16.0), 9551], [HFloat(17.0), 10032], [HFloat(18.0), 10132], [HFloat(19.0), 9464], [HFloat(20.0), 8368], [HFloat(21.0), 6917], [HFloat(22.0), 5520], [HFloat(23.0), 3855], [HFloat(24.0), 2620], [HFloat(25.0), 1624], [HFloat(26.0), 866], [HFloat(27.0), 462], [HFloat(28.0), 202], [HFloat(29.0), 53], [HFloat(30.0), 15]], [[HFloat(5.0), 0.1600000000e-1], [HFloat(6.0), 0.5900000000e-1], [HFloat(7.0), .1860000000], [HFloat(8.0), .4380000000], [HFloat(9.0), .9160000000], [HFloat(10.0), 1.648000000], [HFloat(11.0), 2.624000000], [HFloat(12.0), 3.883000000], [HFloat(13.0), 5.410000000], [HFloat(14.0), 6.834000000], [HFloat(15.0), 8.305000000], [HFloat(16.0), 9.551000000], [HFloat(17.0), 10.03200000], [HFloat(18.0), 10.13200000], [HFloat(19.0), 9.464000000], [HFloat(20.0), 8.368000000], [HFloat(21.0), 6.917000000], [HFloat(22.0), 5.520000000], [HFloat(23.0), 3.855000000], [HFloat(24.0), 2.620000000], [HFloat(25.0), 1.624000000], [HFloat(26.0), .8660000000], [HFloat(27.0), .4620000000], [HFloat(28.0), .2020000000], [HFloat(29.0), 0.5300000000e-1], [HFloat(30.0), 0.1500000000e-1]]

  

 
> 

 

``


 

Download rollDice.mw

 

More an observation really...

tomleslie 13876
January 22 2019
0 0

I tried the DirectSearch() add-on with this problem. First attempt bailed out because the function evaluation limit was exceeded, so I used the output of this first calculation as an 'initialpoint', increased the evaluation limit drastically (and went to bed)

I should have put on a timer on the second attempt - but I forgot: I'm guessing (from the number of iterations) that it took a couple of hours to come up with an "answer" which isn't incredibly good (see attached): the residuals are significant  compared to some of the variable values (particularly h0).

Can't think of any way to do better other than restricting the allowable range(s) of the required parameters.

See the attached


 

> 

gleichungssystem1 := {h[0] = 1-(sum(h_0[i]*((w[i]-T[i]-w[0]+T[0])/(c_0[i]-c_0[0]))^eta[i], i = 1 .. 13)), (T[1]-T[0])/(w[1]-T[1]-w[0]+T[0]) = (h_0[1]*((w[1]-T[1]-w[0]+T[0])/(c_0[1]-c_0[0]))^eta[1]*(1-1/(p*(w[1]-T[1])^nu)-eta[1]*(T[1]-T[0])/(w[1]-T[1]-w[0]+T[0]))+h_0[2]*((w[2]-T[2]-w[0]+T[0])/(c_0[2]-c_0[0]))^eta[2]*(1-1/(p*(w[2]-T[2])^nu)-eta[2]*(T[2]-T[0])/(w[2]-T[2]-w[0]+T[0]))+h_0[3]*((w[3]-T[3]-w[0]+T[0])/(c_0[3]-c_0[0]))^eta[3]*(1-1/(p*(w[3]-T[3])^nu)-eta[3]*(T[3]-T[0])/(w[3]-T[3]-w[0]+T[0]))+h_0[4]*((w[4]-T[4]-w[0]+T[0])/(c_0[4]-c_0[0]))^eta[4]*(1-1/(p*(w[4]-T[4])^nu)-eta[4]*(T[4]-T[0])/(w[4]-T[4]-w[0]+T[0]))+h_0[5]*((w[5]-T[5]-w[0]+T[0])/(c_0[5]-c_0[0]))^eta[5]*(1-1/(p*(w[5]-T[5])^nu)-eta[5]*(T[5]-T[0])/(w[5]-T[5]-w[0]+T[0]))+h_0[6]*((w[6]-T[6]-w[0]+T[0])/(c_0[6]-c_0[0]))^eta[6]*(1-1/(p*(w[6]-T[6])^nu)-eta[6]*(T[6]-T[0])/(w[6]-T[6]-w[0]+T[0]))+h_0[7]*((w[7]-T[7]-w[0]+T[0])/(c_0[7]-c_0[0]))^eta[7]*(1-1/(p*(w[7]-T[7])^nu)-eta[7]*(T[7]-T[0])/(w[7]-T[7]-w[0]+T[0]))+h_0[8]*((w[8]-T[8]-w[0]+T[0])/(c_0[8]-c_0[0]))^eta[8]*(1-1/(p*(w[8]-T[8])^nu)-eta[8]*(T[8]-T[0])/(w[8]-T[8]-w[0]+T[0]))+h_0[9]*((w[9]-T[9]-w[0]+T[0])/(c_0[9]-c_0[0]))^eta[9]*(1-1/(p*(w[9]-T[9])^nu)-eta[9]*(T[9]-T[0])/(w[9]-T[9]-w[0]+T[0]))+h_0[10]*((w[10]-T[10]-w[0]+T[0])/(c_0[10]-c_0[0]))^eta[10]*(1-1/(p*(w[10]-T[10])^nu)-eta[10]*(T[10]-T[0])/(w[10]-T[10]-w[0]+T[0]))+h_0[11]*((w[11]-T[11]-w[0]+T[0])/(c_0[11]-c_0[0]))^eta[11]*(1-1/(p*(w[11]-T[11])^nu)-eta[11]*(T[11]-T[0])/(w[11]-T[11]-w[0]+T[0]))+h_0[12]*((w[12]-T[12]-w[0]+T[0])/(c_0[12]-c_0[0]))^eta[12]*(1-1/(p*(w[12]-T[12])^nu)-eta[12]*(T[12]-T[0])/(w[12]-T[12]-w[0]+T[0]))+h_0[13]*((w[13]-T[13]-w[0]+T[0])/(c_0[13]-c_0[0]))^eta[13]*(1-1/(p*(w[13]-T[13])^nu)-eta[13]*(T[13]-T[0])/(w[13]-T[13]-w[0]+T[0])))/(Zeta[1]*h_0[1]*((w[1]-T[1]-w[0]+T[0])/(c_0[1]-c_0[0]))^eta[1]), (T[2]-T[1])/(w[2]-T[2]-w[1]+T[1]) = (h_0[2]*((w[2]-T[2]-w[0]+T[0])/(c_0[2]-c_0[0]))^eta[2]*(1-1/(p*(w[2]-T[2])^nu)-eta[2]*(T[2]-T[0])/(w[2]-T[2]-w[0]+T[0]))+h_0[3]*((w[3]-T[3]-w[0]+T[0])/(c_0[3]-c_0[0]))^eta[3]*(1-1/(p*(w[3]-T[3])^nu)-eta[3]*(T[3]-T[0])/(w[3]-T[3]-w[0]+T[0]))+h_0[4]*((w[4]-T[4]-w[0]+T[0])/(c_0[4]-c_0[0]))^eta[4]*(1-1/(p*(w[4]-T[4])^nu)-eta[4]*(T[4]-T[0])/(w[4]-T[4]-w[0]+T[0]))+h_0[5]*((w[5]-T[5]-w[0]+T[0])/(c_0[5]-c_0[0]))^eta[5]*(1-1/(p*(w[5]-T[5])^nu)-eta[5]*(T[5]-T[0])/(w[5]-T[5]-w[0]+T[0]))+h_0[6]*((w[6]-T[6]-w[0]+T[0])/(c_0[6]-c_0[0]))^eta[6]*(1-1/(p*(w[6]-T[6])^nu)-eta[6]*(T[6]-T[0])/(w[6]-T[6]-w[0]+T[0]))+h_0[7]*((w[7]-T[7]-w[0]+T[0])/(c_0[7]-c_0[0]))^eta[7]*(1-1/(p*(w[7]-T[7])^nu)-eta[7]*(T[7]-T[0])/(w[7]-T[7]-w[0]+T[0]))+h_0[8]*((w[8]-T[8]-w[0]+T[0])/(c_0[8]-c_0[0]))^eta[8]*(1-1/(p*(w[8]-T[8])^nu)-eta[8]*(T[8]-T[0])/(w[8]-T[8]-w[0]+T[0]))+h_0[9]*((w[9]-T[9]-w[0]+T[0])/(c_0[9]-c_0[0]))^eta[9]*(1-1/(p*(w[9]-T[9])^nu)-eta[9]*(T[9]-T[0])/(w[9]-T[9]-w[0]+T[0]))+h_0[10]*((w[10]-T[10]-w[0]+T[0])/(c_0[10]-c_0[0]))^eta[10]*(1-1/(p*(w[10]-T[10])^nu)-eta[10]*(T[10]-T[0])/(w[10]-T[10]-w[0]+T[0]))+h_0[11]*((w[11]-T[11]-w[0]+T[0])/(c_0[11]-c_0[0]))^eta[11]*(1-1/(p*(w[11]-T[11])^nu)-eta[11]*(T[11]-T[0])/(w[11]-T[11]-w[0]+T[0]))+h_0[12]*((w[12]-T[12]-w[0]+T[0])/(c_0[12]-c_0[0]))^eta[12]*(1-1/(p*(w[12]-T[12])^nu)-eta[12]*(T[12]-T[0])/(w[12]-T[12]-w[0]+T[0]))+h_0[13]*((w[13]-T[13]-w[0]+T[0])/(c_0[13]-c_0[0]))^eta[13]*(1-1/(p*(w[13]-T[13])^nu)-eta[13]*(T[13]-T[0])/(w[13]-T[13]-w[0]+T[0])))/(Zeta[2]*h_0[2]*((w[2]-T[2]-w[0]+T[0])/(c_0[2]-c_0[0]))^eta[2]), (T[3]-T[2])/(w[3]-T[3]-w[2]+T[2]) = (h_0[3]*((w[3]-T[3]-w[0]+T[0])/(c_0[3]-c_0[0]))^eta[3]*(1-1/(p*(w[3]-T[3])^nu)-eta[3]*(T[3]-T[0])/(w[3]-T[3]-w[0]+T[0]))+h_0[4]*((w[4]-T[4]-w[0]+T[0])/(c_0[4]-c_0[0]))^eta[4]*(1-1/(p*(w[4]-T[4])^nu)-eta[4]*(T[4]-T[0])/(w[4]-T[4]-w[0]+T[0]))+h_0[5]*((w[5]-T[5]-w[0]+T[0])/(c_0[5]-c_0[0]))^eta[5]*(1-1/(p*(w[5]-T[5])^nu)-eta[5]*(T[5]-T[0])/(w[5]-T[5]-w[0]+T[0]))+h_0[6]*((w[6]-T[6]-w[0]+T[0])/(c_0[6]-c_0[0]))^eta[6]*(1-1/(p*(w[6]-T[6])^nu)-eta[6]*(T[6]-T[0])/(w[6]-T[6]-w[0]+T[0]))+h_0[7]*((w[7]-T[7]-w[0]+T[0])/(c_0[7]-c_0[0]))^eta[7]*(1-1/(p*(w[7]-T[7])^nu)-eta[7]*(T[7]-T[0])/(w[7]-T[7]-w[0]+T[0]))+h_0[8]*((w[8]-T[8]-w[0]+T[0])/(c_0[8]-c_0[0]))^eta[8]*(1-1/(p*(w[8]-T[8])^nu)-eta[8]*(T[8]-T[0])/(w[8]-T[8]-w[0]+T[0]))+h_0[9]*((w[9]-T[9]-w[0]+T[0])/(c_0[9]-c_0[0]))^eta[9]*(1-1/(p*(w[9]-T[9])^nu)-eta[9]*(T[9]-T[0])/(w[9]-T[9]-w[0]+T[0]))+h_0[10]*((w[10]-T[10]-w[0]+T[0])/(c_0[10]-c_0[0]))^eta[10]*(1-1/(p*(w[10]-T[10])^nu)-eta[10]*(T[10]-T[0])/(w[10]-T[10]-w[0]+T[0]))+h_0[11]*((w[11]-T[11]-w[0]+T[0])/(c_0[11]-c_0[0]))^eta[11]*(1-1/(p*(w[11]-T[11])^nu)-eta[11]*(T[11]-T[0])/(w[11]-T[11]-w[0]+T[0]))+h_0[12]*((w[12]-T[12]-w[0]+T[0])/(c_0[12]-c_0[0]))^eta[12]*(1-1/(p*(w[12]-T[12])^nu)-eta[12]*(T[12]-T[0])/(w[12]-T[12]-w[0]+T[0]))+h_0[13]*((w[13]-T[13]-w[0]+T[0])/(c_0[13]-c_0[0]))^eta[13]*(1-1/(p*(w[13]-T[13])^nu)-eta[13]*(T[13]-T[0])/(w[13]-T[13]-w[0]+T[0])))/(Zeta[3]*h_0[3]*((w[3]-T[3]-w[0]+T[0])/(c_0[3]-c_0[0]))^eta[3]), (T[4]-T[3])/(w[4]-T[4]-w[3]+T[3]) = (h_0[4]*((w[4]-T[4]-w[0]+T[0])/(c_0[4]-c_0[0]))^eta[4]*(1-1/(p*(w[4]-T[4])^nu)-eta[4]*(T[4]-T[0])/(w[4]-T[4]-w[0]+T[0]))+h_0[5]*((w[5]-T[5]-w[0]+T[0])/(c_0[5]-c_0[0]))^eta[5]*(1-1/(p*(w[5]-T[5])^nu)-eta[5]*(T[5]-T[0])/(w[5]-T[5]-w[0]+T[0]))+h_0[6]*((w[6]-T[6]-w[0]+T[0])/(c_0[6]-c_0[0]))^eta[6]*(1-1/(p*(w[6]-T[6])^nu)-eta[6]*(T[6]-T[0])/(w[6]-T[6]-w[0]+T[0]))+h_0[7]*((w[7]-T[7]-w[0]+T[0])/(c_0[7]-c_0[0]))^eta[7]*(1-1/(p*(w[7]-T[7])^nu)-eta[7]*(T[7]-T[0])/(w[7]-T[7]-w[0]+T[0]))+h_0[8]*((w[8]-T[8]-w[0]+T[0])/(c_0[8]-c_0[0]))^eta[8]*(1-1/(p*(w[8]-T[8])^nu)-eta[8]*(T[8]-T[0])/(w[8]-T[8]-w[0]+T[0]))+h_0[9]*((w[9]-T[9]-w[0]+T[0])/(c_0[9]-c_0[0]))^eta[9]*(1-1/(p*(w[9]-T[9])^nu)-eta[9]*(T[9]-T[0])/(w[9]-T[9]-w[0]+T[0]))+h_0[10]*((w[10]-T[10]-w[0]+T[0])/(c_0[10]-c_0[0]))^eta[10]*(1-1/(p*(w[10]-T[10])^nu)-eta[10]*(T[10]-T[0])/(w[10]-T[10]-w[0]+T[0]))+h_0[11]*((w[11]-T[11]-w[0]+T[0])/(c_0[11]-c_0[0]))^eta[11]*(1-1/(p*(w[11]-T[11])^nu)-eta[11]*(T[11]-T[0])/(w[11]-T[11]-w[0]+T[0]))+h_0[12]*((w[12]-T[12]-w[0]+T[0])/(c_0[12]-c_0[0]))^eta[12]*(1-1/(p*(w[12]-T[12])^nu)-eta[12]*(T[12]-T[0])/(w[12]-T[12]-w[0]+T[0]))+h_0[13]*((w[13]-T[13]-w[0]+T[0])/(c_0[13]-c_0[0]))^eta[13]*(1-1/(p*(w[13]-T[13])^nu)-eta[13]*(T[13]-T[0])/(w[13]-T[13]-w[0]+T[0])))/(Zeta[4]*h_0[4]*((w[4]-T[4]-w[0]+T[0])/(c_0[4]-c_0[0]))^eta[4]), (T[5]-T[4])/(w[5]-T[5]-w[4]+T[4]) = (h_0[5]*((w[5]-T[5]-w[0]+T[0])/(c_0[5]-c_0[0]))^eta[5]*(1-1/(p*(w[5]-T[5])^nu)-eta[5]*(T[5]-T[0])/(w[5]-T[5]-w[0]+T[0]))+h_0[6]*((w[6]-T[6]-w[0]+T[0])/(c_0[6]-c_0[0]))^eta[6]*(1-1/(p*(w[6]-T[6])^nu)-eta[6]*(T[6]-T[0])/(w[6]-T[6]-w[0]+T[0]))+h_0[7]*((w[7]-T[7]-w[0]+T[0])/(c_0[7]-c_0[0]))^eta[7]*(1-1/(p*(w[7]-T[7])^nu)-eta[7]*(T[7]-T[0])/(w[7]-T[7]-w[0]+T[0]))+h_0[8]*((w[8]-T[8]-w[0]+T[0])/(c_0[8]-c_0[0]))^eta[8]*(1-1/(p*(w[8]-T[8])^nu)-eta[8]*(T[8]-T[0])/(w[8]-T[8]-w[0]+T[0]))+h_0[9]*((w[9]-T[9]-w[0]+T[0])/(c_0[9]-c_0[0]))^eta[9]*(1-1/(p*(w[9]-T[9])^nu)-eta[9]*(T[9]-T[0])/(w[9]-T[9]-w[0]+T[0]))+h_0[10]*((w[10]-T[10]-w[0]+T[0])/(c_0[10]-c_0[0]))^eta[10]*(1-1/(p*(w[10]-T[10])^nu)-eta[10]*(T[10]-T[0])/(w[10]-T[10]-w[0]+T[0]))+h_0[11]*((w[11]-T[11]-w[0]+T[0])/(c_0[11]-c_0[0]))^eta[11]*(1-1/(p*(w[11]-T[11])^nu)-eta[11]*(T[11]-T[0])/(w[11]-T[11]-w[0]+T[0]))+h_0[12]*((w[12]-T[12]-w[0]+T[0])/(c_0[12]-c_0[0]))^eta[12]*(1-1/(p*(w[12]-T[12])^nu)-eta[12]*(T[12]-T[0])/(w[12]-T[12]-w[0]+T[0]))+h_0[13]*((w[13]-T[13]-w[0]+T[0])/(c_0[13]-c_0[0]))^eta[13]*(1-1/(p*(w[13]-T[13])^nu)-eta[13]*(T[13]-T[0])/(w[13]-T[13]-w[0]+T[0])))/(Zeta[5]*h_0[5]*((w[5]-T[5]-w[0]+T[0])/(c_0[5]-c_0[0]))^eta[5]), (T[6]-T[5])/(w[6]-T[6]-w[5]+T[5]) = (h_0[6]*((w[6]-T[6]-w[0]+T[0])/(c_0[6]-c_0[0]))^eta[6]*(1-1/(p*(w[6]-T[6])^nu)-eta[6]*(T[6]-T[0])/(w[6]-T[6]-w[0]+T[0]))+h_0[7]*((w[7]-T[7]-w[0]+T[0])/(c_0[7]-c_0[0]))^eta[7]*(1-1/(p*(w[7]-T[7])^nu)-eta[7]*(T[7]-T[0])/(w[7]-T[7]-w[0]+T[0]))+h_0[8]*((w[8]-T[8]-w[0]+T[0])/(c_0[8]-c_0[0]))^eta[8]*(1-1/(p*(w[8]-T[8])^nu)-eta[8]*(T[8]-T[0])/(w[8]-T[8]-w[0]+T[0]))+h_0[9]*((w[9]-T[9]-w[0]+T[0])/(c_0[9]-c_0[0]))^eta[9]*(1-1/(p*(w[9]-T[9])^nu)-eta[9]*(T[9]-T[0])/(w[9]-T[9]-w[0]+T[0]))+h_0[10]*((w[10]-T[10]-w[0]+T[0])/(c_0[10]-c_0[0]))^eta[10]*(1-1/(p*(w[10]-T[10])^nu)-eta[10]*(T[10]-T[0])/(w[10]-T[10]-w[0]+T[0]))+h_0[11]*((w[11]-T[11]-w[0]+T[0])/(c_0[11]-c_0[0]))^eta[11]*(1-1/(p*(w[11]-T[11])^nu)-eta[11]*(T[11]-T[0])/(w[11]-T[11]-w[0]+T[0]))+h_0[12]*((w[12]-T[12]-w[0]+T[0])/(c_0[12]-c_0[0]))^eta[12]*(1-1/(p*(w[12]-T[12])^nu)-eta[12]*(T[12]-T[0])/(w[12]-T[12]-w[0]+T[0]))+h_0[13]*((w[13]-T[13]-w[0]+T[0])/(c_0[13]-c_0[0]))^eta[13]*(1-1/(p*(w[13]-T[13])^nu)-eta[13]*(T[13]-T[0])/(w[13]-T[13]-w[0]+T[0])))/(Zeta[6]*h_0[6]*((w[6]-T[6]-w[0]+T[0])/(c_0[6]-c_0[0]))^eta[6]), (T[7]-T[6])/(w[7]-T[7]-w[6]+T[6]) = (h_0[7]*((w[7]-T[7]-w[0]+T[0])/(c_0[7]-c_0[0]))^eta[7]*(1-1/(p*(w[7]-T[7])^nu)-eta[7]*(T[7]-T[0])/(w[7]-T[7]-w[0]+T[0]))+h_0[8]*((w[8]-T[8]-w[0]+T[0])/(c_0[8]-c_0[0]))^eta[8]*(1-1/(p*(w[8]-T[8])^nu)-eta[8]*(T[8]-T[0])/(w[8]-T[8]-w[0]+T[0]))+h_0[9]*((w[9]-T[9]-w[0]+T[0])/(c_0[9]-c_0[0]))^eta[9]*(1-1/(p*(w[9]-T[9])^nu)-eta[9]*(T[9]-T[0])/(w[9]-T[9]-w[0]+T[0]))+h_0[10]*((w[10]-T[10]-w[0]+T[0])/(c_0[10]-c_0[0]))^eta[10]*(1-1/(p*(w[10]-T[10])^nu)-eta[10]*(T[10]-T[0])/(w[10]-T[10]-w[0]+T[0]))+h_0[11]*((w[11]-T[11]-w[0]+T[0])/(c_0[11]-c_0[0]))^eta[11]*(1-1/(p*(w[11]-T[11])^nu)-eta[11]*(T[11]-T[0])/(w[11]-T[11]-w[0]+T[0]))+h_0[12]*((w[12]-T[12]-w[0]+T[0])/(c_0[12]-c_0[0]))^eta[12]*(1-1/(p*(w[12]-T[12])^nu)-eta[12]*(T[12]-T[0])/(w[12]-T[12]-w[0]+T[0]))+h_0[13]*((w[13]-T[13]-w[0]+T[0])/(c_0[13]-c_0[0]))^eta[13]*(1-1/(p*(w[13]-T[13])^nu)-eta[13]*(T[13]-T[0])/(w[13]-T[13]-w[0]+T[0])))/(Zeta[7]*h_0[7]*((w[7]-T[7]-w[0]+T[0])/(c_0[7]-c_0[0]))^eta[7]), (T[8]-T[7])/(w[8]-T[8]-w[7]+T[7]) = (h_0[8]*((w[8]-T[8]-w[0]+T[0])/(c_0[8]-c_0[0]))^eta[8]*(1-1/(p*(w[8]-T[8])^nu)-eta[8]*(T[8]-T[0])/(w[8]-T[8]-w[0]+T[0]))+h_0[9]*((w[9]-T[9]-w[0]+T[0])/(c_0[9]-c_0[0]))^eta[9]*(1-1/(p*(w[9]-T[9])^nu)-eta[9]*(T[9]-T[0])/(w[9]-T[9]-w[0]+T[0]))+h_0[10]*((w[10]-T[10]-w[0]+T[0])/(c_0[10]-c_0[0]))^eta[10]*(1-1/(p*(w[10]-T[10])^nu)-eta[10]*(T[10]-T[0])/(w[10]-T[10]-w[0]+T[0]))+h_0[11]*((w[11]-T[11]-w[0]+T[0])/(c_0[11]-c_0[0]))^eta[11]*(1-1/(p*(w[11]-T[11])^nu)-eta[11]*(T[11]-T[0])/(w[11]-T[11]-w[0]+T[0]))+h_0[12]*((w[12]-T[12]-w[0]+T[0])/(c_0[12]-c_0[0]))^eta[12]*(1-1/(p*(w[12]-T[12])^nu)-eta[12]*(T[12]-T[0])/(w[12]-T[12]-w[0]+T[0]))+h_0[13]*((w[13]-T[13]-w[0]+T[0])/(c_0[13]-c_0[0]))^eta[13]*(1-1/(p*(w[13]-T[13])^nu)-eta[13]*(T[13]-T[0])/(w[13]-T[13]-w[0]+T[0])))/(Zeta[8]*h_0[8]*((w[8]-T[8]-w[0]+T[0])/(c_0[8]-c_0[0]))^eta[8]), (T[9]-T[8])/(w[9]-T[9]-w[8]+T[8]) = (h_0[9]*((w[9]-T[9]-w[0]+T[0])/(c_0[9]-c_0[0]))^eta[9]*(1-1/(p*(w[9]-T[9])^nu)-eta[9]*(T[9]-T[0])/(w[9]-T[9]-w[0]+T[0]))+h_0[10]*((w[10]-T[10]-w[0]+T[0])/(c_0[10]-c_0[0]))^eta[10]*(1-1/(p*(w[10]-T[10])^nu)-eta[10]*(T[10]-T[0])/(w[10]-T[10]-w[0]+T[0]))+h_0[11]*((w[11]-T[11]-w[0]+T[0])/(c_0[11]-c_0[0]))^eta[11]*(1-1/(p*(w[11]-T[11])^nu)-eta[11]*(T[11]-T[0])/(w[11]-T[11]-w[0]+T[0]))+h_0[12]*((w[12]-T[12]-w[0]+T[0])/(c_0[12]-c_0[0]))^eta[12]*(1-1/(p*(w[12]-T[12])^nu)-eta[12]*(T[12]-T[0])/(w[12]-T[12]-w[0]+T[0]))+h_0[13]*((w[13]-T[13]-w[0]+T[0])/(c_0[13]-c_0[0]))^eta[13]*(1-1/(p*(w[13]-T[13])^nu)-eta[13]*(T[13]-T[0])/(w[13]-T[13]-w[0]+T[0])))/(Zeta[9]*h_0[9]*((w[9]-T[9]-w[0]+T[0])/(c_0[9]-c_0[0]))^eta[9]), (T[10]-T[9])/(w[10]-T[10]-w[9]+T[9]) = (h_0[10]*((w[10]-T[10]-w[0]+T[0])/(c_0[10]-c_0[0]))^eta[10]*(1-1/(p*(w[10]-T[10])^nu)-eta[10]*(T[10]-T[0])/(w[10]-T[10]-w[0]+T[0]))+h_0[11]*((w[11]-T[11]-w[0]+T[0])/(c_0[11]-c_0[0]))^eta[11]*(1-1/(p*(w[11]-T[11])^nu)-eta[11]*(T[11]-T[0])/(w[11]-T[11]-w[0]+T[0]))+h_0[12]*((w[12]-T[12]-w[0]+T[0])/(c_0[12]-c_0[0]))^eta[12]*(1-1/(p*(w[12]-T[12])^nu)-eta[12]*(T[12]-T[0])/(w[12]-T[12]-w[0]+T[0]))+h_0[13]*((w[13]-T[13]-w[0]+T[0])/(c_0[13]-c_0[0]))^eta[13]*(1-1/(p*(w[13]-T[13])^nu)-eta[13]*(T[13]-T[0])/(w[13]-T[13]-w[0]+T[0])))/(Zeta[10]*h_0[10]*((w[10]-T[10]-w[0]+T[0])/(c_0[10]-c_0[0]))^eta[10]), (T[11]-T[10])/(w[11]-T[11]-w[10]+T[10]) = (h_0[11]*((w[11]-T[11]-w[0]+T[0])/(c_0[11]-c_0[0]))^eta[11]*(1-1/(p*(w[11]-T[11])^nu)-eta[11]*(T[11]-T[0])/(w[11]-T[11]-w[0]+T[0]))+h_0[12]*((w[12]-T[12]-w[0]+T[0])/(c_0[12]-c_0[0]))^eta[12]*(1-1/(p*(w[12]-T[12])^nu)-eta[12]*(T[12]-T[0])/(w[12]-T[12]-w[0]+T[0]))+h_0[13]*((w[13]-T[13]-w[0]+T[0])/(c_0[13]-c_0[0]))^eta[13]*(1-1/(p*(w[13]-T[13])^nu)-eta[13]*(T[13]-T[0])/(w[13]-T[13]-w[0]+T[0])))/(Zeta[11]*h_0[11]*((w[11]-T[11]-w[0]+T[0])/(c_0[11]-c_0[0]))^eta[11]), (T[12]-T[11])/(w[12]-T[12]-w[11]+T[11]) = (h_0[12]*((w[12]-T[12]-w[0]+T[0])/(c_0[12]-c_0[0]))^eta[12]*(1-1/(p*(w[12]-T[12])^nu)-eta[12]*(T[12]-T[0])/(w[12]-T[12]-w[0]+T[0]))+h_0[13]*((w[13]-T[13]-w[0]+T[0])/(c_0[13]-c_0[0]))^eta[13]*(1-1/(p*(w[13]-T[13])^nu)-eta[13]*(T[13]-T[0])/(w[13]-T[13]-w[0]+T[0])))/(Zeta[12]*h_0[12]*((w[12]-T[12]-w[0]+T[0])/(c_0[12]-c_0[0]))^eta[12]), (T[13]-T[12])/(w[13]-T[13]-w[12]+T[12]) = h_0[13]*((w[13]-T[13]-w[0]+T[0])/(c_0[13]-c_0[0]))^eta[13]*(1-1/(p*(w[13]-T[13])^nu)-eta[13]*(T[13]-T[0])/(w[13]-T[13]-w[0]+T[0]))*(1/(Zeta[13]*h_0[13]*((w[13]-T[13]-w[0]+T[0])/(c_0[13]-c_0[0]))^eta[13])), h[0]*T[0]+sum(h_0[i]*((w[i]-T[i]-w[0]+T[0])/(c_0[i]-c_0[0]))^eta[i]*T[i], i = 1 .. 13)-H = 0, h[0]/(p*(w[0]-T[0]))+sum(h_0[i]*((w[i]-T[i]-w[0]+T[0])/(c_0[i]-c_0[0]))^eta[i]/(p*(w[i]-T[i])), i = 1 .. 13)-1 = 0}
> 

Inputparameter := w[0] = 0, eta[0] = .5, Zeta[0] = 0, h_0[0] = .142, c_0[0] = 6000, w[1] = 2000, eta[1] = .5, Zeta[1] = .25, h_0[1] = 0.33e-1, c_0[1] = 7200, w[2] = 4000, eta[2] = .5, Zeta[2] = .5, h_0[2] = 0.27e-1, c_0[2] = 8400, w[3] = 6000, eta[3] = .5, Zeta[3] = .75, h_0[3] = 0.28e-1, c_0[3] = 9600, w[4] = 8000, eta[4] = .5, Zeta[4] = 1, h_0[4] = 0.3e-1, c_0[4] = 10800, w[5] = 10000, eta[5] = .5, Zeta[5] = 1.25, h_0[5] = 0.48e-1, c_0[5] = 12000, w[6] = 12500, eta[6] = .5, Zeta[6] = 1.25, h_0[6] = 0.52e-1, c_0[6] = 13500, w[7] = 15000, eta[7] = .5, Zeta[7] = 1.5, h_0[7] = 0.65e-1, c_0[7] = 15000, w[8] = 17500, eta[8] = .5, Zeta[8] = 1.75, h_0[8] = 0.47e-1, c_0[8] = 16500, w[9] = 20000, eta[9] = 0, Zeta[9] = 2, h_0[9] = 0.82e-1, c_0[9] = 18000, w[10] = 25000, eta[10] = 0, Zeta[10] = 1.25, h_0[10] = 0.98e-1, c_0[10] = 21000, w[11] = 30000, eta[11] = 0, Zeta[11] = 1.5, h_0[11] = .164, c_0[11] = 24000, w[12] = 50000, eta[12] = 0, Zeta[12] = .625, h_0[12] = .145, c_0[12] = 36000, w[13] = 100000, eta[13] = 0, Zeta[13] = .5, h_0[13] = 0.39e-1, c_0[13] = 66000, H = 5000, nu = 1

w[0] = 0, eta[0] = .5, Zeta[0] = 0, h_0[0] = .142, c_0[0] = 6000, w[1] = 2000, eta[1] = .5, Zeta[1] = .25, h_0[1] = 0.33e-1, c_0[1] = 7200, w[2] = 4000, eta[2] = .5, Zeta[2] = .5, h_0[2] = 0.27e-1, c_0[2] = 8400, w[3] = 6000, eta[3] = .5, Zeta[3] = .75, h_0[3] = 0.28e-1, c_0[3] = 9600, w[4] = 8000, eta[4] = .5, Zeta[4] = 1, h_0[4] = 0.3e-1, c_0[4] = 10800, w[5] = 10000, eta[5] = .5, Zeta[5] = 1.25, h_0[5] = 0.48e-1, c_0[5] = 12000, w[6] = 12500, eta[6] = .5, Zeta[6] = 1.25, h_0[6] = 0.52e-1, c_0[6] = 13500, w[7] = 15000, eta[7] = .5, Zeta[7] = 1.5, h_0[7] = 0.65e-1, c_0[7] = 15000, w[8] = 17500, eta[8] = .5, Zeta[8] = 1.75, h_0[8] = 0.47e-1, c_0[8] = 16500, w[9] = 20000, eta[9] = 0, Zeta[9] = 2, h_0[9] = 0.82e-1, c_0[9] = 18000, w[10] = 25000, eta[10] = 0, Zeta[10] = 1.25, h_0[10] = 0.98e-1, c_0[10] = 21000, w[11] = 30000, eta[11] = 0, Zeta[11] = 1.5, h_0[11] = .164, c_0[11] = 24000, w[12] = 50000, eta[12] = 0, Zeta[12] = .625, h_0[12] = .145, c_0[12] = 36000, w[13] = 100000, eta[13] = 0, Zeta[13] = .5, h_0[13] = 0.39e-1, c_0[13] = 66000, H = 5000, nu = 1

 (1)
> 

gleichungssystem2 := subs({Inputparameter}, gleichungssystem1)

{h[0] = .472-0.33e-1*(5/3-(1/1200)*T[1]+(1/1200)*T[0])^.5-0.27e-1*(5/3-(1/2400)*T[2]+(1/2400)*T[0])^.5-0.28e-1*(5/3-(1/3600)*T[3]+(1/3600)*T[0])^.5-0.3e-1*(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5-0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5-0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5-0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5-0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5, (T[1]-T[0])/(2000-T[1]+T[0]) = 121.2121212*(0.33e-1*(5/3-(1/1200)*T[1]+(1/1200)*T[0])^.5*(1-1/(p*(2000-T[1]))-.5*(T[1]-T[0])/(2000-T[1]+T[0]))+0.27e-1*(5/3-(1/2400)*T[2]+(1/2400)*T[0])^.5*(1-1/(p*(4000-T[2]))-.5*(T[2]-T[0])/(4000-T[2]+T[0]))+0.28e-1*(5/3-(1/3600)*T[3]+(1/3600)*T[0])^.5*(1-1/(p*(6000-T[3]))-.5*(T[3]-T[0])/(6000-T[3]+T[0]))+0.3e-1*(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5*(1-1/(p*(8000-T[4]))-.5*(T[4]-T[0])/(8000-T[4]+T[0]))+0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5*(1-1/(p*(10000-T[5]))-.5*(T[5]-T[0])/(10000-T[5]+T[0]))+0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5*(1-1/(p*(12500-T[6]))-.5*(T[6]-T[0])/(12500-T[6]+T[0]))+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*(1-1/(p*(15000-T[7]))-.5*(T[7]-T[0])/(15000-T[7]+T[0]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/1200)*T[1]+(1/1200)*T[0])^.5, (T[2]-T[1])/(2000-T[2]+T[1]) = 74.07407407*(0.27e-1*(5/3-(1/2400)*T[2]+(1/2400)*T[0])^.5*(1-1/(p*(4000-T[2]))-.5*(T[2]-T[0])/(4000-T[2]+T[0]))+0.28e-1*(5/3-(1/3600)*T[3]+(1/3600)*T[0])^.5*(1-1/(p*(6000-T[3]))-.5*(T[3]-T[0])/(6000-T[3]+T[0]))+0.3e-1*(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5*(1-1/(p*(8000-T[4]))-.5*(T[4]-T[0])/(8000-T[4]+T[0]))+0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5*(1-1/(p*(10000-T[5]))-.5*(T[5]-T[0])/(10000-T[5]+T[0]))+0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5*(1-1/(p*(12500-T[6]))-.5*(T[6]-T[0])/(12500-T[6]+T[0]))+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*(1-1/(p*(15000-T[7]))-.5*(T[7]-T[0])/(15000-T[7]+T[0]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/2400)*T[2]+(1/2400)*T[0])^.5, (T[3]-T[2])/(2000-T[3]+T[2]) = 47.61904761*(0.28e-1*(5/3-(1/3600)*T[3]+(1/3600)*T[0])^.5*(1-1/(p*(6000-T[3]))-.5*(T[3]-T[0])/(6000-T[3]+T[0]))+0.3e-1*(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5*(1-1/(p*(8000-T[4]))-.5*(T[4]-T[0])/(8000-T[4]+T[0]))+0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5*(1-1/(p*(10000-T[5]))-.5*(T[5]-T[0])/(10000-T[5]+T[0]))+0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5*(1-1/(p*(12500-T[6]))-.5*(T[6]-T[0])/(12500-T[6]+T[0]))+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*(1-1/(p*(15000-T[7]))-.5*(T[7]-T[0])/(15000-T[7]+T[0]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/3600)*T[3]+(1/3600)*T[0])^.5, (T[4]-T[3])/(2000-T[4]+T[3]) = 33.33333333*(0.3e-1*(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5*(1-1/(p*(8000-T[4]))-.5*(T[4]-T[0])/(8000-T[4]+T[0]))+0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5*(1-1/(p*(10000-T[5]))-.5*(T[5]-T[0])/(10000-T[5]+T[0]))+0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5*(1-1/(p*(12500-T[6]))-.5*(T[6]-T[0])/(12500-T[6]+T[0]))+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*(1-1/(p*(15000-T[7]))-.5*(T[7]-T[0])/(15000-T[7]+T[0]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5, (T[5]-T[4])/(2000-T[5]+T[4]) = 16.66666667*(0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5*(1-1/(p*(10000-T[5]))-.5*(T[5]-T[0])/(10000-T[5]+T[0]))+0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5*(1-1/(p*(12500-T[6]))-.5*(T[6]-T[0])/(12500-T[6]+T[0]))+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*(1-1/(p*(15000-T[7]))-.5*(T[7]-T[0])/(15000-T[7]+T[0]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5, (T[6]-T[5])/(2500-T[6]+T[5]) = 15.38461538*(0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5*(1-1/(p*(12500-T[6]))-.5*(T[6]-T[0])/(12500-T[6]+T[0]))+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*(1-1/(p*(15000-T[7]))-.5*(T[7]-T[0])/(15000-T[7]+T[0]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5, (T[7]-T[6])/(2500-T[7]+T[6]) = 10.25641026*(0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*(1-1/(p*(15000-T[7]))-.5*(T[7]-T[0])/(15000-T[7]+T[0]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5, (T[8]-T[7])/(2500-T[8]+T[7]) = 12.15805471*(0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5, (T[9]-T[8])/(2500-T[9]+T[8]) = 3.219512195-.5000000000/(p*(20000-T[9]))-.5975609756/(p*(25000-T[10]))-1.000000000/(p*(30000-T[11]))-.8841463415/(p*(50000-T[12]))-.2378048781/(p*(100000-T[13])), (T[10]-T[9])/(5000-T[10]+T[9]) = 3.640816326-.8000000000/(p*(25000-T[10]))-1.338775510/(p*(30000-T[11]))-1.183673469/(p*(50000-T[12]))-.3183673469/(p*(100000-T[13])), (T[11]-T[10])/(5000-T[11]+T[10]) = 1.414634147-.6666666668/(p*(30000-T[11]))-.5894308944/(p*(50000-T[12]))-.1585365854/(p*(100000-T[13])), (T[12]-T[11])/(20000-T[12]+T[11]) = 2.030344828-1.600000000/(p*(50000-T[12]))-.4303448276/(p*(100000-T[13])), (T[13]-T[12])/(50000-T[13]+T[12]) = 2.000000000-2.000000000/(p*(100000-T[13])), h[0]*T[0]+0.33e-1*(5/3-(1/1200)*T[1]+(1/1200)*T[0])^.5*T[1]+0.27e-1*(5/3-(1/2400)*T[2]+(1/2400)*T[0])^.5*T[2]+0.28e-1*(5/3-(1/3600)*T[3]+(1/3600)*T[0])^.5*T[3]+0.3e-1*(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5*T[4]+0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5*T[5]+0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5*T[6]+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*T[7]+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*T[8]+0.82e-1*T[9]+0.98e-1*T[10]+.164*T[11]+.145*T[12]+0.39e-1*T[13]-5000 = 0, -h[0]/(p*T[0])+0.33e-1*(5/3-(1/1200)*T[1]+(1/1200)*T[0])^.5/(p*(2000-T[1]))+0.27e-1*(5/3-(1/2400)*T[2]+(1/2400)*T[0])^.5/(p*(4000-T[2]))+0.28e-1*(5/3-(1/3600)*T[3]+(1/3600)*T[0])^.5/(p*(6000-T[3]))+0.3e-1*(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5/(p*(8000-T[4]))+0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5/(p*(10000-T[5]))+0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5/(p*(12500-T[6]))+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5/(p*(15000-T[7]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5/(p*(17500-T[8]))+0.82e-1/(p*(20000-T[9]))+0.98e-1/(p*(25000-T[10]))+.164/(p*(30000-T[11]))+.145/(p*(50000-T[12]))+0.39e-1/(p*(100000-T[13]))-1 = 0}

 (2)
> 

fsolve(gleichungssystem2)

fsolve({h[0] = .472-0.33e-1*(5/3-(1/1200)*T[1]+(1/1200)*T[0])^.5-0.27e-1*(5/3-(1/2400)*T[2]+(1/2400)*T[0])^.5-0.28e-1*(5/3-(1/3600)*T[3]+(1/3600)*T[0])^.5-0.3e-1*(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5-0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5-0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5-0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5-0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5, (T[1]-T[0])/(2000-T[1]+T[0]) = 121.2121212*(0.33e-1*(5/3-(1/1200)*T[1]+(1/1200)*T[0])^.5*(1-1/(p*(2000-T[1]))-.5*(T[1]-T[0])/(2000-T[1]+T[0]))+0.27e-1*(5/3-(1/2400)*T[2]+(1/2400)*T[0])^.5*(1-1/(p*(4000-T[2]))-.5*(T[2]-T[0])/(4000-T[2]+T[0]))+0.28e-1*(5/3-(1/3600)*T[3]+(1/3600)*T[0])^.5*(1-1/(p*(6000-T[3]))-.5*(T[3]-T[0])/(6000-T[3]+T[0]))+0.3e-1*(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5*(1-1/(p*(8000-T[4]))-.5*(T[4]-T[0])/(8000-T[4]+T[0]))+0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5*(1-1/(p*(10000-T[5]))-.5*(T[5]-T[0])/(10000-T[5]+T[0]))+0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5*(1-1/(p*(12500-T[6]))-.5*(T[6]-T[0])/(12500-T[6]+T[0]))+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*(1-1/(p*(15000-T[7]))-.5*(T[7]-T[0])/(15000-T[7]+T[0]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/1200)*T[1]+(1/1200)*T[0])^.5, (T[2]-T[1])/(2000-T[2]+T[1]) = 74.07407407*(0.27e-1*(5/3-(1/2400)*T[2]+(1/2400)*T[0])^.5*(1-1/(p*(4000-T[2]))-.5*(T[2]-T[0])/(4000-T[2]+T[0]))+0.28e-1*(5/3-(1/3600)*T[3]+(1/3600)*T[0])^.5*(1-1/(p*(6000-T[3]))-.5*(T[3]-T[0])/(6000-T[3]+T[0]))+0.3e-1*(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5*(1-1/(p*(8000-T[4]))-.5*(T[4]-T[0])/(8000-T[4]+T[0]))+0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5*(1-1/(p*(10000-T[5]))-.5*(T[5]-T[0])/(10000-T[5]+T[0]))+0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5*(1-1/(p*(12500-T[6]))-.5*(T[6]-T[0])/(12500-T[6]+T[0]))+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*(1-1/(p*(15000-T[7]))-.5*(T[7]-T[0])/(15000-T[7]+T[0]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/2400)*T[2]+(1/2400)*T[0])^.5, (T[3]-T[2])/(2000-T[3]+T[2]) = 47.61904761*(0.28e-1*(5/3-(1/3600)*T[3]+(1/3600)*T[0])^.5*(1-1/(p*(6000-T[3]))-.5*(T[3]-T[0])/(6000-T[3]+T[0]))+0.3e-1*(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5*(1-1/(p*(8000-T[4]))-.5*(T[4]-T[0])/(8000-T[4]+T[0]))+0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5*(1-1/(p*(10000-T[5]))-.5*(T[5]-T[0])/(10000-T[5]+T[0]))+0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5*(1-1/(p*(12500-T[6]))-.5*(T[6]-T[0])/(12500-T[6]+T[0]))+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*(1-1/(p*(15000-T[7]))-.5*(T[7]-T[0])/(15000-T[7]+T[0]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/3600)*T[3]+(1/3600)*T[0])^.5, (T[4]-T[3])/(2000-T[4]+T[3]) = 33.33333333*(0.3e-1*(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5*(1-1/(p*(8000-T[4]))-.5*(T[4]-T[0])/(8000-T[4]+T[0]))+0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5*(1-1/(p*(10000-T[5]))-.5*(T[5]-T[0])/(10000-T[5]+T[0]))+0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5*(1-1/(p*(12500-T[6]))-.5*(T[6]-T[0])/(12500-T[6]+T[0]))+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*(1-1/(p*(15000-T[7]))-.5*(T[7]-T[0])/(15000-T[7]+T[0]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5, (T[5]-T[4])/(2000-T[5]+T[4]) = 16.66666667*(0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5*(1-1/(p*(10000-T[5]))-.5*(T[5]-T[0])/(10000-T[5]+T[0]))+0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5*(1-1/(p*(12500-T[6]))-.5*(T[6]-T[0])/(12500-T[6]+T[0]))+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*(1-1/(p*(15000-T[7]))-.5*(T[7]-T[0])/(15000-T[7]+T[0]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5, (T[6]-T[5])/(2500-T[6]+T[5]) = 15.38461538*(0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5*(1-1/(p*(12500-T[6]))-.5*(T[6]-T[0])/(12500-T[6]+T[0]))+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*(1-1/(p*(15000-T[7]))-.5*(T[7]-T[0])/(15000-T[7]+T[0]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5, (T[7]-T[6])/(2500-T[7]+T[6]) = 10.25641026*(0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*(1-1/(p*(15000-T[7]))-.5*(T[7]-T[0])/(15000-T[7]+T[0]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5, (T[8]-T[7])/(2500-T[8]+T[7]) = 12.15805471*(0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*(1-1/(p*(17500-T[8]))-.5*(T[8]-T[0])/(17500-T[8]+T[0]))+.528-0.82e-1/(p*(20000-T[9]))-0.98e-1/(p*(25000-T[10]))-.164/(p*(30000-T[11]))-.145/(p*(50000-T[12]))-0.39e-1/(p*(100000-T[13])))/(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5, (T[9]-T[8])/(2500-T[9]+T[8]) = 3.219512195-.5000000000/(p*(20000-T[9]))-.5975609756/(p*(25000-T[10]))-1.000000000/(p*(30000-T[11]))-.8841463415/(p*(50000-T[12]))-.2378048781/(p*(100000-T[13])), (T[10]-T[9])/(5000-T[10]+T[9]) = 3.640816326-.8000000000/(p*(25000-T[10]))-1.338775510/(p*(30000-T[11]))-1.183673469/(p*(50000-T[12]))-.3183673469/(p*(100000-T[13])), (T[11]-T[10])/(5000-T[11]+T[10]) = 1.414634147-.6666666668/(p*(30000-T[11]))-.5894308944/(p*(50000-T[12]))-.1585365854/(p*(100000-T[13])), (T[12]-T[11])/(20000-T[12]+T[11]) = 2.030344828-1.600000000/(p*(50000-T[12]))-.4303448276/(p*(100000-T[13])), (T[13]-T[12])/(50000-T[13]+T[12]) = 2.000000000-2.000000000/(p*(100000-T[13])), h[0]*T[0]+0.33e-1*(5/3-(1/1200)*T[1]+(1/1200)*T[0])^.5*T[1]+0.27e-1*(5/3-(1/2400)*T[2]+(1/2400)*T[0])^.5*T[2]+0.28e-1*(5/3-(1/3600)*T[3]+(1/3600)*T[0])^.5*T[3]+0.3e-1*(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5*T[4]+0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5*T[5]+0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5*T[6]+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5*T[7]+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5*T[8]+0.82e-1*T[9]+0.98e-1*T[10]+.164*T[11]+.145*T[12]+0.39e-1*T[13]-5000 = 0, -h[0]/(p*T[0])+0.33e-1*(5/3-(1/1200)*T[1]+(1/1200)*T[0])^.5/(p*(2000-T[1]))+0.27e-1*(5/3-(1/2400)*T[2]+(1/2400)*T[0])^.5/(p*(4000-T[2]))+0.28e-1*(5/3-(1/3600)*T[3]+(1/3600)*T[0])^.5/(p*(6000-T[3]))+0.3e-1*(5/3-(1/4800)*T[4]+(1/4800)*T[0])^.5/(p*(8000-T[4]))+0.48e-1*(5/3-(1/6000)*T[5]+(1/6000)*T[0])^.5/(p*(10000-T[5]))+0.52e-1*(5/3-(1/7500)*T[6]+(1/7500)*T[0])^.5/(p*(12500-T[6]))+0.65e-1*(5/3-(1/9000)*T[7]+(1/9000)*T[0])^.5/(p*(15000-T[7]))+0.47e-1*(5/3-(1/10500)*T[8]+(1/10500)*T[0])^.5/(p*(17500-T[8]))+0.82e-1/(p*(20000-T[9]))+0.98e-1/(p*(25000-T[10]))+.164/(p*(30000-T[11]))+.145/(p*(50000-T[12]))+0.39e-1/(p*(100000-T[13]))-1 = 0}, {p, T[0], T[1], T[2], T[3], T[4], T[5], T[6], T[7], T[8], T[9], T[10], T[11], T[12], T[13], h[0]})

 (3)
> 

DirectSearch:-SolveEquations([gleichungssystem2[]], initialpoint = [p = 5189.45766634515, T[0] = 4142.73188147455, T[1] = 5165.78474578260, T[2] = 3999.99988961048, T[3] = 3710.61111623718, T[4] = 2264.27637331942, T[5] = 10121.9460737861, T[6] = 14995.3110785673, T[7] = 19086.0219008713, T[8] = -22814.6070170229, T[9] = -20878.7750143556, T[10] = -16909.4648547424, T[11] = 1221.50041381499, T[12] = 4256.97724538737, T[13] = 40072.5597923931, h[0] = 1.64179118970150], evaluationlimit = 10000000)

Warning, complex or non-numeric value encountered; trying to find a feasible point

  

[37.2486614376350, Vector(16, {(1) = HFloat(1.3482006360188898), (2) = HFloat(-0.08296406227007092), (3) = HFloat(0.20777245002313566), (4) = HFloat(-0.18083286644256513), (5) = HFloat(-0.3291429319265933), (6) = HFloat(-0.19652741815664654), (7) = HFloat(0.05623163827525346), (8) = HFloat(0.11857809215140325), (9) = HFloat(-4.809358999641587), (10) = HFloat(0.0900115006818834), (11) = HFloat(-0.6565012592844091), (12) = HFloat(-2.788330148560681), (13) = HFloat(-1.7199202319393911), (14) = HFloat(0.13296887464081936), (15) = HFloat(-3.2455308974022046e-4), (16) = HFloat(-0.9297640904924019)}), [p = 5143.36074091882, T[0] = 4256.31433756824, T[1] = 5181.68565045029, T[2] = 3999.99990046679, T[3] = 4825.26702657791, T[4] = 2317.85576625528, T[5] = 10121.6396223712, T[6] = 14925.5769526497, T[7] = 19199.4262160823, T[8] = -22661.7045370251, T[9] = -20741.8150518266, T[10] = -16996.7358952387, T[11] = 1383.12333975187, T[12] = 6120.89444869930, T[13] = 40161.5912751428, h[0] = 1.50808340875895], 2410254]

 (4)
> 

ans := %

ans := [37.2486614376350, Vector(16, {(1) = HFloat(1.3482006360188898), (2) = HFloat(-0.08296406227007092), (3) = HFloat(0.20777245002313566), (4) = HFloat(-0.18083286644256513), (5) = HFloat(-0.3291429319265933), (6) = HFloat(-0.19652741815664654), (7) = HFloat(0.05623163827525346), (8) = HFloat(0.11857809215140325), (9) = HFloat(-4.809358999641587), (10) = HFloat(0.0900115006818834), (11) = HFloat(-0.6565012592844091), (12) = HFloat(-2.788330148560681), (13) = HFloat(-1.7199202319393911), (14) = HFloat(0.13296887464081936), (15) = HFloat(-3.2455308974022046e-4), (16) = HFloat(-0.9297640904924019)}), [p = 5143.36074091882, T[0] = 4256.31433756824, T[1] = 5181.68565045029, T[2] = 3999.99990046679, T[3] = 4825.26702657791, T[4] = 2317.85576625528, T[5] = 10121.6396223712, T[6] = 14925.5769526497, T[7] = 19199.4262160823, T[8] = -22661.7045370251, T[9] = -20741.8150518266, T[10] = -16996.7358952387, T[11] = 1383.12333975187, T[12] = 6120.89444869930, T[13] = 40161.5912751428, h[0] = 1.50808340875895], 2410254]

 (5)
> 

``
> 

``
> 

``


 

Download solveSys.mw

I admire...

tomleslie 13876
January 20 2019
0 0

yoor confidence about sinewave fitting, since multiple results are always possible unless a reasonable set of initial values or constraints are given.

Since the model is non-linear in (some of) the parameters, the bassic Statistics:-Fit() will not return a 'standard error'. However you can get the residuals, residualsquaresum etc,

See the attached


 

> 

  restart;
> 

#
# Read data: OP will need to change file path to something
# appropriate for his/her installation
#
  dat:= ImportMatrix( "C:/Users/TomLeslie/Desktop/Data.txt",
                      source=delimited,
                      delimiter=" "
                    ):
  p1:= plots:-pointplot(dat):
> 

#
# How not to do it!
#
# Perform a fit with no initial values given for any of
# the parameters
#
# This fit is RUBBISH: the residual mean square is
# comparable with amplitude in the original data
#
  f1:= Statistics:-Fit( A*sin(B*x+C)+D,
                        dat,
                        x,
                        output=[ leastsquaresfunction,
                                 residuals,
                                 residualmeansquare
                               ]
                     );
#
# Plot original data and the fit
#
  pf1:= plot( f1[1],
              x=dat[1,1]..dat[-1,1],
              color=red
            ):
  plots:-display([p1,pf1]);
#
# Plot residual errors
#
  plot( dat[..,1], f1[2], style=point, symbolsize=16);
 

f1 := [-0.896595324468726e-2*sin(.154794517922427*x+4.89292796410930)+.904162514221887, Vector[row](50, {(1) = .17355276786676854, (2) = .1388166842250128, (3) = 0.8547628145823583e-1, (4) = 0.20194022857735083e-1, (5) = -0.47866600382178826e-1, (6) = -.11148529242697358, (7) = -.15706387978344094, (8) = -.18403440545956284, (9) = -.18296710336205735, (10) = -.15469170737461047, (11) = -.10669847494334461, (12) = -0.4316691564404851e-1, (13) = 0.2480845509529317e-1, (14) = 0.894499476886782e-1, (15) = .14242858649223655, (16) = .17235876226563118, (17) = .17840932949246657, (18) = .16030084502194675, (19) = .11942205799728856, (20) = 0.5771773311369932e-1, (21) = -0.7044155674130859e-2, (22) = -0.7624861213931877e-1, (23) = -.13296280450065434, (24) = -.17246499600800158, (25) = -.18613865774705696, (26) = -.17427389770373436, (27) = -.13631130757376364, (28) = -0.8141725561453717e-1, (29) = -0.12915625409010367e-1, (30) = 0.5169628916013913e-1, (31) = .11353008183496782, (32) = .1553713549868101, (33) = .17555260820448415, (34) = .17157620311115251, (35) = .14316325293245458, (36) = 0.917043175143244e-1, (37) = 0.27747969737985212e-1, (38) = -0.3954095122220003e-1, (39) = -.10267152033372928, (40) = -.15163081824728475, (41) = -.17894092059178912, (42) = -.18182289983143507, (43) = -.15999643715013467, (44) = -.1159505637063808, (45) = -0.5552396816093863e-1, (46) = 0.14060730077227634e-1, (47) = 0.7809130223568705e-1, (48) = .13045889284123446, (49) = .1636663075490522, (50) = .17355024015868326}, datatype = float[8]), 0.1805624272e-1]

  

 

 
> 

#
# Perform a fit with initial guesses for the parameters
#
# Note that different initial guesses may lead to different
# (equally good) fits. For example
#
# 1. Shifting the phase (parameter C) by 2*Pi
# 2. Shifting the phase (parameter C) by Pi and simultaneously
#    flipping the sign of the amplitude (parameter A)
#
  f2:= Statistics:-Fit( A*sin(B*x+C)+D,
                        dat,
                        x,
                        initialvalues=[A=0.2, B=4, C=-1.8, D=0.9],
                        output=[ leastsquaresfunction,
                                 residuals,
                                 residualmeansquare
                               ]
                      );
#
# Plot original data and fit
#
  pf2:= plot( f2[1],
              x=dat[1,1]..dat[-1,1],
              color=red
            ):
  plots:-display([p1,pf2]);
#
# Plot residual errors
#
  plot( dat[..,1], f2[2], style=point, symbolsize=16);

f2 := [.181832870288035*sin(3.79156978183654*x-1.64292522411902)+.914540108643829, Vector[row](50, {(1) = 0.18064757013760424e-2, (2) = -0.25030056814112633e-3, (3) = -0.5037066035783466e-3, (4) = -0.6254668019197718e-3, (5) = 0.19911776577297147e-3, (6) = -0.12638325211224988e-2, (7) = -0.4570764810623018e-3, (8) = -0.31113765386350245e-2, (9) = -0.3107631003915534e-2, (10) = -0.1112769449802098e-2, (11) = -0.10263760860307336e-2, (12) = 0.3088329869911366e-3, (13) = -0.3563127313199743e-5, (14) = 0.27292506106602055e-3, (15) = 0.18384613754475199e-2, (16) = 0.21894416995771326e-3, (17) = -0.4020317795802031e-3, (18) = 0.453289568317361e-3, (19) = 0.18437545375006303e-2, (20) = -0.7680093687201239e-3, (21) = 0.2895184817535368e-2, (22) = 0.7323072474725301e-3, (23) = 0.31674141317883375e-3, (24) = -0.14580436457982682e-2, (25) = -0.14961492330960446e-2, (26) = -0.16961367397563265e-2, (27) = 0.5710821220028528e-4, (28) = -0.5590114140607838e-4, (29) = 0.2635388301808095e-2, (30) = -0.3047162673959214e-3, (31) = 0.146270374773394e-2, (32) = -0.3394668268715817e-3, (33) = -0.13696499061508494e-2, (34) = -0.8811062563960359e-3, (35) = 0.18362565507579198e-3, (36) = -0.7334450159202088e-3, (37) = -0.5118503367194105e-3, (38) = 0.9561667306395849e-3, (39) = 0.14275743036642474e-2, (40) = 0.185211020547138e-2, (41) = 0.2392956429152271e-2, (42) = 0.2265009682278052e-2, (43) = 0.11776565963290686e-2, (44) = 0.34977566714822217e-3, (45) = 0.2515166751094755e-3, (46) = 0.17931571925111633e-2, (47) = 0.3175816187949154e-3, (48) = -0.10148050425294874e-2, (49) = -0.25387481550006585e-2, (50) = -0.29762067955680926e-2}, datatype = float[8]), 2.237096241*10^(-6)]

  

 

 
> 

 


 

Download DataFit.mw

First 113 114 115 116 117 118 119 Last Page 115 of 207