tomleslie

Some observations...

Answered: tomleslie 13729

June 10 2023

0 0

Your system, is third-order in f(), second-order in g(), and second-order in h(). Thus it can only support a total of 7 boundary/initial conditions. You have seven initial conditions (ie at x=0) and three boundary conditions (ie at x=infinity). That is three too many.
In such a case, I would suggext that you delete the three initial conditions given by , f''(0)=a[1], g'(0)=a[2], h'(0)=a[3], then include three boundary conditions f(inf)=0, g(inf)=0, h(inf)=0, where 'inf' is an "approximation of infinity" suitable for the scal of your problem, a value of 10 is not a bad start.
With the correct number of initial/boundary conditions, solve the system numerically. This will require assigning a value to 'alpha'. An example solution is shown in the attached where I have use inf=10 and alpha=Pi/3. The values of f''(0), g'(0), and h'(0) can be obtained from this solution, if required.

>

  restart;
  with(plots):
  inf:=10:
  alpha:=Pi/3:
  eq1 := diff(f(x), x, x, x)+(1/2)*cos(alpha)*x*(diff(f(x), x, x))+(1/2)*sin(alpha)*f(x)*(diff(f(x), x, x)) = 0:

  eq2 := diff(g(x), x, x)+diff(g(x), x)+(diff(g(x), x))*(diff(h(x), x))+cos(alpha)*x*(diff(g(x), x))+sin(alpha)*f(x)*g(x) = 0:

  eq3 := diff(g(x), x, x)+diff(h(x), x, x)+1/2*(cos(alpha)*x+sin(alpha)*f(x)) = 0:

  ics:=f(0)=0, D(f)(0)=1, f(inf)=0, g(0)=0, g(inf)=0, h(0)=1, h(inf)=0:

  #bcs:=f(x) , g(x), h(x) tends to 0 ad x tends to infinity

  sol:=dsolve( [eq1, eq2, eq3, ics], numeric):
  sol(0);
  odeplot( sol,
           [ [ x, f(x) ],
             [ x, g(x) ],
             [ x, h(x) ]
           ],
           x=0..inf,
           color=[ red, blue, black],
           axes=boxed
         );

[x = 0., f(x) = HFloat(0.0), diff(f(x), x) = HFloat(1.0000000000000002), diff(diff(f(x), x), x) = HFloat(-0.6408049672682483), g(x) = HFloat(0.0), diff(g(x), x) = HFloat(-1.4342282007671064e-30), h(x) = HFloat(1.0000000000000002), diff(h(x), x) = HFloat(5.5011987078123825)]

>

Download tooManyCond.mw

Running correctly in Maple 18...

Answered: tomleslie 13729

June 08 2023

0 0

see the attached.

I think you have to accept that using software releases which are (more-or-less) ten years out of date, is asking for trouble.

In future please use the big green up-arrow in the Mapleprimes toolbar to upload worksheets. A "picture" of a worksheet is about as much use as a chocolate teapot

>

  restart;
  interface(version);
  with(plots):
  with(plottools):
  sol7:= dsolve({D(N)(t) = N(t), N(0)=1},
                 numeric
               ):
#
# The scale, shifted, reversed, plot of the ODE solution
#
  display
  ( textplot([0.6,120,N(t)]),
    odeplot
    ( sol7,
      [5-t, (N(t)-rhs(sol7(0)[2]))*150/( rhs(sol7(5)[2])-rhs(sol7(0)[2]))],
      0..5,
      color = red,
      thickness = 2,
      linestyle = dot
    )
  );

>

Download oddPlots5.mw

Or...

Answered: tomleslie 13729

June 08 2023

0 0

just as weird, but more succinctly

>

  restart;
  with(plots):
  with(plottools):
  sol7:= dsolve({D(N)(t) = N(t), N(0)=1},
                 numeric
               ):
#
# The scale, shifted, reversed, plot of the ODE solution
#
  p1:=display
      ( textplot([0.6,120,N(t)]),
        odeplot
        ( sol7,
          [5-t, (N(t)-rhs(sol7(0)[2]))*150/( rhs(sol7(5)[2])-rhs(sol7(0)[2]))],
          0..5,
          color = red,
          thickness = 2,
          linestyle = dot
        )
      );

>

Download oddPlots4.mw

It gets weirder...

Answered: tomleslie 13729

June 08 2023

0 0

@jmakinde

The attached reflects the plot about x=2.5, then scales and shifts in the y-direction.

>

  restart;
  with(plots):
  with(plottools):
  sol7:= dsolve({D(N)(t) = N(t), N(0)=1},
                 numeric
               ):
#
# The "normal" plot
  p1:=display
      ( textplot([4.4,120,N(t)]),
        odeplot
        ( sol7,
          [t, N(t)],
          0..5,
          color = red,
          thickness = 2,
          linestyle = dot
        )
      );
#
# The normal plot, scaled and shifted in the y-direction
# and spoofed tickmarks
#
  f1:=transform( (x, y)-> [5-x, (y-rhs(sol7(0)[2]))*150/( rhs(sol7(5)[2])-rhs(sol7(0)[2])) ] ):
  display( f1(p1),
           tickmarks=[[seq( j=5-j, j=0..5)], default]
         );

>

Download oddPlots3.mw

I should probably have mentioned...

Answered: tomleslie 13729

June 08 2023

0 0

@AHSAN

that you can achieve the same thing without the view option just by plotting data[12..-1, [1, 4]], rather than data[2..-1, [1, 4]]

Use the view option...

Answered: tomleslie 13729

June 08 2023

0 0

@AHSAN

as shown in the attached

>

>	with(plots):

>

p1 := display
      ( plot
         ( data[2..-1, [1, 4]],
           color=blue,
           linestyle=3
         ),
         plot
         ( data[2..-1, [1, 4]],
            style=point,
            color=blue,
            symbol=asterisk,
            symbolsize=20,
            color=red
         ),
         axes=boxed,
         labels=[y, "Absolute Error"],
         labeldirections=[default, vertical],
         view=[0..1, default]
       );

>

Download basicPlot2.mw

I agree with Acer...

Answered: tomleslie 13729

June 07 2023

0 0

that this is a *weird* thing to do, although I can think of a couple of ways it can be achieved

Either

Use the original ODE
Plot "normally"
Reflect the "normal" plot about x=2.5
Spoof the tickmarks on the reflected plot

Or

Reflect the original ODE
Plot normally
Spoof the tickmarks

Both methods are shown in the attached (NB unwanted gridlines do not appear in the actual Maple worksheet - they are a *feature* of the rendering on this site)

>

  restart;
  with(plots):
  with(plottools):
  sol7:= dsolve({D(N)(t) = N(t), N(0)=1},
                 numeric
               ):
#
# The "normal" plot
#
  display
  ( textplot([4.4,120,N(t)]),
    odeplot
    ( sol7,
      [t, N(t)],
      0..5,
      color = red,
      thickness = 2,
      linestyle = dot
    )
  );
