## Error when specifying the optimality tolerance opt...

Hi everyone,

I have a strange problem with the LSSolve function. I get an error when specifying an optimality tolerance as an option for LSSolve function : "unexpected parameters: optimalitytolerance = 1/1000" (see attached). When I remove the option, the LSSolve function works and returns something. Why ?

This is a problem for me because in my program I automatically generate lots of equations and I need to solve them using the same parameters.  "list1" is an example of a list that makes LSSolve to return an error and therefore it makes my program stop.
Is it possible, when LSSolve returns an Error, to re-run the function without the optimalitytolerance option ?

 (1)

 >

worksheet_help.mw

Lilian

## Weird results from LSSolve when adding constraints...

Hi everyone,

I am desperatly trying to find a reason to those weird results I get using LSSolve. It could really help me to understand, maybe I am using the function the wrong way.
I have a system of equations which is overdetermined that I wrote using an electrical simulation and kirchoff's laws.
I am trying to resolve it using the LSSolve function. Here is the list of residuals :

`list := [-0.444299277411586e-2+(270.100000000000-Phi12_18)*D18, -.264819908561346+(627.030000000000-Phi23_18)*D18, .191242220011840+(-259.080000000000-Phi34_18)*D18, 0.269723795794403e-1+(-40.5060000000000-Phi45_18)*D18, 0.674200455699644e-2+(-10.1270000000000-Phi56_18)*D18, .109534122562258+(-197.290000000000-Phi67_18)*D18, 0.481462872723211e-3+(-2.41420000000000-Phi78_18)*D18, -0.346014532189641e-4+(-2.53290000000000-Phi89_18)*D18, -0.402474969346295e-4+(-2.94150000000000-Phi910_18)*D18, -0.632005430249463e-3+(-8.57100000000000-Phi1011_18)*D18, -0.105749265697549e-1+(-37.6580000000000-Phi1112_18)*D18, -0.116305497595306e-1+(-55.3250000000000-Phi1213_18)*D18, -0.581547498854927e-3+(-2.76630000000000-Phi1314_18)*D18, -0.371408130367776e-2+(-22.0900000000000-Phi1415_18)*D18, -0.886173700610320e-2+(-56.4810000000000-Phi1516_18)*D18, -0.478846208996643e-1+(262.447651185421-Phi12_18)*D29+(262.447651185421-Phi12_24)*D36, .348429199898355+(62.3165310883292-Phi23_18)*D29+(62.3165310883292-Phi23_24)*D36, .237294781239637+(41.8563477700905-Phi34_18)*D29+(41.8563477700905-Phi34_24)*D36, 0.356987380524040e-1+(6.12136413036823-Phi45_18)*D29+(6.12136413036823-Phi45_24)*D36, 0.892515544035472e-2+(1.53042068810978-Phi56_18)*D29+(1.53042068810978-Phi56_24)*D36, .163733792213247+(26.7554245920538-Phi67_18)*D29+(26.7554245920538-Phi67_24)*D36, 0.917897899527287e-3+(-0.110562085900856e-3-Phi78_18)*D29+(-0.110562085900856e-3-Phi78_24)*D36, 0.242480164562623e-4+(-.283316330467957-Phi89_18)*D29+(-.283316330467957-Phi89_24)*D36, 0.281967728090880e-4+(-.329007391842407-Phi910_18)*D29+(-.329007391842407-Phi910_24)*D36, -0.812318100863302e-3+(-1.22850243118112-Phi1011_18)*D29+(-1.22850243118112-Phi1011_24)*D36, -0.174002698946928e-1+(-9.57006175329410-Phi1112_18)*D29+(-9.57006175329410-Phi1112_24)*D36, -.125540933056649+(-44.2197489328973-Phi1213_18)*D29+(-44.2197489328973-Phi1213_24)*D36, -0.627722694977691e-2+(-2.21106159188713-Phi1314_18)*D29+(-2.21106159188713-Phi1314_24)*D36, -0.739424545575381e-1+(-24.8403831529913-Phi1415_18)*D29+(-24.8403831529913-Phi1415_24)*D36, -.203976357415920+(-68.0132712014090-Phi1516_18)*D29+(-68.0132712014090-Phi1516_24)*D36, 0.196522429267177e-1+(197.940000000000-Phi12_24)*D27, 0.368371276889244e-2+(57.8900000000000-Phi23_24)*D27, 0.144256702539785e-2+(48.4450000000000-Phi34_24)*D27, -0.115630146715321e-3+(10.-Phi45_24)*D27, -0.283028527731083e-4+(2.50010000000000-Phi56_24)*D27, -0.300476205822746e-2+(66.2640000000000-Phi67_24)*D27, -0.653509876948917e-3+(2.69040000000000-Phi78_24)*D27, -0.126753046978926e-2+(4.44790000000000-Phi89_24)*D27, -0.147212636486122e-2+(5.16530000000000-Phi910_24)*D27, -0.484316181019253e-2+(16.6000000000000-Phi1011_24)*D27, -0.298854531528585e-1+(96.8770000000000-Phi1112_24)*D27, -.120604432493978+(315.410000000000-Phi1213_24)*D27, -0.603334119632106e-2+(15.7700000000000-Phi1314_24)*D27, -0.664471982996522e-1+(167.170000000000-Phi1415_24)*D27, 0.786913003105101e-1+(-326.760000000000-Phi1516_24)*D27] `

