tomleslie

In principle...

Answered: tomleslie 13876

May 25 2018

0 4

the best you are probably going to be able to do is to use

solve({eq1, eq3, eq5}, {a, b, c});

which will attempt to solve for the variables {a, b, c} of the equations {eq1, eq3, eq5} in terms of the 'other' unknowns. This may or may not be possible.

fsolve() will be of no use in an underdetermined system, since it only works for cases where a purely numerical solution is possible.

Further assistance/suggestions may be possible for your specific problem. This will almost certainly require you to upload your worksheet using the big green up-arrow in the Mapleprimes toolbar.

Badly placed restart;...

Answered: tomleslie 13876

May 25 2018

1 0

with(LinearAlgebra) is used to load the LinearAlgebra() package, in order to access the Determinant() command. However you subsequently issue a 'restart' command, wihc will unload this package, meaning that the Determinant() command is undefined - this helps the solve() command to fail.

If I just swap round these two commands, then your code will execute and produce solutions - although all of those I got were of the form a+I*b, where both 'a' and 'b' are very close to zero - so much so that maybe Determinant(A)=0, only when omega=0.

Anyhow you can investigate this at your leisure

>

restart*with(LinearAlgebra); delta := 0.407e-9; C := 2023.2; L := 5*delta

(1)

>

eq := (1/12)*delta^2*S^4+S^2+omega^2/C^2

(2)

>

(3)

>

U := _C1*exp(sqrt(-2*C*(3*C-sqrt(-3*delta^2*omega^2+9*C^2)))*x/(C*delta))+_C2*exp(-sqrt(-2*C*(3*C-sqrt(-3*delta^2*omega^2+9*C^2)))*x/(C*delta))+_C3*exp(sqrt(-2*C*(3*C+sqrt(-3*delta^2*omega^2+9*C^2)))*x/(C*delta))+_C4*exp(-sqrt(-2*C*(3*C+sqrt(-3*delta^2*omega^2+9*C^2)))*x/(C*delta))

>

(4)

>

Matrix(%id = 18446744074835891374)

(5)

>

Warning, solutions may have been lost

(6)

>

Download detProb.mw

"premature evaluation", rather than "sco...

Answered: tomleslie 13876

May 24 2018

1 0

In a plot command (like many others), the arguments are evaluated before the plot is attempted.

Forget procedures for the moment, thta's just an extral layer of complication: at the top level

x:=10;
plot(sin(x),x=-Pi..Pi);

The arguments will evaluate to sin(10), and 10=-Pi..Pi : sin(10) won't be what you want, but is valid! However 10=-Pi..Pi is not a valid option fpr the plot command.

You can circumvent this problem, by using unevaluation quotes. For example, in

x:=10;
plot(sin(x),'x=-Pi..Pi');

the arguments will evaluate to sin(10) and x=-Pi..Pi, which is valid, but will plot a constant value, which again is probably not what you want!

Unevaluating both arguments, will provide what you (probably) want, ie

x:=10;
plot('sin(x)','x=-Pi..Pi');

I've just seen vv's answer. It is true in the above that you don't have to unevaluate the complete argument(s) - just so long as you unevaluate x - so rather than 'x=-Pi..Pi',the version 'x'=-Pi..Pi will work equally well. Ditto sin('x') rather than 'sin(x)'

Actually seems a bit of a waste...

Answered: tomleslie 13876

May 23 2018

0 1

not to post my version of the solution so FWIW.

Same caveat as Rouben: OP is using Maple 13, and I haven't got anything earlier than Maple 18 - si I might have inadvertently introduced some kind of version issue. Unlike Rouben I' didn't update the with(linalg), and associated commands to with(LinearAlgebra), even although the former has been deprectaed since about Maple V

>	restart: with(linalg): #deprecated even in Maple 13 with(plots):