#
# "reflect" the "normal" plot and spoof the
# tickmarks
#
  display
  ( reflect
    ( display
      ( textplot([4.4,120,N(t)]),
        odeplot
        ( sol7,
          [t, N(t)],
          0..5,
          color = red,
          thickness = 2,
          linestyle = dot
        )
      ),
      [[2.5,0], [2.5,140]]
    ),
    tickmarks=[[seq( j=5-j, j=0..5)], default]
  );
#
# "reflect" the ODE, plot normally, and spoof
# the tickmarks
#
  sol8:=dsolve( [PDEtools:-dchange
                           ( t=5-tau,
                             D(N)(t) = N(t)
                           ),
                N(0)=rhs( sol7(5)[2])],
                numeric
              ):
  display
  ( textplot([0.6,120,N(t)]),
    odeplot
    ( sol8,
      [tau, N(tau)],
      0..5,
      color = red,
      thickness = 2,
      linestyle = dot
    ),
    tickmarks=[[seq( j=5-j, j=0..5)], default]
  );

>

Download oddPlots.mw

Both files...

Answered: tomleslie 13729

June 06 2023

0 0

suggested by Acer use the IterativeMaps() package, which was introduced in Maple 2015 (released 2015). If you really have Maple 18, that was released in 2014, and was the last release not to have the iterativeMaps capability

Can only suggest that you upload yoour worksheet here (using the big green up-arrow in the Mapleprimes toolbar) - then somebody *might* be able to solve your problem without the IterativeMaps() package

Shorter and neater...

Answered: tomleslie 13729

June 02 2023

0 0

I have edited this reply, because the attached is even shorter, neater and easier to understand

>

restart:

Test code written by Dr. Venkat Subramanian at UT Austin, 05/31/2023. This code uses CVP approach (piecwise constant) to perform optimal control. NLPSolve combination with dsolve numeric parametric form is buggy and fails for some values of nvar, and works for some values of nvar. Ideally increasing nvar should show convergence with respect to the objective function.

>

restart:

>	Digits:=15;

(1)

>	eqodes:=[diff(ca(t),t)=-(u+u^2/2)1.0ca(t),diff(cb(t),t)=1.0uca(t)-0.1*cb(t)];

(2)

>	soln:=dsolve({op(eqodes),ca(0)=alpha,cb(0)=beta},type=numeric,'parameters'=[alpha,beta,u],compile=true,savebinary=true):

>

  ss:= proc()
            interface(warnlevel=0):
            if   type( [_passed][1], numeric)
            then local z1,n1,i,c10,c20,dt,u, x:=[_passed];
                 global soln;
                 dt:=evalf(1.0/_npassed):
                 c10:=1.0:c20:=0.0:
                 for i from 1 to _npassed do
                     u:=x[i]:
                     soln('parameters'=[c10,c20,u]):
                     z1:=soln(dt):
                     c10:=subs(z1,ca(t)):c20:=subs(z1,cb(t)):
                 od:
                 -c20;
            else 'thisproc(_passed)'
            fi;
       end proc:

  getSol:= proc(N)
                local j;
                uses DirectSearch;
                Search( ss,
                        [seq( [z[j]>=0, z[j]<=5][], j=1..N)],
                        initialpoint=[0.1 $ N],
                        variables=[seq( z[j], j=1..N)]
                      )[1..2];
           end proc:
  getSol(2);
  getSol(3);
  getSol(4);
  getSol(5);

[-.5239011644866959, Vector(2, {(1) = .9099304633675794, (2) = 1.927217449401273})]

[-.5314533818996774, Vector(3, {(1) = .8114216459678198, (2) = 1.189408026035539, (3) = 2.427684833791241})]

[-.5351539426577706, Vector(4, {(1) = .7679420450967919, (2) = .9914999815821715, (3) = 1.4234427560343161, (4) = 2.850217883990661})]

[HFloat(-0.537338871817855), Vector[column](%id = 36893488148167978148)]

(3)

>

Download usingDS3.mw

And if you don't believe me...

Answered: tomleslie 13729

June 01 2023

0 0

consider the attached where the solution surface is plotted, together with the "result" which NLPSolve() finds (see previous post), but doesn't report - choosing instead to display an error?????

>

restart:

Test code written by Dr. Venkat Subramanian at UT Austin, 05/31/2023. This code uses CVP approach (piecwise constant) to perform optimal control. NLPSolve combination with dsolve numeric parametric form is buggy and fails for some values of nvar, and works for some values of nvar. Ideally increasing nvar should show convergence with respect to the objective function.

>

restart:

>	Digits:=15;

(1)

>	eqodes:=[diff(ca(t),t)=-(u+u^2/2)1.0ca(t),diff(cb(t),t)=1.0uca(t)-0.1*cb(t)];

(2)

>	soln:=dsolve({op(eqodes),ca(0)=alpha,cb(0)=beta},type=numeric,'parameters'=[alpha,beta,u],compile=true,savebinary=true):

>

ss:=proc(x)
interface(warnlevel=0):
#if type(x[1],numeric)
if type(x,Vector)
then

local z1,n1,i,c10,c20,dt,u;
global soln,nvar;
dt:=evalf(1.0/nvar):
c10:=1.0:c20:=0.0:
for i from 1 to nvar do
u:=x[i]:
soln('parameters'=[c10,c20,u]):
z1:=soln(dt):
c10:=subs(z1,ca(t)):c20:=subs(z1,cb(t)):
od:
#printf("%18.16f %18.16f %18.16f\n", x[1], x[2], -c20);
-c20;
else 'procname'(args):
end if:
end proc:

>	nvar:=2;#code works for nvar:=3, but not for nvar:=2

(3)

>	ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float[8]);

Vector(2, {(1) = .1, (2) = .1})

(4)

>	bl := Vector(nvar,[seq(0.,i=1..nvar)]);bu := Vector(nvar,[seq(5.,i=1..nvar)]);

Vector(2, {(1) = 0., (2) = 0.})

Vector[column](%id = 36893488148165423460)

(5)

>	infolevel[all]:=1;

(6)

>

plots:-display( [ plot3d( 'ss(<x,y>)', x=0..5, y=0..5),
                  plots:-pointplot3d( [0.9098284012521642, 1.9272152108626920, -0.5239011634950818],
                                      symbol=solidsphere,
                                      symbolsize=16,
                                      color=red
                                    )
                ]
              );

>

Download badNLPSolve2.mw

Sorry but I'm not buying this...

Answered: tomleslie 13729

June 01 2023

0 0

I accept that I screwed up on my early responses to this question - mea culpa etc, but I'm not buying this explanation of "nested tolerances" Something more fundamentqal is going on. It is relatively trivial to add a printf() statement to the main function evaluation procedure ss(). This allows one to see what values are being passed by NLPSolve() to the procedure ss(), and what values are being returned.

The attached shows the case for nvar=2. The fprintf() statement shows that for passed values of 0.9098284012521642 1.9272152108626920 then ss() will return -0.5239011634950818, which is a much "better" answer than the original supplied starting point