I know that all D values must be positive. When resolving the system without any constraints (D >= 0), i get the values I expected (knowing the input I used in the simulation), with a really low error :

`result := LSSolve(list);`

`[1.82130325886306*10^(-8), [D10 = 0.200009334740825e-2, D11 = 0.666620509302803e-3, D14 = 0.222215208246154e-2, D15 = 0.128202791383597e-2, D17 = 0.499886140344411e-2, D19 = 0.302925526676043e-3, D2 = 0.100002349341980e-2, D20 = 0.142849446596938e-3, D22 = 0.111121127122156e-1, D23 = 0.222228054119820e-2, D25 = 0.714293621502836e-3, D26 = 0.833326349912537e-3, D28 = 0.217396531719902e-3, D3 = 0.400217567900069e-3, D6 = 0.166878862202449e-3, D7 = 0.999969828547956e-2, Phi1011_17 = -1.22850243118112, Phi1011_19 = -20.5335193736012, Phi1011_21 = -104.090964313150, Phi1011_23 = 19.2144499395683, Phi1112_17 = -9.57006175329410, Phi1112_19 = -81.6848630234903, Phi1112_21 = -242.149849175388, Phi1112_23 = 109.001351349915, Phi1213_17 = -44.2197489328973, Phi1213_19 = -92.8267195929548, Phi1213_21 = -204.444165890808, Phi1213_23 = -61.4447612788985, Phi12_17 = 262.447651185421, Phi12_19 = 262.149192406679, Phi12_21 = 256.248405276737, Phi12_23 = 246.521172863223, Phi1314_17 = -2.21106159188713, Phi1314_19 = -4.64255435896474, Phi1314_21 = -10.2212158757032, Phi1314_23 = -3.07798400495386, Phi1415_17 = -24.8403831529913, Phi1415_19 = -30.5944507718603, Phi1415_21 = -45.5847025259923, Phi1415_23 = -77.3297680041818, Phi1516_17 = -68.0132712014090, Phi1516_19 = -74.2023324471993, Phi1516_21 = -95.1952296374558, Phi1516_23 = -132.328467080565, Phi23_17 = 62.3165310883292, Phi23_19 = 200.804225452845, Phi23_21 = 130.018791598707, Phi23_23 = 73.7043262431720, Phi34_17 = 41.8563477700905, Phi34_19 = 343.409987932231, Phi34_21 = 159.593996060841, Phi34_23 = 62.6564757702407, Phi45_17 = 6.12136413036823, Phi45_19 = 12.3839171939746, Phi45_21 = 46.0005281797016, Phi45_23 = 13.1665796516893, Phi56_17 = 1.53042068810978, Phi56_19 = 3.16614687399595, Phi56_21 = 11.4998114891963, Phi56_23 = 3.29093394692614, Phi67_17 = 26.7554245920538, Phi67_19 = -244.288977944524, Phi67_21 = 376.351493538080, Phi67_23 = 88.4830465193635, Phi78_17 = -0.110562085900856e-3, Phi78_19 = -6.28061380389266, Phi78_21 = 43.7035845962372, Phi78_23 = 3.35123473697264, Phi89_17 = -.283316330467957, Phi89_19 = -6.18811507913178, Phi89_21 = -13.9258224376815, Phi89_23 = 5.20325572546379, Phi910_17 = -.329007391842407, Phi910_19 = -7.18580970783931, Phi910_21 = -16.1669897128450, Phi910_23 = 6.04291224185087]] `