>

  xo(t):=u*x(t):
  yo(t):=u*y(t):
  xm(t):=(1-u)*x(t):
  ym(t):=(1-u)*y(t):
  W:=Vector([p(t),q(t),r(t)]):
  Ibody:=Matrix( [ [A+M*yo(t)^2,-M*xo(t)*yo(t),0],
                   [-M*xo(t)*yo(t),B+M*xo(t)^2,0],
                   [0,0,C+M*(xo(t)^2+yo(t)^2)]]):
  Kbody:=multiply(Ibody,W):
  Iem:=Matrix( [ [m*ym(t)^2,-m*xm(t)*ym(t),0],
                 [-m*xm(t)*ym(t),m*xm(t)^2,0],
                 [0,0,m*(xm(t)^2+ym(t)^2)]]):
  Kem:=multiply(Iem,W):
  Kbm:=evalm(Kbody+Kem):
  dK1:=Vector([diff(Kbm[1],t),diff(Kbm[2],t),diff(Kbm[3],t)]):
  dK2:=crossprod(W,Kbm):
  dK:=evalm(dK1+dK2):
  dK[1]:
  dK[2]:
  dK[3]:
  x(t):=cxp*p(t)+cxq*q(t)+cxr*r(t)+cx0:
  y(t):=cyp*p(t)+cyq*q(t)+cyr*r(t)+cy0:
  pars:= [ A=8, B=6, C=4, M=60, m=6, u=0.09, P=1, Mz=1,
           cxp=0.00025, cxq=0.0625, cxr=-0.0375, cx0=0,
           cyp=0.075, cyq=0.00025, cyr=0.0375, cy0=0
          ]:
  Sys:=dK[1]+P*u*y(t)=0,
       dK[2]-P*u*x(t)=0,dK[3]-Mz=0,
       diff(theta(t),t)-p(t)*cos(phi(t))+q(t)*sin(phi(t))=0,
       diff(psi(t),t)-(p(t)*sin(phi(t))+q(t)*cos(phi(t)))/sin(theta(t))=0,
       diff(phi(t),t)-r(t)+cot(theta(t))*(p(t)*sin(phi(t))+q(t)*cos(phi(t)))=0:
  inc:=p(0)=1,q(0)=0,r(0)=0,theta(0)=1,psi(0)=2,phi(0)=1:
  SOL:=dsolve(eval( {Sys,inc},pars), numeric, maxfun=0):

>

  odeplot(SOL,[t,p(t)],0..100);
  odeplot(SOL,[t,q(t)],0..100);
  odeplot(SOL,[t,r(t)],0..100);
  odeplot(SOL,[p(t),q(t),r(t)],0..3500,numpoints=50000);
  odeplot(SOL,[t,eval(cxp*p(t)+cxq*q(t)+cxr*r(t)+cx0, pars)],0..700);
  odeplot(SOL,[t,eval(cyp*p(t)+cyq*q(t)+cyr*r(t)+cy0, pars)],0..700);
  odeplot(SOL, [t,phi(t)], t=0..100);
  odeplot(SOL, [t,theta(t)], t=0..100);
  odeplot(SOL, [t,psi(t)], t=0..100);

>

Download fixODE.mws

An alternative,...

Answered: tomleslie 13876

May 23 2018

2 0

just because alternatives are good.

Consider using a (simple) module, which offers both encapsulation and the speed you require, because all of the commands which only need to be performed once, will only be performed once, during the module definition.

Contemplate the attached.

>	restart; with(CodeTools): # # Check the length of time to run the ifactor() command # Usage( op([1,1], ifactor(45!+1)));

memory used=139.67MiB, alloc change=72.00MiB, cpu time=1.65s, real time=1.65s, gc time=62.40ms

(1)

>

  restart;
  with(CodeTools):
#
# Check the length of time to define a module containing
# the ifactor command, and a 'simple' procedure which uses
# the result of the ifactor() command. This is (almost)
# exactly the same as the above, because the ifactor()
# command is executed and output assigned to 'myfun' during
# the module definition
#
  m3:= Usage( module( )
                     local myfun;
                     export foo;
                     myfun:= op([1,1], ifactor(45!+1)):
                     foo:= proc(n)
                                n*myfun;
                     end proc;
              end module
            );

memory used=139.78MiB, alloc change=72.00MiB, cpu time=1.65s, real time=1.65s, gc time=78.00ms

(2)

>	# # Verify that subsequent calls to the procedure 'foo' do # not involve reevaluation of 'myfun' - so this should be # blindingly fast # Usage(m3:-foo(2), iterations=10000);

memory used=32 bytes, alloc change=0 bytes, cpu time=1.50us, real time=7.0e+02ns, gc time=0ns

(3)

>

#
# So for the test case of your integral, you would want
# to define a module something like
#
  restart;
  with(CodeTools):
  m1:= module( )
              local myfun;
              export foo;
              myfun:= unapply( int(x*sin(a*x),x=0..Pi), a) assuming a::posint;
              foo:= proc(n::posint)
                         myfun(n)
                    end proc;
       end module:
#
# Or even,
#
  m2:= module( )
              local myfun;
              export foo;
              myfun:= int(x*sin(a*x),x=0..Pi) assuming a::posint;
              foo:= proc(n::posint)
                         eval( myfun, a=n);
                    end proc;
      end module:
#
# Rather than
#
  foo:=proc(n::posint)
            local r;
            r:= int(x*sin(n*x),x=0..Pi);
  end proc:
#
# Check the speeds for these three options: using either
# of the 'modules' defined above is waaaay faster cos the
# integral was evaluated (just the once), during the
# module definition.
#
  Usage(m1:-foo(2), iterations=10000);
  Usage(m2:-foo(2), iterations=10000);
  Usage(foo(2), iterations=10000);

memory used=208 bytes, alloc change=0 bytes, cpu time=3.10us, real time=3.00us, gc time=0ns

memory used=312 bytes, alloc change=0 bytes, cpu time=3.10us, real time=3.30us, gc time=0ns

memory used=6.91KiB, alloc change=33.00MiB, cpu time=81.10us, real time=79.90us, gc time=4.68us

(4)

>

Download useModule.mw

limit versus Limit...

Answered: tomleslie 13876

May 21 2018

0 5

your procedure for g() uses Limit() rather than limit(). The former is the inert form of the limit command. If you change this to limit() you will get a very different answer!

err..... Joseph Lagrange (maybe)...

Answered: tomleslie 13876

May 21 2018

0 1

The Wikipedia article https://en.wikipedia.org/wiki/Lagrange_polynomial is quite informative

I think the basic idea that you could fit (n+1)-points with an n-th order polynomial was known before Lagrange.

Maybe he just got the credit because he showed that the process always worked, by constructing the interpolation polynomial in a specific way?

Dataframes...

Answered: tomleslie 13876

May 20 2018

0 3

when you Import() a 'csv' file as a dataframe, then default 'labels' will be added for each row and column - generally these are just 1..m and 1..n. These 'labels' will appear along the top and down the left edge of the displayed data structure (assuming the overall size is small enough to be displayable).

It is important to realise that DataFrames are not matrixes or Arrays, but are in fact implemented as modules, with a bunch of associated 'methods'. Thus the 'labels' are not really zeroth rows/columns in an Array.

If you enter A[2], then you will get all the row labels, as well as the data in the column whose label is 2

Observations...

Answered: tomleslie 13876

May 20 2018

0 1

According to my simple-minded calibrated eyeball, you need three boundary conditions on u(x,y) + three boundary conditions on theta(x,y) + one boundary condtion on v(x,y): ie a total of seven. You supply four. Hence there is no way that you are going to get a numeric solution (or even an analytic solution without arbitrary constants).

This, for the moment, ignores the fact that not only do you need seven BCs/ICs for a numeric solution, they have to be of the correct type - eg derivative boundary conditions cannot be 'normal' to the boundary

So lets just forget BCs/ICs altogether: The command pdsolve(sys) will return a couple of possible solutions, in terms of (some) arbitrary functions and constants - and these are worth some thought.

In the first solution u(x,y)=0, which isn't very 'interesting'.
In the second solution u(x,y)=C, and on examination of the boundary conditions which you do have, 'C' would have to be zero: so still not very interesting

I can only suggest that you double-check the PDEs which you have supplied. If these are OK, then come up with seven BCs/ICs whihc do not involve tangential derivatives

Well...

Answered: tomleslie 13876

May 20 2018

1 4

I'm not even going to seriously try to "improve" the code you supplied. Just a few comments

Every time you issue the divisors(n) command, then this will be re-evaluated. How many times will your two nested seq() commands re-evaluate divisors(n)??
Why do you do this unnecessary re-evaluation?
One problem with seq() commands is that you can't "break" when a solution meets a particular criterion

You state that your immediate requirement is to "obtain the least divisor of a number that is greater than it's square root".The attached provides two methods, both of which are reasonably quick. The second, using the select() function is marginally quicker and more memory-efficient.

>

restart;
getDiv1:=proc( n::posint)
               uses numtheory:
               local sqrtN:= evalf(sqrt(n)),
                     p:= divisors(n),
                     j:
               for j in p do
                   if   j > sqrtN
                   then break:
                   fi
               od;
               return j;
        end proc:
getDiv2:= proc( n::posint)
               uses numtheory:
               local sqrtN:= evalf(sqrt(n)),
                     p:= divisors(n):
               return min(select( i->i>sqrtN, p));
         end proc:
CodeTools:-Usage(getDiv1(20000000));
CodeTools:-Usage(getDiv2(20000000));

memory used=93.20KiB, alloc change=0 bytes, cpu time=0ns, real time=7.00ms, gc time=0ns