I have absolutely no idea why NLPSolve() reports this outcome as

Error, (in Optimization:-NLPSolve) no improved point could be found

when the initial point/return value was for passed values 0.1000000000000000, 0.1000000000000000 and return value -0.0902579301181078.

This anomaly has nothing to do with tolerances!! Looks more like a "bug" with NLPSolve

>

restart:

Test code written by Dr. Venkat Subramanian at UT Austin, 05/31/2023. This code uses CVP approach (piecwise constant) to perform optimal control. NLPSolve combination with dsolve numeric parametric form is buggy and fails for some values of nvar, and works for some values of nvar. Ideally increasing nvar should show convergence with respect to the objective function.

>

restart:

>	Digits:=15;

(1)

>	eqodes:=[diff(ca(t),t)=-(u+u^2/2)1.0ca(t),diff(cb(t),t)=1.0uca(t)-0.1*cb(t)];

(2)

>	soln:=dsolve({op(eqodes),ca(0)=alpha,cb(0)=beta},type=numeric,'parameters'=[alpha,beta,u],compile=true,savebinary=true):

>

ss:=proc(x)
interface(warnlevel=0):
#if type(x[1],numeric)
if type(x,Vector)
then

local z1,n1,i,c10,c20,dt,u;
global soln,nvar;
dt:=evalf(1.0/nvar):
c10:=1.0:c20:=0.0:
for i from 1 to nvar do
u:=x[i]:
soln('parameters'=[c10,c20,u]):
z1:=soln(dt):
c10:=subs(z1,ca(t)):c20:=subs(z1,cb(t)):
od:
printf("%18.16f %18.16f %18.16f\n", x[1], x[2], -c20);
-c20;
else 'procname'(args):
end if:
end proc:

>	nvar:=2;#code works for nvar:=3, but not for nvar:=2

(3)

>	ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float[8]);

Vector(2, {(1) = .1, (2) = .1})

(4)

>

ss(ic0);

0.1000000000000000 0.1000000000000000 -0.0902579301181078

(5)

>	bl := Vector(nvar,[seq(0.,i=1..nvar)]);bu := Vector(nvar,[seq(5.,i=1..nvar)]);

Vector(2, {(1) = 0., (2) = 0.})

Vector[column](%id = 36893488147938281100)

(6)

>	infolevel[all]:=1;

(7)

>	Optimization:-NLPSolve(nvar,ss,[],initialpoint=ic0,[bl,bu]);

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 2

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