When adding the constraints that D should be positive (and that are actually positive in the previous result), I get a worse result in term of precisions :

`LSSolve(list, {D10 >= 0, D11 >= 0, D14 >= 0, D15 >= 0, D17 >= 0, D19 >= 0, D2 >= 0, D20 >= 0, D22 >= 0, D23 >= 0, D25 >= 0, D26 >= 0, D28 >= 0, D3 >= 0, D6 >= 0, D7 >= 0});`

`[0.667302976414869964e-1, [D10 = 0.240199442379079e-2, D11 = 0.666577572133538e-3, D14 = 0.222218786790062e-2, D15 = 0.128192441757651e-2, D17 = 0.278678889056743e-2, D19 = 0.200473317719685e-3, D2 = 0.109938538155804e-2, D20 = 0.840721762649974e-4, D22 = 0.685770482726534e-3, D23 = -1.387530857*10^(-312), D25 = 0.714397733627028e-3, D26 = 0.833201232339238e-3, D28 = 0.204319731851617e-3, D3 = 0.419994015872111e-3, D6 = 0.191996909862889e-3, D7 = 0.103884505319047e-1, Phi1011_17 = -.709707335593168, Phi1011_19 = -15.7863975896827, Phi1011_21 = -151.171843708558, Phi1011_23 = 19.2211409030343, Phi1112_17 = -8.90604676283968, Phi1112_19 = -75.8627539382983, Phi1112_21 = -311.423930967299, Phi1112_23 = 109.002880650927, Phi1213_17 = -54.9212365194647, Phi1213_19 = -89.9790565093006, Phi1213_21 = -250.971671756001, Phi1213_23 = -61.5160003335629, Phi12_17 = 251.480872515883, Phi12_19 = 255.977573006508, Phi12_21 = 254.397100891354, Phi12_23 = 246.524672366158, Phi1314_17 = -2.74614386328796, Phi1314_19 = -4.48401822538664, Phi1314_21 = -12.5381572344771, Phi1314_23 = -3.08154491280567, Phi1415_17 = -31.8947514252141, Phi1415_19 = -30.8090400512349, Phi1415_21 = -51.0499196769535, Phi1415_23 = -77.3268969229600, Phi1516_17 = -87.7947790488482, Phi1516_19 = -75.5403005246575, Phi1516_21 = -101.763771364478, Phi1516_23 = -132.314524393221, Phi23_17 = 94.0093590848714, Phi23_19 = 86.4429757025976, Phi23_21 = 108.554765004168, Phi23_23 = 73.7072279431268, Phi34_17 = 87.6938924370977, Phi34_19 = 82.3922347753764, Phi34_21 = 88.5582078636840, Phi34_23 = 62.6604078051191, Phi45_17 = 13.1910198060107, Phi45_19 = 69.3008595136787, Phi45_21 = 15.5530983566712, Phi45_23 = 13.1677681559684, Phi56_17 = 3.29792072169498, Phi56_19 = 17.4003349272078, Phi56_21 = 3.88187632917493, Phi56_23 = 3.29123115133383, Phi67_17 = 60.0045707036166, Phi67_19 = 54.3070868626015, Phi67_21 = 87.9421288802858, Phi67_23 = 88.4929920125095, Phi78_17 = .279952186311827, Phi78_19 = -3.50632712699693, Phi78_21 = -20.3872167203319, Phi78_23 = 3.35213748642018, Phi89_17 = -0.991299169828910e-1, Phi89_19 = -4.44636683843093, Phi89_21 = 297.888811926331, Phi89_23 = 5.20500671661437, Phi910_17 = -.115101700555720, Phi910_19 = -5.16603761776826, Phi910_21 = 346.033351291632, Phi910_23 = 6.04494564310825]] `

I also get the warning "limiting number of major iterations has been reached".
Can someone explain me?

It may not seem important at first sight, but sometimes when using my program I get wrong values and a negative D, which is not possible. Therefore I try to add a positive constraint, but the LLSolve function doesn't return anything except the error "no improved point could be found", which is weird because when I manually substitute the value I consider correct, i get a really low error. I can show you the related list of equations if you are interested...

Lilian

## fitting data on ode system...

hello agian i have the following problem. i need to fit data using the following model:

ode_sub := diff(S(t), t) = -k1*S(t)-S(t)/T1_s;

ode_P1 := diff(P1(t), t) = 2*k1*S(t)-k2*(P1(t)-P2(t)/keq)-P1(t)/T1_p1;