memory used=11.21KiB, alloc change=0 bytes, cpu time=0ns, real time=0ns, gc time=0ns

(1)

>

Download getDiv.mw

It would be much easier...

Answered: tomleslie 13876

May 17 2018

0 32

if you posted the relevant code using the biog green up-arrow in the MaplePrimes toolbar, but you need to set up a function/procedure which does the following

accepts as input values for the parameters p1, p2, ..p6
generates the vector V1 which is a function of p1, p2, .....p6
subtracts the vector V1 from your 'experimental' vector V2, ie compute V2-V1
square and add all elements of V2-V1, this will be a single number to be returned by the function/procedure.
Use the Optimization:-Minimize() command to find the values of p1,,p6 which minimizes the value returned by (4) above

Doesn't sound too difficult!

Syntax and more!...

Answered: tomleslie 13876

May 17 2018

0 15

When you are fitting a multivariate 'model function' (in your case with the variables 'x' and 't') the points to be fitted have to be supplied in a matrix. If your function has n variables, this matrix will have n+1 columns. The first n columns are values for the variables and the final column is the associated value of the function. So for your case, each row in this matrix has to be in the form

x-value, t-value, p1(x-value, t-value)

I have generated this matrix for the 21 x-values and 11 t-values which you specified: so the resulting matrix is 231x3.

I also had to change the syntax of the NonLinearFit() command to use this matrix

You were inconsistent in your definition of the parameters to be fitted - sometimes these were given as a[1], a[2], b[1], b[2] and sometimes a_1, a_2, b_1, b_2. I have used the latter throughout

I probably fixed a few other syntax errors as and when they occurred. The attached now produces the 'correct' answer

>

  restart;
  with(Statistics):
#
# Define function
#
  p1 := proc (x, t)
              return (sin(x)+sin(x)*a_1*x*t-sin(x)*t+cos(x)*sin(x)*b_1*t^2*x+(1/2)*sin(x)*t^2-sin(x)*((1/2)*cos(x)*a_1*x*t^2+(1/2)*sin(x)*a_1*t^2-(1/2)*cos(x)*t^2)-cos(x)*((1/2)*sin(x)*b_1*x*t^2-(1/2)*sin(x)*t^2)+cos(x)*((1/2)*sin(x)*a_1*x*t^2-(1/2)*sin(x)*t^2)-sin(x)*((1/2)*cos(x)*b_1*x*t^2+(1/2)*sin(x)*b_1*t^2-(1/2)*cos(x)*t^2)+a_1*x*((1/2)*sin(x)*a_1*x*t^2-(1/2)*sin(x)*t^2)-cos(x)*sin(x)*t^2+sin(x)*a_2*t^2*x^2+((1/2)*sin(x)*a_1*x-(1/2)*sin(x))*a_1*t^2*x-sin(x)*a_1^2*x^2*t^2+(1/2)*sin(x)^2*a_1*t^2+(1/2)*sin(x)^2*b_1*t^2)/(1+a_1*x*t+a_2*x^2*t^2);
        end proc:
#
# Set some values for the 'test' parameters.
# These will become the 'unknowns' which the
# nonLinearFit() command will try to determine
#
# NB p1() does not depend on b_2
#
  a_1 := 1: a_2 := 5: b_1:= 6: b_2 := 3:
#
# Generate a matrix of points which will be used
# for the fit. This matrix is organised with three
# columns
#
# column1 = x_value
# column2 = t_value
# column3 = p1( x_value, t_value)
#
  p1Vals := Matrix([seq( seq( [i, j, evalf(p1(i,j))], i=0..20), j=0..10)]);
#
# Produce a plot of the function with the assumed values
# of the parameters, for comparison with the yet-to-be-
# obtained plot of the function with the fitted parameters
#
  plt1:=plot3d( p1, 0..20, 0..10, color=red);
#
# Reset the parameters to be unknown
#
  a_1 := 'a_1': a_2 := 'a_2': b_1:= 'b_1': b_2 := 'b_2':
#
# Generate the values for these parameters from the
# NonLinearFit() command. Note that b_2 will not be
# obtained because it does not exist in p1()
#
  sol:=NonlinearFit(p1(x,t), p1Vals, [x, t], output = parametervalues);
#
# Note that a_1 and a_2 are very close to the 'assumed'
# values, but b_1 seems to be off by a mile!
#
# Produce a plot of the target function with these fitted parameters
# Looks pretty similar to the original plot - even with a grossly
# different value for b__1
#
  plt2:=plot3d( 'eval'(p1(x,t),sol), x=0..20, t=0..10, color=blue);