0.1000000000000000 0.1000000000000000 -0.0902579301181078
0.1000000000000000 0.1000000000000000 -0.0902579301181078
0.1000000050000000 0.1000000000000000 -0.0902579321909171
0.0999999950000000 0.1000000000000000 -0.0902579280426463
0.1000000000000000 0.1000000000000000 -0.0902579301181078
0.1000000000000000 0.1000000050000000 -0.0902579323091734
0.1000000000000000 0.0999999950000000 -0.0902579279258143
0.5148270796300000 0.5383359092400000 -0.3644408224142451
0.5148270796300000 0.5383359092400000 -0.3644408224142451
0.5148270846300000 0.5383359092400000 -0.3644408234190942
0.5148270746300000 0.5383359092400000 -0.3644408214127478
0.5148270796300000 0.5383359092400000 -0.3644408224142451
0.5148270796300000 0.5383359142400000 -0.3644408236125698
0.5148270796300000 0.5383359042400000 -0.3644408212158754
0.9319439239912755 1.0363622309131530 -0.4885391290938756
0.9319439239912755 1.0363622309131530 -0.4885391290937617
0.9319439289912755 1.0363622309131530 -0.4885391292439379
0.9319439189912755 1.0363622309131530 -0.4885391289443778
0.9319439239912755 1.0363622309131530 -0.4885391290938755
0.9319439239912755 1.0363622359131530 -0.4885391295499261
0.9319439239912755 1.0363622259131530 -0.4885391286378250
1.0654644218475708 1.2730000574186295 -0.5064149446467445
1.0654644218475710 1.2730000574186295 -0.5064149446468069
1.0654644268475710 1.2730000574186295 -0.5064149445909414
1.0654644168475710 1.2730000574186295 -0.5064149447025125
1.0654644218475708 1.2730000574186300 -0.5064149446467442
1.0654644218475708 1.2730000624186300 -0.5064149449095958
1.0654644218475708 1.2730000524186300 -0.5064149443838926
1.1307175233927365 1.4869677975103834 -0.5136906299609928
1.1307175233927360 1.4869677975103834 -0.5136906299609961
1.1307175283927360 1.4869677975103834 -0.5136906298091304
1.1307175183927360 1.4869677975103834 -0.5136906301127701
1.1307175233927365 1.4869677975103830 -0.5136906299609926
1.1307175233927365 1.4869678025103830 -0.5136906301083968
1.1307175233927365 1.4869677925103830 -0.5136906298135882
1.1356641651835215 1.6702375806251024 -0.5175948676348219
1.1356641651835220 1.6702375806251024 -0.5175948676347549
1.1356641701835220 1.6702375806251024 -0.5175948674573128
1.1356641601835220 1.6702375806251024 -0.5175948678115854
1.1356641651835215 1.6702375806251020 -0.5175948676348220
1.1356641651835215 1.6702375856251020 -0.5175948677116450
1.1356641651835215 1.6702375756251020 -0.5175948675579989
1.0866483401716573 1.8748814282519424 -0.5210518726129031
1.0866483401716570 1.8748814282519424 -0.5210518726128130
1.0866483451716570 1.8748814282519424 -0.5210518724619215
1.0866483351716570 1.8748814282519424 -0.5210518727634627
1.0866483401716573 1.8748814282519420 -0.5210518726129032
1.0866483401716573 1.8748814332519420 -0.5210518726273581
1.0866483401716573 1.8748814232519420 -0.5210518725984485
1.0024722840877023 2.0052178656925972 -0.5229426723170228
1.0024722840877020 2.0052178656925972 -0.5229426723168223
1.0024722890877020 2.0052178656925972 -0.5229426722342294
1.0024722790877020 2.0052178656925972 -0.5229426723997439
1.0024722840877023 2.0052178656925970 -0.5229426723170227
1.0024722840877023 2.0052178706925970 -0.5229426722955273
1.0024722840877023 2.0052178606925970 -0.5229426723385182
0.9430320481867387 2.0086060956910257 -0.5236011432252972
0.9430320481867387 2.0086060956910257 -0.5236011432250253
0.9430320531867387 2.0086060956910257 -0.5236011431950941
0.9430320431867387 2.0086060956910257 -0.5236011432554187
0.9430320481867387 2.0086060956910260 -0.5236011432252974
0.9430320481867387 2.0086061006910260 -0.5236011432015695
0.9430320481867387 2.0086060906910260 -0.5236011432490251
0.9081145521282933 1.9530482453329765 -0.5238798233753532
0.9081145521282933 1.9530482453329765 -0.5238798233754096
0.9081145571282933 1.9530482453329765 -0.5238798233772498
0.9081145471282933 1.9530482453329765 -0.5238798233738003
0.9081145521282933 1.9530482453329760 -0.5238798233753530
0.9081145521282933 1.9530482503329760 -0.5238798233672611
0.9081145521282933 1.9530482403329760 -0.5238798233834446
0.9082511241039157 1.9278697744231208 -0.5239008814957404
0.9082511241039157 1.9278697744231208 -0.5239008814956164
0.9082511291039157 1.9278697744231208 -0.5239008814971667
0.9082511191039157 1.9278697744231208 -0.5239008814940147
0.9082511241039157 1.9278697744231210 -0.5239008814957403
0.9082511241039157 1.9278697794231210 -0.5239008814955319
0.9082511241039157 1.9278697694231210 -0.5239008814959485
0.9098284200385571 1.9272070230422622 -0.5239011634931299
0.9098284200385571 1.9272070230422622 -0.5239011634929417
0.9098284250385571 1.9272070230422622 -0.5239011634927924
0.9098284150385571 1.9272070230422622 -0.5239011634927948
0.9098284200385571 1.9272070230422620 -0.5239011634931297
0.9098284200385571 1.9272070280422620 -0.5239011634931323
0.9098284200385571 1.9272070180422620 -0.5239011634931271
0.9098272144282054 1.9272153377624019 -0.5239011634716886
0.9098272144282054 1.9272153377624019 -0.5239011634720514
0.9098272194282054 1.9272153377624019 -0.5239011634720205
0.9098272094282054 1.9272153377624019 -0.5239011634716544
0.9098272144282054 1.9272153377624020 -0.5239011634716885
0.9098272144282054 1.9272153427624020 -0.5239011634716885
0.9098272144282054 1.9272153327624020 -0.5239011634716886
0.9098283236129462 1.9272076880597515 -0.5239011634915934
0.9098283236129462 1.9272076880597515 -0.5239011634915982
0.9098283286129462 1.9272076880597515 -0.5239011634916982
0.9098283186129462 1.9272076880597515 -0.5239011634916264
0.9098283236129462 1.9272076880597520 -0.5239011634915935
0.9098283236129462 1.9272076930597520 -0.5239011634915961
0.9098283236129462 1.9272076830597520 -0.5239011634915909
0.9098284163373421 1.9272070485683923 -0.5239011634931292
0.9098284163373421 1.9272070485683923 -0.5239011634930038
0.9098284213373421 1.9272070485683923 -0.5239011634927385
0.9098284113373421 1.9272070485683923 -0.5239011634931530
0.9098284163373421 1.9272070485683920 -0.5239011634931291
0.9098284163373421 1.9272070535683920 -0.5239011634931318
0.9098284163373421 1.9272070435683920 -0.5239011634931264
0.9098284181879497 1.9272070358053273 -0.5239011634929194
0.9098284181879497 1.9272070358053273 -0.5239011634928628
0.9098284231879497 1.9272070358053273 -0.5239011634931987
0.9098284131879497 1.9272070358053273 -0.5239011634928153
0.9098284181879497 1.9272070358053270 -0.5239011634929196
0.9098284181879497 1.9272070408053270 -0.5239011634929222
0.9098284181879497 1.9272070308053270 -0.5239011634929166
0.9098284200381022 1.9272070230453999 -0.5239011634929444
0.9098284200381022 1.9272070230453999 -0.5239011634931870
0.9098284250381022 1.9272070230453999 -0.5239011634934626
0.9098284150381022 1.9272070230453999 -0.5239011634927948
0.9098284200381022 1.9272070230454000 -0.5239011634929444
0.9098284200381022 1.9272070280454000 -0.5239011634929475
0.9098284200381022 1.9272070180454000 -0.5239011634929419
0.9098284200384796 1.9272070230427967 -0.5239011634933074
0.9098284200384796 1.9272070230427967 -0.5239011634930013
0.9098284250384796 1.9272070230427967 -0.5239011634930425
0.9098284150384796 1.9272070230427967 -0.5239011634928509
0.9098284200384796 1.9272070230427970 -0.5239011634933077
0.9098284200384796 1.9272070280427970 -0.5239011634933104
0.9098284200384796 1.9272070180427970 -0.5239011634933047
0.9098281054318165 1.9273442248026260 -0.5239011629636382
0.9098281054318165 1.9273442248026260 -0.5239011629636305
0.9098281104318165 1.9273442248026260 -0.5239011629640344
0.9098281004318165 1.9273442248026260 -0.5239011629636130
0.9098281054318165 1.9273442248026260 -0.5239011629636383
0.9098281054318165 1.9273442298026260 -0.5239011629635976
0.9098281054318165 1.9273442198026260 -0.5239011629636792
0.9098284012400943 1.9272152211258320 -0.5239011634949601
0.9098284012400943 1.9272152211258320 -0.5239011634949032
0.9098284062400943 1.9272152211258320 -0.5239011634951821
0.9098283962400943 1.9272152211258320 -0.5239011634949272
0.9098284012400943 1.9272152211258320 -0.5239011634949600
0.9098284012400943 1.9272152261258320 -0.5239011634949601
0.9098284012400943 1.9272152161258320 -0.5239011634949597
0.9098284021037307 1.9272148445322757 -0.5239011634948513
0.9098284021037307 1.9272148445322757 -0.5239011634950216
0.9098284071037307 1.9272148445322757 -0.5239011634951234
0.9098283971037307 1.9272148445322757 -0.5239011634950456
0.9098284021037307 1.9272148445322760 -0.5239011634948513
0.9098284021037307 1.9272148495322760 -0.5239011634948515
0.9098284021037307 1.9272148395322760 -0.5239011634948513
0.9098284012521618 1.9272152158637716 -0.5239011634950859
0.9098284012521618 1.9272152158637716 -0.5239011634950884
0.9098284062521618 1.9272152158637716 -0.5239011634951252
0.9098283962521618 1.9272152158637716 -0.5239011634948711
0.9098284012521618 1.9272152158637720 -0.5239011634950861
0.9098284012521618 1.9272152208637720 -0.5239011634950859
0.9098284012521618 1.9272152108637720 -0.5239011634950856
0.9098284012521648 1.9272152158624476 -0.5239011634949724
0.9098284012521648 1.9272152158624476 -0.5239011634950814
0.9098284062521648 1.9272152158624476 -0.5239011634951218
0.9098283962521648 1.9272152158624476 -0.5239011634951661
0.9098284012521648 1.9272152158624480 -0.5239011634949721
0.9098284012521648 1.9272152208624480 -0.5239011634949726
0.9098284012521648 1.9272152108624480 -0.5239011634949721
0.9098284012521630 1.9272152158632316 -0.5239011634951454
0.9098284012521630 1.9272152158632316 -0.5239011634951454
0.9098284062521630 1.9272152158632316 -0.5239011634950614
0.9098283962521630 1.9272152158632316 -0.5239011634948064
0.9098284012521630 1.9272152158632320 -0.5239011634951457
0.9098284012521630 1.9272152208632320 -0.5239011634951450
0.9098284012521630 1.9272152108632320 -0.5239011634951454
0.9098284044033246 1.9272138416656097 -0.5239011634949304
0.9098284044033246 1.9272138416656097 -0.5239011634950476
0.9098284094033246 1.9272138416656097 -0.5239011634952592
0.9098283994033246 1.9272138416656097 -0.5239011634950145
0.9098284044033246 1.9272138416656100 -0.5239011634949305
0.9098284044033246 1.9272138466656100 -0.5239011634949308
0.9098284044033246 1.9272138366656100 -0.5239011634949298
0.9098284015080507 1.9272151042725862 -0.5239011634950941
0.9098284015080507 1.9272151042725862 -0.5239011634949092
0.9098284065080507 1.9272151042725862 -0.5239011634951867
0.9098283965080507 1.9272151042725862 -0.5239011634951776
0.9098284015080507 1.9272151042725860 -0.5239011634950940
0.9098284015080507 1.9272151092725860 -0.5239011634950941
0.9098284015080507 1.9272150992725860 -0.5239011634950937
0.9098284012600648 1.9272152124173210 -0.5239011634950286
0.9098284012600648 1.9272152124173210 -0.5239011634949610
0.9098284062600648 1.9272152124173210 -0.5239011634950114
0.9098283962600648 1.9272152124173210 -0.5239011634948073
0.9098284012600648 1.9272152124173210 -0.5239011634950287
0.9098284012600648 1.9272152174173210 -0.5239011634950291
0.9098284012600648 1.9272152074173210 -0.5239011634950286
0.9098284012521642 1.9272152158626916 -0.5239011634950819
0.9098284012521642 1.9272152158626916 -0.5239011634953270
0.9098284062521642 1.9272152158626916 -0.5239011634954255
0.9098283962521642 1.9272152158626916 -0.5239011634949881
0.9098284012521642 1.9272152158626920 -0.5239011634950818
0.9098284012521642 1.9272152208626920 -0.5239011634950820
0.9098284012521642 1.9272152108626920 -0.5239011634950818