ode_P2 := diff(P2(t), t) = -k2*(-keq*P1(t)+P2(t))/keq-k4*P2(t)-P2(t)/T1_p2;

ode_P2e := diff(P2_e(t), t) = k4*P2(t)-P2_e(t)/T1_p2_e;

ode_system := ode_sub, ode_P1, ode_P2, ode_P2e

known paramters: s0 := 10000; k2 := 1000; T1_s := 14; T1_p2_e := 35; T1_p2 := T1_p1

initial conditions: init := S(0) = s0, P1(0) = 0, P2(0) = 0, P2_e(0) = 0

using the following data the fitting is fine:

T := [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 72, 74, 76, 78, 80, 82, 84, 86, 88, 90, 92, 94, 96, 98, 100]

data_s := [10000.00001377746, 7880.718836217512, 6210.572917625112, 4894.377814704496, 3857.121616806618, 3039.689036293312, 2395.493468824973, 1887.821030424934, 1487.738721378254, 1172.445015238064, 923.970970533911, 728.1555148262845, 573.8389058920359, 452.2263410725434, 356.3868133709972, 280.8584641247961, 221.3366863178145, 174.4291648589732, 137.4627967029932, 108.3305246342751, 85.37230370803576, 67.2794974831867, 53.0210755540487, 41.78436556408739, 32.92915589156605, 25.95049906181449, 20.45092590987601, 16.11678944747619, 12.70115104646253, 10.009439492356, 7.888058666357939, 6.216464754698855, 4.899036441885205, 3.860793948052506, 3.042584031379331, 2.397778731385364, 1.889601342133927, 1.489145259940784, 1.173591417634974, .9248694992255094, .7288443588090404, .5743879613705721, .4526024635132348, .3567687570029392, .2810508404011898, .2215650452813882, .174603181767829, .1376243372528356, .1084232889842753, 0.8542884707952822e-1, 0.6738282660463157e-1]

data_p1 := [0.1194335334401124e-4, 244.8746046039949, 374.8721199398692, 430.5392805383767, 439.6598364813143, 421.0353424914179, 387.1842556288343, 346.2646897222593, 303.4377508746471, 261.8283447091155, 223.1996547160051, 188.4213144493491, 157.7924449350029, 131.2622073983344, 108.5771635112278, 89.37951009190863, 73.26979150957087, 59.84578653950572, 48.72563358658898, 39.56010490461378, 32.03855466968536, 25.88922670933322, 20.87834763772145, 16.80708274458702, 13.50774122768974, 10.84014148654258, 8.687656394505874, 6.954093898245485, 5.560224107929433, 4.441209458524726, 3.544128529596104, 2.825755811619965, 2.251247757181308, 1.792233305651086, 1.425861347838012, 1.133566009019768, .90081361320016, .7153496336919163, .5678861241754847, .4505952916932289, .3572989037753538, .2832489239941939, .2244289868248577, .1778450590752305, .1408633578784151, .1114667192753896, 0.8826814044702111e-1, 0.6979315954603076e-1, 0.5526502783606788e-1, 0.4370354298880999e-1, 0.3456334307662573e-1]

data_p2_p2e := [-0.1821397630630296e-4, 1000.40572909871, 1568.064904416198, 1848.900129881268, 1944.147939710225, 1923.352973299286, 1833.705314342611, 1706.726235937363, 1563.036042115902, 1415.741121363331, 1272.825952816517, 1138.833091575137, 1016.03557293539, 905.2470623752856, 806.3754051707843, 718.7979384094563, 641.6091822001032, 573.7825966275556, 514.2718125966452, 462.0710416682647, 416.2491499591012, 375.9656328260581, 340.4748518743264, 309.124348579787, 281.3484612079911, 256.6599441255977, 234.6416876403397, 214.9378461256197, 197.2453333305823, 181.3059798685686, 166.90059218783, 153.8422174795751, 141.9717783871606, 131.1529759058513, 121.2692959115991, 112.2199678247825, 103.9187616370328, 96.29007054510998, 89.26877137303566, 82.79743967408479, 76.82586273439793, 71.30944943617081, 66.2085909515396, 61.48838150744805, 57.11714763242225, 53.06666224006544, 49.31130738219119, 45.82807853990728, 42.59597194910467, 39.59575450632008, 36.81013335261527]

the fitting is done the following way:

P1fu, P2fu, P2e_fu, Sfu := op(subs(res, [P1(t), P2(t), P2_e(t), S(t)]))

