## 6543 Reputation

10 years, 142 days

## Didn't I fix this a while ago?...

I thiught my response in the following thread

https://www.mapleprimes.com/questions/229996-The-Problem-With-Updating-The-Physics-Package

Did it, or didn't it?

I did warn you in that thread, every time you try to update the Physics Package, you are probably going to have to repeat the  process I outlined. And all because the MapleSoft update process is completely incapable of making a couple of files "appear" in the correct places

Still, it does give one a lot of confidence about MapleSoft's ability to produce reliable software - NOT

## Can't check anything in Maple 13...

because, I don't have anything so OLD.

Just a suggestion - the problem may be with the and() function, so I might be tempted to try something involving a simple algebraic expression, rather than a boolean - such as

event1 := [ gln^2+gln^2, halt]]

## Well...

I seem to have fixed most of the problems in your original question, and the techniques I have used are generally applicable to any problem.

But now you come up with a new issue - that it is "not a general procedure yet and making one could be complicated?"

Actually I disagree - it is trivial to convert my original worksheet into a more "general" procedure, but a lot depends on exactly what kind of input data you want to handle.

Consider the attached

 > restart;   with(plots):   doPlot:= proc( A )                  display                  ( [ seq                      ( seq                        ( plot3d                          ( A[i,j],                            x=i-1..i,                            y=j-1..j,                            shading=zhue,                            style=surface,                            axes=normal,                            view=[ 0..op([2, 1, 2], A),                                   0..op([2, 2, 2], A),                                   min(0, A[..,2] )..max(A[..,2])                                 ]                          ),                          i=1..op([2,1, 2], A)                        ),                        j=1..2                      )                    ],                    scaling=constrained,                    size=[1000, 1000]                  );           end proc:   T1:=Array( [ [ 0.5, -1  ],                [ 1,    2  ]              ]            ):   T2:=Array( [ [ 0.5, -1  ],                [ 1,    2  ],                [ 1.5,  1  ],                [ 1.75, 2  ],                [ 2,    2.5]             ]           ):   doPlot(T1);   doPlot(T2);  >

## Something like...

the attached perhaps?

 > restart;   with(plots): with(plottools): with(ColorTools): with(combinat):
 > # # Number of terms in the fibonacci sequence. Ajust # as required #   N:=24: # # A couple of utilities #   g:= k-> round~( convert( op(1, a[k])[-1,..], list)):   h:= k-> round~( convert( op(1, a[k])[1,..], list)): # # Initialise the first arc and rectangle #   a:= arc( [0,0], fibonacci(1),0..Pi/2):   r:= rectangle( g(1), h(1)):
 > # # Loop through all subsequent arc sections and # rectangles #   for j from 2 by 1 to N do       fj:=fibonacci(j):     #     # Generate the arc section     #       a[j]:= arc              ( [ g(j-1)+cos((j+1)*Pi/2)*fj,                  g(j-1)+sin((j+1)*Pi/2)*fj                ],                             fj,                (j-1)*Pi/2..j*Pi/2              );     #     # and the corresponding rectangle     #       r[j]:= rectangle              ( g(j),                h(j),                color=Color( [ seq( rand()/10^12, k=1..3 ) ] )              );   od:
 > # # Animate both the arc and the rectangles of the spiral #   display   ( [ seq       ( display         ( [ seq( [a[j],r[j]][], j=1..i)],            caption=typeset( "Fibonacci number is %1",                             fibonacci(i)                           ),            captionfont=[times, bold, 24]         ),         i=1..N       )     ],     titlefont=[times, bold, 24],     scaling=constrained,     axes=none,     insequence=true,     size=[1200,1200]   ); >

## Use ImportMatrix...

which at least seems to make both data sets the same. I'd upload an example, but when I try this site is giving me a "secure server error"

## Be aware...

Next time you update your Physics - you may have to play a similar game - unless Maplesoft actually get their act together!

## Need answers to basic questions...