Error, (in Optimization:-NLPSolve) no improved point could be found

>

Download badNLPSolve.mw

across...

Answered: tomleslie 13729

May 31 2023

0 0

all versions of Maple I can check, the first argument to the NLPSolve() command is the objective function whose extremum is required. In your worksheet, the first argument to the NLPSolve() command is either 2, 3 or 5 - none of which make much sense as an objectiev function!!!

I converted this to a question...

Answered: tomleslie 13729

May 31 2023

0 0

Before I start digging into this, what version of Maple are you running?

When I attempt to run the worksheet which you have posted, I get the error message

Error, (in Optimization:-NLPSolve) no improved point could be found

every time.

Furthermore I can find no call to NLPSolve(), which makes sense when the first argument is a "simple" number, rather than an expression/function whose minimum is being sought. In all three cases, you have

Optimization:-NLPSolve(nvar, ss........

where 'nvar' has been assigned to either 2, 3, or 5 - this makes no sense!!!

Provided...

Answered: tomleslie 13729

May 29 2023

0 0

that you can come up with a "good" starting point, then it is never a bad idea to give fsolve() as much help as possible. In you case fsolve() has no problems with calibrations 1,2,3 - but by giving it a "good" starting point, it will execute faster.

The reason...

Answered: tomleslie 13729

May 29 2023

0 0

that data "goes missing" is that the fsolve() command is failing to find solution for certain loop index values, at which point all subsequent code "breaks". In order to assist the fsolve() command as much as possible, there are a couple of possibilities

one strategy is to restrict the region over which it has to search for a solution
second strategy is to provide fsolve with a "starting point" which is "close" to the desired answer

The attached uses (2) above, where the starting point for the i-th loop iteration is the result of the (i-1)-th loop iteration. The initial gues for the first iteration in each calibration, I took as {lambda__1 = 0.25, lambda__2 = 0.25, lambda__3 = 0.25}. This strategy generates all values for calibrations 7 and 8. (See the attached)

Interestingly in calibration 7, lambda__1 and lambda__3 vary quite a lot over the range of the loop (ie ~0..~100), whereas lambda__2 doesn't vary much at all. In calibration 8, lambda__2 and lambda__3 vary quite a lot over the range of the loop (ie ~0..~100), whereas lambda__1 doesn't vary much at all.

>

restart;

>	_local(gamma);

(1)

#Define the assumptions ex-ante (variances as real and positive, correlations in between -1 and +1 and so on...) - or Maple wouldn't know

>

assume(gamma::real, lambda__1::real, lambda__2::real, lambda__3::real, sigma__epsilon1::real, sigma__epsilon2::real, nu__0[1]::real, nu__0[2]::real, rho__u[1, 2]::real, rho__u[1, 3]::real, rho__u[2, 3]::real, rho__v[1, 2]::real, sigma__u[1]::real, sigma__u[2]::real, sigma__u[3]::real, sigma__v[1]::real, sigma__v[2]::real);

>

assume(0 <= gamma, 0 <= lambda__1, 0 <= lambda__2, 0 <= lambda__3, 0 <= sigma__epsilon1, 0 <= sigma__epsilon2, 0 <= nu__0[1], 0 <= nu__0[2], -1 <= rho__u[1, 2] and rho__u[1, 2] <= 1, -1 <= rho__u[1, 3] and rho__u[1, 3] <= 1, -1 <= rho__u[2, 3] and rho__u[2, 3] <= 1, -1 <= rho__v[1, 2] and rho__v[1, 2] <= 1, 0 <= sigma__u[1], 0 <= sigma__u[2], 0 <= sigma__u[3], 0 <= sigma__v[1], 0 <= sigma__v[2]);

1. FIRST PART

>	E := bExp_nu1_S + cExp_nu2_S + a;

(2)

>	V := 2Cov_u12de + 2Cov_u13df + 2Cov_u23ef + Var_u1d^2 + Var_u2e^2 + Var_u3f^2 + b^2Var_nu1_S + 2bcCov_nu12_S + c^2*Var_nu2_S;

(3)

>	argmin := E - gamma/2*V;

b*Exp_nu1_S+c*Exp_nu2_S+a-(1/2)*gamma*(2*Cov_u12*d*e+2*Cov_u13*d*f+2*Cov_u23*e*f+Var_u1*d^2+Var_u2*e^2+Var_u3*f^2+b^2*Var_nu1_S+2*b*c*Cov_nu12_S+c^2*Var_nu2_S)

(4)

>

a := X__1*(-X__1*lambda__1 - nu__0[1]) + X__2*(-X__2*lambda__2 - nu__0[2]) + X__3*(-X__3*lambda__3 + (-nu__0[1] - nu__0[2]));
b := X__1 + X__3;
c := X__2 + X__3;
d := -X__1*lambda__1;
e := -X__2*lambda__2;
f := -X__3*lambda__3;

Exp_nu1_S := SIG1_m_nu01*theta__11 + SIG2_m_nu02*theta__12 + nu__0[1];
Exp_nu2_S := SIG1_m_nu01*theta__21 + SIG2_m_nu02*theta__22 + nu__0[2];

(5)

>	FOC_1 := diff(argmin, X__1) = 0;

SIG1_m_nu01*theta__11+SIG2_m_nu02*theta__12-2*X__1*lambda__1-(1/2)*gamma*(2*Cov_u12*lambda__1*X__2*lambda__2+2*Cov_u13*lambda__1*X__3*lambda__3+2*Var_u1*X__1*lambda__1^2+2*(X__1+X__3)*Var_nu1_S+2*(X__2+X__3)*Cov_nu12_S) = 0

(6)

>	FOC_2 := diff(argmin, X__2) = 0;

SIG1_m_nu01*theta__21+SIG2_m_nu02*theta__22-2*X__2*lambda__2-(1/2)*gamma*(2*Cov_u12*X__1*lambda__1*lambda__2+2*Cov_u23*lambda__2*X__3*lambda__3+2*Var_u2*X__2*lambda__2^2+2*(X__1+X__3)*Cov_nu12_S+2*(X__2+X__3)*Var_nu2_S) = 0

(7)

>	FOC_3 := diff(argmin, X__3) = 0;

SIG1_m_nu01*theta__11+SIG2_m_nu02*theta__12+SIG1_m_nu01*theta__21+SIG2_m_nu02*theta__22-2*X__3*lambda__3-(1/2)*gamma*(2*Cov_u13*X__1*lambda__1*lambda__3+2*Cov_u23*X__2*lambda__2*lambda__3+2*Var_u3*X__3*lambda__3^2+2*(X__1+X__3)*Var_nu1_S+2*(X__2+X__3)*Cov_nu12_S+2*(X__1+X__3)*Cov_nu12_S+2*(X__2+X__3)*Var_nu2_S) = 0

(8)

>	$strat := convert(solve({FOC_1, FOC_2, FOC_3}, {X__1, X__2, X__3}), list)$

>	# First strategy

${Cov_u12, Cov_u13, Cov_u23, Var_u1, Var_u2, Var_u3, gamma, Cov_nu12_S, SIG1_m_nu01, SIG2_m_nu02, Var_nu1_S, Var_nu2_S, lambda__1, lambda__2, lambda__3, theta__11, theta__12, theta__21, theta__22}$

(9)

>	# Second strategy

${Cov_u12, Cov_u13, Cov_u23, Var_u1, Var_u2, Var_u3, gamma, Cov_nu12_S, SIG1_m_nu01, SIG2_m_nu02, Var_nu1_S, Var_nu2_S, lambda__1, lambda__2, lambda__3, theta__11, theta__12, theta__21, theta__22}$

(10)

>	# Third strategy

${Cov_u12, Cov_u13, Cov_u23, Var_u1, Var_u2, Var_u3, gamma, Cov_nu12_S, SIG1_m_nu01, SIG2_m_nu02, Var_nu1_S, Var_nu2_S, lambda__1, lambda__2, lambda__3, theta__11, theta__12, theta__21, theta__22}$

(11)

2. SECOND PART

>

# Check linear projection theorem's calculations first:
eql1 := (beta__11*Var_nu1+beta__12*Cov_nu12)/(((beta__11)^2)*Var_S1+((beta__12)^2)*Var_S2+2*beta__11*beta__12*Cov_S12+Var_u1);

# Now plug in the betas:
beta__11 := simplify(coeff(Strategy1, SIG1_m_nu01)):
beta__12 := simplify(coeff(Strategy1, SIG2_m_nu02)):

# Check linear projection theorem's calculations first:
eql2 := (beta__21*Cov_nu12+beta__22*Var_nu2)/(((beta__21)^2)*Var_S1+((beta__22)^2)*Var_S2+2*beta__21*beta__22*Cov_S12+Var_u2);

# Now plug in the betas:
beta__21 := simplify(coeff(Strategy2, SIG1_m_nu01)):
beta__22 := simplify(coeff(Strategy2, SIG2_m_nu02)):

# Check linear projection theorem's calculations first:
eql3 := (beta__31*Var_nu1+(beta__31+beta__32)*Cov_nu12+beta__32*Var_nu2)/(((beta__31)^2)*Var_S1+((beta__32)^2)*Var_S2+2*beta__31*beta__32*Cov_S12+Var_u3);

# Now plug in the betas:
beta__31 := simplify(coeff(Strategy3, SIG1_m_nu01)):
beta__32 := simplify(coeff(Strategy3, SIG2_m_nu02)):

(Var_nu1*beta__11+Cov_nu12*beta__12)/(2*Cov_S12*beta__11*beta__12+Var_S1*beta__11^2+Var_S2*beta__12^2+Var_u1)

(Var_nu2*beta__22+Cov_nu12*beta__21)/(2*Cov_S12*beta__21*beta__22+Var_S1*beta__21^2+Var_S2*beta__22^2+Var_u2)

(beta__31*Var_nu1+(beta__31+beta__32)*Cov_nu12+beta__32*Var_nu2)/(2*Cov_S12*beta__31*beta__32+Var_S1*beta__31^2+Var_S2*beta__32^2+Var_u3)

(12)

>	MyEqs := {eql1 - lambda__1, eql2 - lambda__2, eql3 - lambda__3}: MyVars := {lambda__1, lambda__2, lambda__3}: # [i = {....}] means: the equation number i contains this set of variables print~([ seq([i=indets(MyEqs[i], name) intersect MyVars], i=1..nops(MyEqs)) ]):

$[1 = {lambda__1, lambda__2, lambda__3}]$

$[2 = {lambda__1, lambda__2, lambda__3}]$

$[3 = {lambda__1, lambda__2, lambda__3}]$

(13)

>

# General Calibration:

P := indets(MyEqs, name) minus MyVars;

d := 1/(((sigma__v[1])^2+(sigma__epsilon1)^2)*((sigma__v[2])^2+(sigma__epsilon2)^2)-(sigma__v[1])^2*(rho__v[1, 2])^2*(sigma__v[2])^2);

cal := [
  Cov_S12 = sigma__v[1]*rho__v[1, 2]*sigma__v[2],
  Cov_u12 = sigma__u[1]*rho__u[1, 2]*sigma__u[2],
  Cov_u13 = sigma__u[1]*rho__u[1, 3]*sigma__u[3],
  Cov_u23 = sigma__u[2]*rho__u[2, 3]*sigma__u[3],
  Var_S1 = (sigma__v[1])^2+(sigma__epsilon1)^2,
  Var_S2 = (sigma__v[2])^2+(sigma__epsilon2)^2,
  Var_nu1 = (sigma__v[1])^2,
  Var_nu2 = (sigma__v[2])^2,
  Var_u1 = (sigma__u[1])^2,
  Var_u2 = (sigma__u[2])^2,
  Var_u3 = (sigma__u[3])^2,
  Cov_nu12 = sigma__v[1]*rho__v[1, 2]*sigma__v[2],
  Cov_nu12_S = d*sigma__v[1]*rho__v[1, 2]*sigma__v[2]*(sigma__epsilon1)^2*(sigma__epsilon2)^2,
  Var_nu1_S = d*(sigma__epsilon1)^2*((sigma__v[1])^2*((sigma__v[2])^2+(sigma__epsilon2)^2)-(sigma__v[1])^2*(rho__v[1, 2])^2*(sigma__v[2])^2),
  Var_nu2_S = d*(sigma__epsilon2)^2*((sigma__v[2])^2*((sigma__v[1])^2+(sigma__epsilon1)^2)-(sigma__v[1])^2*(rho__v[1, 2])^2*(sigma__v[2])^2),
  theta__11 = 1-d*((sigma__epsilon1)^2*((sigma__v[2])^2+(sigma__epsilon2)^2)),
  theta__12 = d*sigma__v[1]*rho__v[1, 2]*sigma__v[2]*(sigma__epsilon1)^2,
  theta__21 = d*sigma__v[1]*rho__v[1, 2]*sigma__v[2]*(sigma__epsilon2)^2,
  theta__22 = 1-d*((sigma__epsilon2)^2*((sigma__v[1])^2+(sigma__epsilon1)^2))
];

MyEqs_cal := eval(MyEqs,cal):
length(MyEqs_cal);

P_cal := indets(MyEqs_cal, name) minus MyVars;

${Cov_S12, Cov_u12, Cov_u13, Cov_u23, Var_S1, Var_S2, Var_nu1, Var_nu2, Var_u1, Var_u2, Var_u3, gamma, Cov_nu12, Cov_nu12_S, Var_nu1_S, Var_nu2_S, theta__11, theta__12, theta__21, theta__22}$

1/((sigma__epsilon1^2+sigma__v[1]^2)*(sigma__epsilon2^2+sigma__v[2]^2)-sigma__v[1]^2*rho__v[1, 2]^2*sigma__v[2]^2)

${gamma, sigma__epsilon1, sigma__epsilon2, rho__u[1, 2], rho__u[1, 3], rho__u[2, 3], rho__v[1, 2], sigma__u[1], sigma__u[2], sigma__u[3], sigma__v[1], sigma__v[2]}$

(14)

>

# Numerical Calibrations:

P_cal_cor := P_cal minus {gamma,sigma__epsilon1,sigma__epsilon2,sigma__v[1],sigma__v[2],sigma__u[1],sigma__u[2],sigma__u[3]};

P_cal_var := P_cal minus {rho__v[1, 2],rho__u[1, 2],rho__u[1, 3],rho__u[2, 3]};

zerocor := {rho__v[1, 2],rho__u[1, 2],rho__u[1, 3],rho__u[2, 3]}=~0:

normvar := {gamma,sigma__epsilon1,sigma__epsilon2,sigma__v[1],sigma__v[2],sigma__u[1],sigma__u[2],sigma__u[3]}=~1:

ncal1 := ({gamma = INDEX} union (P_cal_var minus {gamma} =~ 1)) union zerocor;
ncal2 := ({sigma__epsilon1 = INDEX} union (P_cal_var minus {sigma__epsilon1} =~ 1)) union zerocor;
ncal3 := ({sigma__epsilon2 = INDEX} union (P_cal_var minus {sigma__epsilon2} =~ 1)) union zerocor;
ncal4 := ({sigma__u[1] = INDEX} union (P_cal_var minus {sigma__u[1]} =~ 1)) union zerocor;
ncal5 := ({sigma__u[2] = INDEX} union (P_cal_var minus {sigma__u[2]} =~ 1)) union zerocor;
ncal6 := ({sigma__u[3] = INDEX} union (P_cal_var minus {sigma__u[3]} =~ 1)) union zerocor;
ncal7 := ({sigma__v[1] = INDEX} union (P_cal_var minus {sigma__v[1]} =~ 1)) union zerocor;
ncal8 := ({sigma__v[2] = INDEX} union (P_cal_var minus {sigma__v[2]} =~ 1)) union zerocor;
ncal9 := ({rho__u[1, 2] = INDEX} union (P_cal_cor minus {rho__u[1, 2]} =~ 0)) union normvar;
ncal10 := ({rho__u[1, 3] = INDEX} union (P_cal_cor minus {rho__u[1, 3]} =~ 0)) union normvar;
ncal11 := ({rho__u[2, 3] = INDEX} union (P_cal_cor minus {rho__u[2, 3]} =~ 0)) union normvar;
ncal12 := ({rho__v[1, 2] = INDEX} union (P_cal_cor minus {rho__v[1, 2]} =~ 0)) union normvar;

${rho__u[1, 2], rho__u[1, 3], rho__u[2, 3], rho__v[1, 2]}$

${gamma, sigma__epsilon1, sigma__epsilon2, sigma__u[1], sigma__u[2], sigma__u[3], sigma__v[1], sigma__v[2]}$

${gamma = INDEX, sigma__epsilon1 = 1, sigma__epsilon2 = 1, rho__u[1, 2] = 0, rho__u[1, 3] = 0, rho__u[2, 3] = 0, rho__v[1, 2] = 0, sigma__u[1] = 1, sigma__u[2] = 1, sigma__u[3] = 1, sigma__v[1] = 1, sigma__v[2] = 1}$

${gamma = 1, sigma__epsilon1 = INDEX, sigma__epsilon2 = 1, rho__u[1, 2] = 0, rho__u[1, 3] = 0, rho__u[2, 3] = 0, rho__v[1, 2] = 0, sigma__u[1] = 1, sigma__u[2] = 1, sigma__u[3] = 1, sigma__v[1] = 1, sigma__v[2] = 1}$

${gamma = 1, sigma__epsilon1 = 1, sigma__epsilon2 = INDEX, rho__u[1, 2] = 0, rho__u[1, 3] = 0, rho__u[2, 3] = 0, rho__v[1, 2] = 0, sigma__u[1] = 1, sigma__u[2] = 1, sigma__u[3] = 1, sigma__v[1] = 1, sigma__v[2] = 1}$

${gamma = 1, sigma__epsilon1 = 1, sigma__epsilon2 = 1, rho__u[1, 2] = 0, rho__u[1, 3] = 0, rho__u[2, 3] = 0, rho__v[1, 2] = 0, sigma__u[1] = INDEX, sigma__u[2] = 1, sigma__u[3] = 1, sigma__v[1] = 1, sigma__v[2] = 1}$

${gamma = 1, sigma__epsilon1 = 1, sigma__epsilon2 = 1, rho__u[1, 2] = 0, rho__u[1, 3] = 0, rho__u[2, 3] = 0, rho__v[1, 2] = 0, sigma__u[1] = 1, sigma__u[2] = INDEX, sigma__u[3] = 1, sigma__v[1] = 1, sigma__v[2] = 1}$

${gamma = 1, sigma__epsilon1 = 1, sigma__epsilon2 = 1, rho__u[1, 2] = 0, rho__u[1, 3] = 0, rho__u[2, 3] = 0, rho__v[1, 2] = 0, sigma__u[1] = 1, sigma__u[2] = 1, sigma__u[3] = INDEX, sigma__v[1] = 1, sigma__v[2] = 1}$

${gamma = 1, sigma__epsilon1 = 1, sigma__epsilon2 = 1, rho__u[1, 2] = 0, rho__u[1, 3] = 0, rho__u[2, 3] = 0, rho__v[1, 2] = 0, sigma__u[1] = 1, sigma__u[2] = 1, sigma__u[3] = 1, sigma__v[1] = INDEX, sigma__v[2] = 1}$

${gamma = 1, sigma__epsilon1 = 1, sigma__epsilon2 = 1, rho__u[1, 2] = 0, rho__u[1, 3] = 0, rho__u[2, 3] = 0, rho__v[1, 2] = 0, sigma__u[1] = 1, sigma__u[2] = 1, sigma__u[3] = 1, sigma__v[1] = 1, sigma__v[2] = INDEX}$

${gamma = 1, sigma__epsilon1 = 1, sigma__epsilon2 = 1, rho__u[1, 2] = INDEX, rho__u[1, 3] = 0, rho__u[2, 3] = 0, rho__v[1, 2] = 0, sigma__u[1] = 1, sigma__u[2] = 1, sigma__u[3] = 1, sigma__v[1] = 1, sigma__v[2] = 1}$

${gamma = 1, sigma__epsilon1 = 1, sigma__epsilon2 = 1, rho__u[1, 2] = 0, rho__u[1, 3] = INDEX, rho__u[2, 3] = 0, rho__v[1, 2] = 0, sigma__u[1] = 1, sigma__u[2] = 1, sigma__u[3] = 1, sigma__v[1] = 1, sigma__v[2] = 1}$

${gamma = 1, sigma__epsilon1 = 1, sigma__epsilon2 = 1, rho__u[1, 2] = 0, rho__u[1, 3] = 0, rho__u[2, 3] = INDEX, rho__v[1, 2] = 0, sigma__u[1] = 1, sigma__u[2] = 1, sigma__u[3] = 1, sigma__v[1] = 1, sigma__v[2] = 1}$

${gamma = 1, sigma__epsilon1 = 1, sigma__epsilon2 = 1, rho__u[1, 2] = 0, rho__u[1, 3] = 0, rho__u[2, 3] = 0, rho__v[1, 2] = INDEX, sigma__u[1] = 1, sigma__u[2] = 1, sigma__u[3] = 1, sigma__v[1] = 1, sigma__v[2] = 1}$

(15)

>	betseq:=["beta__11", "beta__12", "beta__21", "beta__22", "beta__31", "beta__32"]: lamseq:=["lambda1", "lambda2", "lambda3"]:

>

  lam:={lambda__1 = 0.25, lambda__2 = 0.25, lambda__3 = 0.25}:
  lamRes1:=Array(1..1000, 1..3):
  betRes1:=Array(1..1000, 1..6):
  for INDEX from 1 by 1 to 1000 do
      lam := fsolve(eval(MyEqs_cal, ncal1), lam);
      betRes1[INDEX]:= Array([seq(eval(eval(j, cal), lam union ncal1), j in parse~(betseq))]);
      lamRes1[INDEX]:= Array(rhs~(convert(lam,list))):
  od:

>

  lam:={lambda__1 = 0.25, lambda__2 = 0.25, lambda__3 = 0.25}:
  lamRes2:=Array(1..1000, 1..3):
  betRes2:=Array(1..1000, 1..6):
  for INDEX from 1 by 1 to 1000 do
      lam := fsolve(eval(MyEqs_cal, ncal2),lam);
      betRes2[INDEX]:= Array([seq(eval(eval(j, cal), lam union ncal2), j in parse~(betseq))]);
      lamRes2[INDEX]:= Array(rhs~(convert(lam,list))):
  od:

>

  lam:={lambda__1 = 0.25, lambda__2 = 0.25, lambda__3 = 0.25}:
  lamRes3:=Array(1..1000, 1..3):
  betRes3:=Array(1..1000, 1..6):
  for INDEX from 1 by 1 to 1000 do
      lam := fsolve(eval(MyEqs_cal, ncal3),lam);
      betRes3[INDEX]:= Array([seq(eval(eval(j, cal), lam union ncal3), j in parse~(betseq))]);
      lamRes3[INDEX]:= Array(rhs~(convert(lam,list))):
  od:

>

  lam:={lambda__1 = 0.25, lambda__2 = 0.25, lambda__3 = 0.25}:
  lamRes7:=Array(1..1000, 1..3):
  betRes7:=Array(1..1000, 1..6):
  for INDEX from 1 by 1 to 1000 do
      lam := fsolve(eval(MyEqs_cal, ncal7), lam);
      betRes7[INDEX]:= Array([seq(eval(eval(j, cal), lam union ncal7), j in parse~(betseq))]);
      lamRes7[INDEX]:= Array(rhs~(convert(lam,list))):
  od:

>

  lam:={lambda__1 = 0.25, lambda__2 = 0.25, lambda__3 = 0.25}:
  lamRes8:=Array(1..1000, 1..3):
  betRes8:=Array(1..1000, 1..6):
  for INDEX from 1 by 1 to 1000 do
      lam := fsolve(eval(MyEqs_cal, ncal8), lam);
      betRes8[INDEX]:= Array([seq(eval(eval(j, cal), lam union ncal8), j in parse~(betseq))]);
      lamRes8[INDEX]:= Array(rhs~(convert(lam,list))):
  od:

>

  cols:=[black, blue, cyan, green, yellow, red]:
  plots:-display( seq( plot( [$ 20], lamRes1[1..20,j], color=cols[j]), j=1..3), legend=lamseq);
  plots:-display( seq( plot( [$ 20], lamRes2[1..20,j], color=cols[j]), j=1..3), legend=lamseq);
  plots:-display( seq( plot( [$ 20], lamRes3[1..20,j], color=cols[j]), j=1..3), legend=lamseq);
  plots:-display( seq( plot( [$ 20], betRes1[1..20,j], color=cols[j]), j=1..6), legend=betseq);
  plots:-display( seq( plot( [$ 20], betRes2[1..20,j], color=cols[j]), j=1..6), legend=betseq);
  plots:-display( seq( plot( [$ 20], betRes3[1..20,j], color=cols[j]), j=1..6), legend=betseq);

>

  lamRes7[1000, 1..3];
  lamRes8[1000, 1..3];
  plots:-display( seq( plot( [$ 1000], lamRes7[1..-1,j], color=cols[j]), j=1..3), legend=lamseq);
  plots:-display( seq( plot( [$ 1000], lamRes8[1..-1,j], color=cols[j]), j=1..3), legend=lamseq);
  plots:-display( seq( plot( [$ 100], betRes7[1..100,j], color=cols[j]), j=1..6), legend=betseq);
  plots:-display( seq( plot( [$ 100], betRes8[1..100,j], color=cols[j]), j=1..6), legend=betseq);