making residuals:

Q := proc (T1_p1, k1, keq, k4) local P1v, P2v, P2e_v, Sv, resid; option remember;

res(parameters = [T1_p1, k1, keq, k4]);

try P1v := `~`[P1fu](T); P2v := `~`[P2fu](T); P2e_v := `~`[P2e_fu](T); Sv := `~`[Sfu](T); resid := [P1v-data_p1, P2v+P2e_v-data_p2_p2e, Sv-data_s]; return [seq(seq(resid[i][j], i = 1 .. 3), j = 1 .. nops(T))] catch: return [1000000\$3*nops(T)] end try end proc;
q := [seq(subs(_nn = n, proc (T1_p1, k1, keq, k4) Q(args)[_nn] end proc), n = 1 .. 3*nops(T))];

finding inital point for the LSsolve:

L := fsolve(q[2 .. 5], [10, 0.2e-1, 4, 4])

fitting the data with the intial point:

solLS := Optimization:-LSSolve(q, initialpoint = L)

this is all good however, when i used the following data it did not turn out so well (using the same approch as above):

data_s := [96304.74567, 77385.03700, 62621.83067, 51239.94333, 42663.82367, 35084.74100, 28480.28367, 23066.01467, 18774.73700, 15179.13700, 12278.50767, 9937.652000, 8046.848333, 6521.242000, 5287.811667, 4277.779000, 3466.518333, 2835.467000, 2297.796333, 1861.249667, 1529.654000, 1235.353000, 999.6626667, 826.2343333, 667.9480000, 559.9230000, 449.2790000, 376.4860000, 289.1203333, 236.1483333]

data_p1 := [0.86e-1, 3.904, 26.975, 31.719, 41.067, 46.779, 52.115, 43.101, 44.344, 41.094, 36.523, 27.742, 26.543, 28.062, 22.178, 21.303, 14.951, 17.871, 11.422, 12.051, 9.232, 6.817, 6.1, .717, 1.215, 6.146, .772, .375, 2.595, .518]

data_p2_p2e := [-3.024, 22.238, 61.731, 103.816, 132.695, 159.069, 167.302, 160.188, 158.398, 152.943, 146.745, 135.22, 132.145, 120.413, 107.864, 95.339, 90.775, 81.828, 71.065, 70.475, 62.872, 49.955, 40.858, 42.938, 41.311, 35.583, 31.573, 29.841, 29.558, 21.762]

the known parameters is the case are:

s0 := 96304.74567; k2 := 10^5; T1_s := 14; T1_p2_e := 35; T1_p2 := T1_p1

additionally the following fuction affects the solution:

K:=t->cos((1/180)*beta*Pi)^(t/Tr)

i included this by doing this:

P1fu_K := proc (t) options operator, arrow; P1fu(t)*K(t) end proc;

P2fu_K := proc (t) options operator, arrow; P2fu(t)*K(t) end proc;

P2e_fu_K := proc (t) options operator, arrow;

P2e_fu(t)*K(t) end proc;

Sfu_K := proc (t) options operator, arrow; Sfu(t)*K(t) end proc

resulting in the following residuals:

Q := proc (T1_p1, k1, keq, k4) local P1v, P2v, P2e_v, Sv, resid; option remember;

res(parameters = [T1_p1, k1, keq, k4]);

try P1v := `~`[P1fu_K](T); P2v := `~`[P2fu_K](T); P2e_v := `~`[P2e_fu_K](T); Sv := `~`[Sfu_K](T); resid := [P1v-data_p1, P2v+P2e_v-data_p2_p2e, Sv-data_s]; return [seq(seq(resid[i][j], i = 1 .. 3), j = 1 .. nops(T))] catch: return [1000000\$3*nops(T)] end try end proc;
q := [seq(subs(_nn = n, proc (T1_p1, k1, keq, k4) Q(args)[_nn] end proc), n = 1 .. 3*nops(T))];

i think my problem is that the inital point in this case is not known. all i know is that all the fitted parameters should be positive and that k1<1, k4>10, keq>1 and T1_p1>100 (more i do not know) - is there a way to determin the inital point without guessing?

i also know the results, which should be close to these values:

k1=0.000438, k4=0.0385, keq=2.7385 and T1_p1=36.8 the output fit should look something like this

where the red curve is (PP2(t)+PP2_e(t))*K(t)

the blue cure is PP1(t)*K(t)

