Items tagged with lssolve

Feed

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:)

 

 

 

 


 

 

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

my idea was this:

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:)

 

 

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?

Thanks for your attention in advance

Amir

I am confused by the below results.

Why does LSSolve (given an arbitrary expected return) produce a higher risk
adjusted return than QPSolve which explicity is given an objective to maximize
risk adjusted returns? ie minimize Transpose(W).Cov.W-Transpose(W).ev
This to me seems very strange?!

Also, how do you specify in the objective function for LSSolve to maximize
risk adjusted returns? Now we have simply provide LSSolve with some user specified

Page 1 of 1