1. In your worksheet, you imply that the lists aca 'times' and 'C__f'' contain experimetal data. So why are these lists of different lengths? In the optimisation process you are only using the first ~20 entries in the C__f list, ie the values 1, 1, [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 8]. Is this intentional!!!!!
2. Let us consider the first entries in these two lists, from the 'times' list we get 5 and from the C__f list we get 1. This suggest that for small values of the independent variable 'T' yoiu expect small (ie~1) values of the independent varaible C(T). However your ODE system contains the initial condition C(0)=160000.
3. In other words you are taking the values of the independent variable given by the list [5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100], with corresponding values of the dependent variable [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 8] and you want to 'fit' these to a function C(T) which has the condition C(0)=160000. Does this sound even remotely sensible to you?
4. In the attached I have used the (free) DirectSearch() package to investigate this issue further. You (probably) will not have the DirectSearch package installed - so you will not be able to execute this worksheet. This worksheet demonstrates two things
1. with the initial values of the parameters l, m, n, q, r, u, v, sigma, iota, nu, phi, upsilon, delta, g which you supply, the fit  is absoutely terrible, which can be seen in the attached both from the residual error and by comparing graphs of the "experimental data" and the ODE solution for C(T)
2. with the optimised values of all the parameters returned by DirectSearch(), the value of the residual error is greatly reduced - but the overall fit is still terrible, as can be seen from the graphs.

I think you seriously have to address the issue of the mismatch between the initial condition for the ODE system, ie C(0)=160000, and the values of your 'experimental data' for low values of T

Bearing in mnd all of the above - I would say that what the attached produces is pretty much rubbish - although it does 'execute' provided you have the DirectSearch package installed

 > > Parameterizing (with respect to ) the numerical solution of the model deq       #Calculating the sum of the square of the errors between the model predictions and experimental data      #Minimizing the sum of the square of the errors to find the best fit values of   (1)
 >  (2)
 > > >  (3)
 > > # # Use the OP's list of initial parameter values #   initVars:=[ l = 0.05, m = 0.02, n = 0.017, q = 0.004,               r = 0.098710, u = 0.0100, v = 0.08543,               sigma = 0.0349, iota = 0.0032, nu = 0.0014,               phi = 0.931, upsilon = 0.0019, delta = 0.01,               g = 0.3]:   res( parameters= rhs~(initVars) ); # # Plot the value of the ODE solution with these # initial parameters #   plots:-odeplot( res, [T,C(T)], T=5..100, axes=boxed); # # plot the pairs [times[j], C__f[j]] - ie the OP's # "experimental" values - for comparison with the # ODE solution above #   plot( [seq( [times[j],C__f[j]], j=1..numelems(times))]);   sse( rhs~(initVars)[]);    (4)
 > # # Use DirectSearch() to obtain an "optimal" set of parameters # NB OP *probably* won't be able to do this because (s)he wont # have the DirectSearch() package installed. #   DSOpts:=DirectSearch:-Search( sse,initialpoint = initVars,assume=positive); (5)
 > # # Double-check the error value returned when using the # "optimised" parameter values found by DirectSearch() # (Should return the first value in the DSOpts() Array # - which it does   sse( entries(DSOpts, nolist)); # # Plot the value of the ODE solution with these # initial parameters #   res(parameters= [entries(DSOpts, nolist)]):   plots:-odeplot( res, [T,C(T)], T=5..100, axes=boxed);  > ## Don't...

post "pictures" of code.

They cannot be executed, cannot be debugged and in fact are almost completely useless - and before you ask, no-one on this site is going to read your code and retype it.

After all why should they when you can upload executable, debuggable code by simply using the big green up-arrow in the Mapleprimes toolbar

Obviously I don't know the complexity of the problem you are trying to address, but for reasonably simple circuits, I would have thought that it would be adequate.

Assuming you have Syrup installed, then you can access thethe general Syrup help by typing

?Syrup

at the Maple command prompt.

Since you ask for examples, you can pull up a pagefull of these of these just by typing

?Syrup/Examples

at the Maple command prompt

## Then I can only suggest...

that you re-read the above-referenced table again

if you could read the response you have already been given - it means that I doin't have to keep restating the blindingly obvious. So when you state

Looking again to the code.

I aspected that yv  -values of f  should be stored in a array, but that happened here.

yv:=unapply( f, indets(f, 'name')[])~(xv);

But there is not Array made for this at forehand , well perhaps it has to do with the  ~ (xv) operator
It depend then from what datacontainer type( source container)  is used by ~(xv) ?

In this case it was a Array ,but was it a list then you get probabably a target datacontainer  list?   Precisley whihc part of this question is NOT covered by the comment from my previous post which states

 # # In a similar way, this (nameless) function can be # applied 'elementwise' to every element in a "container" # - ie a list, an array, vector etc, by using the '~' # operator. This will return a container of the same type, # but with updated entries.

And if you really want to understand this in detail then you read the help page at

help(elementwise);

## Depend on where yu want to export to...

and in what format - and only you know that.

If Excel is the desired target, then the simplest way is probably to use the ExcelTools:-Export() command. So adding the command

ExcelTools:-Export(M, "C:/Users/TomLeslie/Desktop/test.xlsx") ;