anyone able to help - i've tried for 2 days now. it might be that  ode_P1 := diff(P1(t), t) = 2*k1*S(t)-k2*(P1(t)-P2(t)/keq)-P1(t)/T1_p1 should be changed into ode_P1 := diff(P1(t), t) = k1*S(t)-k2*(P1(t)-P2(t)/keq)-P1(t)/T1_p1;

i've tried this but it didnt seem to do much

anyone able to help?:)

NB stiff=true can be used within the dsolve to speed up the process if needed:)

## LSSolve on a large ODE system...

hello i have the following set of ode's:

ode_sub := diff(S(t), t) = -k1*S(t)-S(t)/T1_s;
ode_P1 := diff(P1(t), t) = k1*S(t)-k2*(P1(t)-P2(t)/keq)-P1(t)/T1_p1;
ode_P2 := diff(P2(t), t) = -k2*(-keq*P1(t)+P2(t))/keq-k4*P2(t)-P2(t)/T1_p2;
ode_P2e := diff(P2_e(t), t) = k4*P2(t)-P2_e(t)/T1_p2_e;

ode_system := ode_sub, ode_P1, ode_P2, ode_P2e;

with these parameters:
s0 := 10000;
k2 := 1000; T1_s := 14; T1_p2_e := 35; T1_p2 := T1_p1;

i want to find the unkown parameters : T1_p1, k1, keq and k4

init:=S(0)=s0,P1(0)=0,P2(0)=0,P2_e(0)=0

dsolve({ode_system,init})

sol := combine(expand(%));
PS := subs(sol, [S(t), P1(t), P2(t), P2_e(t)]);

P1fu := unapply(PS[2],t);
Sfu := unapply(PS[1],t);
P2fu := unapply(PS[3],t);
P2e_fu := unapply(PS[4],t);
P2_total := unapply(P2fu+P2e_fu, t);

the following data is given:

T:=<0,2,4,6,8>

S:=<9999.99913146527,8328.870587730016,6937.009129218748,5777.745632133724,4812.209983843559>

P1:=<0.0,67.86790056712294,114.88787098501874,145.95438088662502,164.85650644237887>

P2_P2e:=<0.0,271.68492651947497,461.9130396605823,589.3710176125417,668.9967533337124> # data from P2(t)+P2_e(t)

making the rediduals:

RP1 := convert(P1-P1fu~(T), list);
RS := convert(S-Sfu~(T), list);
RP2_P2e := convert(P2_P2_e-P2_total~(T), list);

RPs := [op(RS), op(RP2_P2_e), op(RP1)]

res := Optimization:-LSSolve(RPs, k1 = 0 .. 1, keq = 0 .. 10, k4 = 0 .. 1, T1_p1 = 0 .. 100)

i dont know wheter or not the last step work to get the parameters becuase it takes to long to compute. is there a smarter way to obtain the parameters of the ode's? a numeric approch ?

i tried with dsolve({ode_sysytem,init},numeric,'parameters'=[k1,keq,k4,T1_p1]) however it doesnt seem to get my anywhere since i need to know the parameters to use this (i think)

hope someone can help:)

## A Problem with LSSOLVE for ODE...

Following my previous question

http://www.mapleprimes.com/questions/200627-Lssolve-Midpoint

I wrote the following code

restart:
Phiavg:=0.06;
lambda:=0.05;
Ha:=0;
NBT:=0.5;
Nr:=500;
#N[bt]:=cc*NBT+(1-cc)*4; ## cc between 0 and 1
N[bt]:=cc*NBT+(1-cc^2)*0.75;

0.06
0.05
0
0.5
500
2
0.5 cc + 0.75 - 0.75 cc
eq1:=diff(u(eta),eta,eta)+1/(mu(eta)/mu1[w])*(sigma-Nr*(phi(eta)-phi1[w])-(1-phi(eta))*T(eta)-Ha^2*u(eta))+((1/mu(eta)*(mu_phi*diff(phi(eta),eta)))*diff(u(eta),eta));
eq2:=diff(T(eta),eta)-1/(k(eta)/k1[w]);
eq3:=diff(phi(eta),eta)-phi(eta)/(N[bt]*(1-gama1*T(eta))^2)*diff(T(eta),eta);
/  d   /  d         \\      1
|----- |----- u(eta)|| + ------- (mu1[w] (sigma - 500 phi(eta)
\ deta \ deta       //   mu(eta)

+ 500 phi1[w] - (1 - phi(eta)) T(eta)))