#
# Generate a few values of p1(xt) with [a_1=1, a_2=5, b_1=6],
# ie the original assumed values, and then with [a_1=1, a_2=5,
# b_1=-154326.052324557] ie the fitted values: just to see
# if the difference is significant - and it doesn't seem to be.
# I guess the original function just isn't that dependent on
# b_1!
#
  evalf(eval(p1(5,5), [a_1=1.0, a_2=5.0, b_1=6.0])); #assumed
  evalf(eval(p1(5,5), sol)); #fitted

  evalf(eval(p1(10,3), [a_1=1.0, a_2=5.0, b_1=6.0])); #assumed
  evalf(eval(p1(10,3), sol)); #fitted

  evalf(eval(p1(7,2), [a_1=1.0, a_2=5.0, b_1=6.0])); #assumed
  evalf(eval(p1(7,2), sol)); #fitted
#
# The reason why there is little/no dependence on the b_1
# parameter can be seen if one just applies a simplify()
# command to the expression with the procedure p1(). The
# simplified expression does not contain b_1!!!
#
  expr:=(sin(x)+sin(x)*a_1*x*t-sin(x)*t+cos(x)*sin(x)*b_1*t^2*x+(1/2)*sin(x)*t^2-sin(x)*((1/2)*cos(x)*a_1*x*t^2+(1/2)*sin(x)*a_1*t^2-(1/2)*cos(x)*t^2)-cos(x)*((1/2)*sin(x)*b_1*x*t^2-(1/2)*sin(x)*t^2)+cos(x)*((1/2)*sin(x)*a_1*x*t^2-(1/2)*sin(x)*t^2)-sin(x)*((1/2)*cos(x)*b_1*x*t^2+(1/2)*sin(x)*b_1*t^2-(1/2)*cos(x)*t^2)+a_1*x*((1/2)*sin(x)*a_1*x*t^2-(1/2)*sin(x)*t^2)-cos(x)*sin(x)*t^2+sin(x)*a_2*t^2*x^2+((1/2)*sin(x)*a_1*x-(1/2)*sin(x))*a_1*t^2*x-sin(x)*a_1^2*x^2*t^2+(1/2)*sin(x)^2*a_1*t^2+(1/2)*sin(x)^2*b_1*t^2)/(1+a_1*x*t+a_2*x^2*t^2);
  simplify(expr);

(2+(2*a_2*x^2-2*a_1*x+1)*t^2+(2*a_1*x-2)*t)*sin(x)/(2*a_2*t^2*x^2+2*a_1*t*x+2)

(1)

>

Download fitProb.mw

Post code not pictures...

Answered: tomleslie 13876

May 16 2018

0 7

Nobody can do anything with a picture.

Do you mean something like

restart;
M1:=Matrix(6,6,(i,j)->cos((i-1)*(j-1)*2*Pi/6)-I*sin((i-1)*(j-1)*2*Pi/6));
M2:=Matrix(6,6,(i,j)->cos((i-1)*(j-1)*2*Pi/6)+I*sin((i-1)*(j-1)*2*Pi/6));
evalc~(M1.M2);

Scanning order for list...

Answered: tomleslie 13876

May 16 2018

1 10

Compare the output of the final two commands in order to see the 'scanning' order for the array as used by Maple's convert/list command

Download arrScan.mw

Syntax problems...

Answered: tomleslie 13876

May 16 2018

0 0

Your first example fails because

the second argument is a list, not a vector
even if the second argument were a vector, it would fail because the arguments are in the wrong order

Your second example fails because there is no command 'multtiply' - but there is a command 'Multiply' - note capitalisation.

The following 'toy' examples all work correctly

with(LinearAlgebra):
M := RandomMatrix(4,4);
c:=Vector[row](4,[0,0,0,1]);
d:=Vector[column](4,[0,0,0,1]);

VectorMatrixMultiply(c, M^(-1));
MatrixVectorMultiply(M^(-1), d);
Multiply(c, M^(-1));
Multiply( M^(-1), d);

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

13876 Reputation

20 Badges

MaplePrimes Activity

These are answers submitted by tomleslie

In principle...

Badly placed restart;...

"premature evaluation", rather than "sco...

Actually seems a bit of a waste...

An alternative,...

limit versus Limit...

err..... Joseph Lagrange (maybe)...

Dataframes...

Observations...

Well...

It would be much easier...

Syntax and more!...

Post code not pictures...

Scanning order for list...

Syntax problems...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again