to the bottom of the worksheet I posted previously will export the matrix 'M' to the Excel file "test.xlsx" at the location "C:/Users/TomLeslie/Desktop". The last of these happens to be my default "desktop" folder. You will have to change this location to something appropriate for your installation

## Step by step...

through the command

unapply( f, indets(f, 'name')[])~(..)

with appropriate help pages

 > # # Now suppose you have an expression #   f:=x^2+2*x+1; (1)
 > # # indets(f, 'name') returns all of the unevaluated # names in the expression 'f' as a set. In this example # the set only containes one entry. #   help(indets);   indets(f, 'name'); (2)
 > # # On this occasion, we don't want this command # to return a 'set', we just want the 'entries'. # There are (at least) a couple of ways to do this # # So all this part does is to return the 'unknown' # variable in the original expression 'f'. #   help(list);   op(indets(f, 'name')); #or   indets(f, 'name')[];  (3)
 > # # By using uanpply(f , indets(f, 'name')[]) the # original supplied expression is converted to an # appliable (nameless) function, as in #   help(unapply);   unapply( f, indets(f, 'name')[]); (4)
 > # # This (nameless) function can be applied to an # argument as in #   unapply( f, indets(f, 'name')[])(2);   unapply( f, indets(f, 'name')[])(a+b);  (5)
 > # # In a similar way, this (nameless) function can be # applied 'elementwise' to every element in a "container" # - ie a list, an array, vector etc, by using the '~' # operator. This will return a container of the same type, # but with updated entries. # # So for the list with entries 1,2,3 and check the returned # type #   help(elementwise);   unapply( f, indets(f, 'name')[])~([1,2,3]);   whattype(%);  (6)
 > # # And for the Array( 1,,2, 1..3, [[1,2,3],[4,5,6]]), # and check the returned #   unapply( f, indets(f, 'name')[])~(Array( 1..2, 1..3, [[ 1,2,3],[4,5,6]]));   whattype(%);  (7)
 >

which contains

1. an expression for the the sum f[i](x) for i=0..N
2. a plot of the expression in (1) above for x=0..5
3. a matrix containing values of x and sum f[i](x) for i=0..N, where x varies from 0..5 in steps of 0.1

see the attached

 > restart;
 > N := 4: # # Change lhs of the assignment from f(x) to F(x) # to avoid potential conflict arising from using # the same name in both indexed and unindexed # contexts # # Also changed sum() to add() #   F(x) :=  add(p^i*f[i](x), i = 0..N); # # Changed dependent variable throughout from f(x) # to F(x) #   HPMEq := (1 - p)*diff(F(x), x \$ 3) + p*(diff(F(x), x \$ 3) + 1/2*diff(F(x), x, x)*F(x)); # # Initialise sol as an empty list #   sol:=[]:   for i from 0 to N do       sol:= [ sol[],               dsolve               ( [ eval                   ( coeff(HPMEq, p, i) = 0,                     sol                   ),                   f[i](0) = 0,                   D(f[i])(0) = 0,                   D(f[i])(5) = 1                 ]               )             ];   end do:   sol;   (1)
 > # # Idle curiosity - what do these functions # look like #   plot( [ seq           ( rhs(sol[j]),             j=1..N+1           )         ],         x = 0..5,         legend = [ seq                    ( lhs(sol[j]),                      j=1..N+1                    )                  ]       ); > # # Expression for the sum of the functions f[i](x) # for i from 0..N #   fsum:= add          ( rhs(sol[j]),            j=1..N+1          ); # # And what does the plot of the sum of these # functions look like #   plot( fsum,         x=0..5       ); # # And generate a matrix whose columns are the # independent variable 'x' and sum(f[i](x)) for # x=0..5 in steps of 0.1 #   interface(rtablesize=60):   M:= Matrix( [ seq                 ( [ i, eval                        ( add                          ( rhs(sol[j]),                            j=1..N+1                          ),                          x=i                        )                   ],                   i=0..5, 0.1                 )               ]             );   interface(rtablesize=10):   (2)
 >
 >

## Hmm Maple 11...

Funny how often version numbers appear after I complain that they did not exist

Now I can't debug anything in Maple 11 ( released in 2007 ) because I only have the last seven Maple versions (going back to Maple 18 released in 2014). So the only thing I can suggest is to "break down" the offending command into "simple" steps, and see which one fails.