/  d           \ /  d         \
mu_phi |----- phi(eta)| |----- u(eta)|
\ deta         / \ deta       /
+ --------------------------------------
mu(eta)
/  d         \   k1[w]
|----- T(eta)| - ------
\ deta       /   k(eta)
/  d         \
phi(eta) |----- T(eta)|
/  d           \                       \ deta       /
|----- phi(eta)| - ----------------------------------------------
\ deta         /   /                       2\                   2
\0.5 cc + 0.75 - 0.75 cc / (1 - gama1 T(eta))
mu:=unapply(mu1[bf]*(1+a[mu1]*phi(eta)+b[mu1]*phi(eta)^2),eta):
k:=unapply(k1[bf]*(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2),eta):
rhop:=3880:
rhobf:=998.2:
cp:=773:
cbf:=4182:
rho:=unapply(  phi(eta)*rhop+(1-phi(eta))*rhobf ,eta):
c:=unapply(  (phi(eta)*rhop*cp+(1-phi(eta))*rhobf*cbf )/rho(eta) ,eta):
mu_phi:=mu1[bf]*(a[mu1]+2*b[mu1]*phi(eta)):
gama1:=0.00:
a[mu1]:=39.11:
b[mu1]:=533.9:
mu1[bf]:=9.93/10000:
a[k1]:=7.47:
b[k1]:=0:
k1[bf]:=0.597:
zet:=1:
phi1[w]:=phi0:
mu1[w]:=mu(0):
k1[w]:=k(0):

eq1:=subs(phi(0)=phi0,eq1);
eq2:=subs(phi(0)=phi0,eq2);
eq3:=subs(phi(0)=phi0,eq3);
/  d   /  d         \\   //
|----- |----- u(eta)|| + \\0.0009930000000 + 0.03883623000 phi0
\ deta \ deta       //

2\
+ 0.5301627000 phi0 / (sigma - 500 phi(eta) + 500 phi0

\//
- (1 - phi(eta)) T(eta))/ \0.0009930000000

2\
+ 0.03883623000 phi(eta) + 0.5301627000 phi(eta) / +

/                                       /  d           \ /  d
|(0.03883623000 + 1.060325400 phi(eta)) |----- phi(eta)| |-----
\                                       \ deta         / \ deta

\\//
u(eta)|| \0.0009930000000 + 0.03883623000 phi(eta)
//

2\
+ 0.5301627000 phi(eta) /
/  d         \     0.597 + 4.45959 phi0
|----- T(eta)| - ------------------------
\ deta       /   0.597 + 4.45959 phi(eta)
/  d         \
1. phi(eta) |----- T(eta)|
/  d           \               \ deta       /
|----- phi(eta)| - --------------------------
\ deta         /                           2
0.5 cc + 0.75 - 0.75 cc
Q:=proc(pp2,fi0) option remember; local res,F0,F1,F2,a,INT0,INT10,B;
print(pp2,fi0);
if not type([pp2,fi0],list(numeric)) then return 'procname(_passed)' end if;
res := dsolve(subs(sigma=pp2,phi0=fi0,{eq1=0,eq2=0,eq3=0,u(1)=-lambda*D(u)(1),u(0)=lambda*D(u)(0),phi(0)=phi0,T(0)=0}), numeric,output=listprocedure,initmesh=10, continuation=cc);
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)]));
INT0:=evalf(Int((abs(F0(eta)),eta=0..1)));
INT10:=evalf(Int(abs(F0(eta))*F1(eta),eta=0..1));
a[1]:=evalf(Int(F0(eta)*(F1(eta)*rhop+(1-F1(eta))*rhobf),eta=0..1));
#a[1]:=evalf(Int((F0(eta),eta=0..1)));
a[2]:=(INT10/INT0-Phiavg)/Phiavg; #relative
[a[1],a[2]]
end proc:
Q1:=proc(pp2,fi0) Q(_passed)[1] end proc;
Q2:=proc(pp2,fi0) Q(_passed)[2] end proc;
proc(pp2, fi0)  ...  end;
proc(pp2, fi0)  ...  end;
#Q(116,0.0041);
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[130,0.01]);
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[43.55,0.39]);
tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[5.65,0.00036]);
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[12,0.003]); # khoob ba 1
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[5,0.01]);
HFloat(5.65), HFloat(3.6e-4)
HFloat(5.650000070103341), HFloat(3.6e-4)
HFloat(5.65), HFloat(3.600105456508193e-4)
HFloat(29.63242379055208), HFloat(0.0205927592420527)
HFloat(12.803902258015825), HFloat(0.006395385884750864)
HFloat(12.803902403534572), HFloat(0.006395385884750864)
HFloat(12.803902258015825), HFloat(0.00639539649402585)
HFloat(12.804004931505949), HFloat(0.0063954867657199386)
HFloat(12.804107604996073), HFloat(0.006395587646689013)
HFloat(12.80400483062498), HFloat(0.006498160255844027)
HFloat(12.803902157134855), HFloat(0.006498059374874952)
HFloat(-1.0206939292143726), HFloat(-3.32764179807047e-4)
HFloat(-1.0206939079125088), HFloat(-3.32764179807047e-4)
HFloat(-1.0206939292143726), HFloat(-3.327536344433438e-4)
HFloat(18.749500943683863), HFloat(0.01993840615828979)
HFloat(3.9953780262640484), HFloat(0.00481041471606933)
HFloat(6.166152606930136), HFloat(0.00703619658484674)
HFloat(7.3193201827812295), HFloat(0.008218585352824569)
Error, (in Optimization:-LSSolve) complex value encountered
sigma:=tempe[2](1);
tempe[2](1)
phi0:=tempe[2](2);
tempe[2](2)
with(plots):

res2 := dsolve({eq1=0,eq2=0,eq3=0,u(1)=-lambda*D(u)(1),u(0)=lambda*D(u)(0),phi(0)=phi0,T(0)=0}, numeric,output=listprocedure,continuation=cc);
Error, (in dsolve/numeric/process_input) boundary conditions specified at too many points: {0, 1, 2}, can only solve two-point boundary value problems
G0,G1,G2:=op(subs(res2,[u(eta),phi(eta),T(eta)])):
ruu:=evalf((Int(abs(G0(eta))*(G1(eta)*rhop+(1-G1(eta))*rhobf ),eta=0..zet)))/(Phiavg*rhop+(1-Phiavg)*rhobf);
phb:=evalf((Int(abs(G0(eta))*G1(eta),eta=0..1))) / evalf((Int(abs(G0(eta)),eta=0..1))) ;
TTb:=evalf(Int(abs(G0(eta))*G2(eta)*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1))/evalf(Int(abs(G0(eta))*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1));
Error, invalid input: subs received res2, which is not valid for its 1st argument
/  /1.                                        \
| |                                           |
0.0008538922115 | |    |G0(eta)| (2881.8 G1(eta) + 998.2) deta|
| |                                           |
\/0.                                          /
/1.
|
|    |G0(eta)| G1(eta) deta
|
/0.
----------------------------
/1.
|
|
|    |G0(eta)| deta
/
0.
/Int(
1                               |
------------------------------------------------------------- |
/1.                                                         |
|                                                            \
|              /             6                       6\
|    |G0(eta)| \-1.1752324 10  G1(eta) + 4.1744724 10 / deta
/
0.

/             6                       6\ , eta = 0. .. 1.)
|G0(eta)| G2(eta) \-1.1752324 10  G1(eta) + 4.1744724 10 /

\
|
|
|
/
#rhouu:=evalf((Int((G1(eta)*rhop+(1-G1(eta))*rhobf)*G0(eta),eta=0..1)));

odeplot(res2,[[eta,u(eta)/ruu],[eta,phi(eta)/phb],[eta,T(eta)/TTb]],0..1);
#odeplot(res2,[[eta,u(eta)],[eta,phi(eta)],[eta,T(eta)]],0..1);
Error, (in plots/odeplot) input is not a valid dsolve/numeric solution
rhou:=evalf((Int(abs(G0(eta))*(G1(eta)*rhop+(1-G1(eta))*rhobf ),eta=0..zet))):

Nub:=(1/G2(1))*(((1+a[k1]*abs(G1(0))+b[k1]*abs(G1(0))^2)/(1+a[k1]*Phiavg+b[k1]*Phiavg^2)));
0.6905123602 (1 + 7.47 |G1(0)|)
-------------------------------
G2(1)
(rhs(res2(0.0000000000001)[3])-rhs(res2(0)[3]))/0.0000000000001;
Error, invalid input: rhs received res2(0.1e-12)[3], which is not valid for its 1st argument, expr
sigma;
tempe[2](1)
NBT;
0.5
>

the above code has been worked for NBT=0.6 and higher, whereas as NBT decreases, the code doesnt converge easily.

How can I fix this problem?

Amir