See the highlighted-green execution group in the attached -and report at which stage does it throw an error message in Maple 11. This worksheet (like the one I originally posted) works perfectly in all Maple version from Maple 18 onwards

 > restart; # # Assorted parameters #   T:= 0.1: beta:= 0.6: k:= 3.5:   A1:= (2*k-1)/(2*k-3):   A2:= 8*sqrt(2/Pi)*(beta-1)*k*GAMMA(k)/(3*GAMMA(k-0.5)*(2*k-3)**(3/2)):   A3:= (4*k**2-1)/(2*(2*k-3)**2):   M1:= 0.1+sqrt(T+(1/A1)): # # define ODES and ICs #   odes:= diff(U1(x),x)=-diff(phi(x),x)/(U1(x)-T/U1(x)),          diff(phi(x),x\$2)=(1+A1*phi(x)+A2*phi(x)**(3/2)+A3*phi(x)**2)-(M1/U1(x)):   bcs:= U1(0)=M1, phi(0)=0, D(phi)(0)=0.001: # # Solution and plots #   sol:= dsolve( [odes, bcs], numeric):   plots:-odeplot( sol,                   [ [x, U1(x)],                     [x, phi(x)],                     [x, diff(phi(x),x)]                   ],                   x=0..20,                   color= [red ,blue, green]                 ); > # # Some tests OP should apply # What does sol(1) return? # It should return a list of the values of x (=1) # U1(x), phi(x) and diff(phi(x),x), all evaluated at x=1 #   sol(1); # # so sol(1)[1..2] should return the first two elements of # the above list, as in #   sol(1)[1..2]; # # And rhs~( sol(1)[1..2]) ought to return the right hand # side of the expression in the above list, as in #   rhs~( sol(1)[1..2]); # # And all the seq() command does, is the above calculation # for values of x from 0..20, as in #   seq( rhs~(sol(j)[1..2]), j=0..20); # # And the Matrix() wrapper on this command inserts the data # into a Matrix wher the firsdt row has been explicitly inserted # as in x U1(x) #   Matrix( [ [ x, U1(x)],                   seq( rhs~(sol(j)[1..2]), j=0..20)                 ]               );     (1)
 > # # Generate the values of x and U1(x) for # x=0..20 in unit steps #   interface(rtablesize=25):   res:= Matrix( [ [ x, U1(x)],                   seq( rhs~(sol(j)[1..2]), j=0..20)                 ]               );   interface(rtablesize=10): # # Now you could use Maple's export() command # to output the Matrix 'res' in any number of # different formats - dependin on which format # the target application understands # # But why bother? # # After what kind of data processing is available # is available in the target application, which # isn't available in Maple? # # If I had to guess-the answer would be none at all # (2)
 > restart; # # Convert the whole of the above to a procedure # will accept any values of T, beta, k amd return # the desired plot #   getODESol:= proc( T, beta, k)                     local A1:= (2*k-1)/(2*k-3),                           A2:= 8*sqrt(2/Pi)*(beta-1)*k*GAMMA(k)/(3*GAMMA(k-0.5)*(2*k-3)**(3/2)),                           A3:= (4*k**2-1)/(2*(2*k-3)**2),                           M1:= 0.1+sqrt(T+(1/A1)),                         #                         # define ODES and ICs                         #                           odes:= [ diff(U1(x),x)=-diff(phi(x),x)/(U1(x)-T/U1(x)),                                    diff(phi(x),x\$2)=(1+A1*phi(x)+A2*phi(x)**(3/2)+A3*phi(x)**2)-(M1/U1(x))                                  ],                           bcs:= [ U1(0)=M1, phi(0)=0, D(phi)(0)=0.001],                         #                         # Solution and plots                         #                           sol:= dsolve( [odes[], bcs[]], numeric):                     return plots:-odeplot( sol,                                            [ [x, U1(x)],                                              [x, phi(x)],                                              [x, diff(phi(x),x)]                                            ],                                            x=0..20,                                            color= [red, blue, green]                                          );                end proc: # # Now plot all sorts of combinations - after all who # knows what the OP wants/needs #   plots:-display( [seq( getODESol( j, 0.6, 3.5), j=0.1..0.5,0.1)]);   plots:-display( [seq( getODESol( 0.1, j, 3.5), j=0.3..1.0,0.1)]);   plots:-display( [seq( getODESol( 0.1, 0.6, j), j=3..4,0.1)]); # # This one gernartes some singularity warnings which # I haven't bothered to nail down - and I won't until # the OP provides details of which combinations of # parameters (s)he is interested in #   plots:-display( [ seq                     ( seq                       ( seq                         ( getODESol( i, j, k),                           i=0.5..0.7, 0.1                         ),                         j=0.5..0.7, 0.1                       ),                       k=3.4..3.6, 0.1                     )                    ]                   );    >
 >

 1 2 3 4 5 6 7 Last Page 1 of 160