tomleslie

13876 Reputation

20 Badges

15 years, 174 days

MaplePrimes Activity


These are answers submitted by tomleslie

You can (nearly) always use showstat() to examine the procedure which Maple is running, so

restart;
with(inttrans):
showstat(invlaplace);

will provide a listing of the invlaplace() procedure.

For a specific calculation, you can play with infolevel(). Be aware that if you set this on a 'complicated' example (such as yours), you may get a lot of output. Probably safer (at least initially) to try with a simple example, as in

restart;
with(inttrans):
infolevel[all]:=3:
invlaplace(1/(s-a), s, t);;

Higher infolevels (up to 5) will provide even more.

 

 

carefully(!!) the following toy example

   restart;
   eq1:=diff(y(t),t);
#
# The following will perform a "syntactic'
# substitution - ie each occurrence of the
# operand 't' will be replaced, whether or
# not the expression which results is
# syntactically correct. In this case the
# substitution produces an expression which
# is syntactically incorrect.
#
   eq2:=subs(t=t/T, eq1);
#
# By default the subs command does not force
# an evaluation. This evaluation will be forced
# simply by calling the assignment again, in
# which case the error is manifest
#
   eq2;
#
# The error message does a good job of explaining
# the error. The second argument to a diff()
# command has to be a 'name' (or evaluate to a
# 'name') and t/T is not a name
#
# You may be able to achieve the substitution you
# desire by using subsop(), as in
#
   subsop(1=t/T,eq1);

 

When you write f(1/0), the argument 1/0 is evaluated - and throws the Error.

The "evaluation" process is two-step -automatic simplification followed by evaluation. Once can control the latter process (using uneval quotes, for example). However one cannot (AFAIK) intercept the automatic simplification process, and it seems to be this stage which is generating the error. My justification for this statement is that I cannot 'catch' a simple divide_by_zero error using a try-catch construct.

As you have observed the NumericEventHandler will allow you to redefine what happens on division by zero - but you only have one choice: so 1/0 and -1/0 will either both be infinity of -infinity whihc may no be what you want. You can define your own procedure for the NumericEventHandler (as you have done) but you will still not be able to get the correct sign the operation '1/0' is interpreted as 1*0^(-1), and when using the NumericEventHandler() the only part of the operation one hsa access to is the 0^(-1) part, so one cannot access the sign of the numerator

You cannot use the NumericEventHandler() within the procedure, because (as observed above) the arguments will be evaluated, and the Error thrown, before you enter the procedure.

My conclusion is that (since you don't like the "global" NumericEventHandler(), and in case it always seems to have a sign ambiguity), you basically have two choices

  1.  "trap" the situation 1/0, before invoking the procedure ought to be fairly trivial to do
  2. rewrite the procedure to accept an optional second argument, with the default value of the second argument being unity. This way, f(infinity) and f(1,0), will deliver the same result. As in

restart;
f:=proc(a, b:=1)
             local u, x;
             if       b=0
             then x:=csgn(a)*infinity
             else x:=a/b
             end if:
             Digits:=15;
             if        x>10^9
             then  u:=1/x;
                       evalf(Pi/2*arccos((u*u-1)/(1+u*u)));
            else    evalf(Pi/2*arccos((1-x*x)/(1+x*x)))
    end if
end proc:

This works for f(infinity) and f(1, 0). It returns undefined for f(-infinity) or f(-1,0), becuase your original procedure does not handle negative infinity

Corrected a typo in the above code - Sorry

I produced the attached to fit your data in a variety of different ways. Problems occurrd because

  1. Your data is incredibly "noisy" so you can fit almost any function you want
  2. The fact that the values of your independent variable ranges from 1988 to 2001 causes some numerical probelms. These can be circumvented simply by subtracting 1988 from the x-data and fitting the result. (One could put this back in after performing the fit -but I haven't bothered)

fits.mw

probaby ought to be Statistics:Fit

This will accept a "model function", containing parameters: a list/vector/whatever of independent variable values, a list/vector/whatever of corresponding dependent variable values, and a variety of output options.

Default output option is to return the model function with parameters set to values which provides the best overall "fir". If the correct output options are specified then you can get such things as the "residuals" at any of the values of the independent variables. (other output options are available).

If you suppied , independent variable values, corresponding dependent variable values, model function, etc then we could try it for you (and maybe make other suggestions, cos other possiblities are available) - but since you don't, we can't

Assuming I understand the requirement correctly, try

restart;
Y:= -L*epsilon^epsilon*(1-epsilon)^(1-epsilon)*(-omega^(-1+epsilon)*k*theta+(1-theta)*omega^epsilon)/(1-epsilon-theta);
z:=freeze(omega^epsilon):
Y1:=thaw(simplify(collect(subs( omega^epsilon=z, expand(Y)), z)));

 

 

Assuming that

  1. you only ever have to deal with square "Design" matrices
  2. you "know" the matrix form you are trying to achieve

then you can do the check on an element by element basis as follows


 

  restart;
  with(LinearAlgebra):
#
# Set X-matrix to be square of dimension Dim
#
  Dim:=3;
  X:=Matrix(Dim, Dim, symbol=x);
  Y:=Matrix(Dim, 1,   symbol=y);
  B:=Matrix(Dim, 1,   symbol=b);

#
# Form the estimator. Note that X.B-Y will result
# in a matrix which is Dim x 1. So the transpose
# will be 1 x Dim, so the estimator will end up as
# 1 x 1
#
  S:= Transpose(X.B-Y).(X.B-Y);

#
# Peform the differentiations on the estimator,
# solve for all b-values, and construct a Dim x 1
# katrix from the results
#
  ans1:= Matrix( Dim,
                 1,
                 rhs~( solve
                       ( [ seq( diff(S[1,1], b[j,1]),
                                j=1..Dim
                              )
                         ],
                         [ seq
                           ( b[j,1],
                             j=1..Dim
                           )
                         ]
                       )[]
                     )
               );

#
# Produce the explicit form of the "expected" answer
#
  ans2:=simplify~(MatrixInverse(Transpose(X).X).Transpose(X).Y);

#
# Compare elements in the "calculated" and  "expected"
# should return nothing but 'true
#
  seq( is(ans1[j,1]=ans2[j,1]), j=1..Dim);

 


 

Download estimator.mw

This works for Dim =2,3,4 although the MatrixInverse() operation is already getting a bit sluggish

You don't provide any data (ie depth, veloc) for me to perform this fit - makes things  a bit tricky

I generated my own data using x-values ('corresponding to 'depth' or u-values in your worksheet) 1..10, and then calculating the values of the dependent variable (veloc), by using your model function, with

  1. arbitrary value for the parameters a, b. I chose a=2.5, b=0.75, since these seem to be within the ranges you expect
  2. an additive random variable in the range -0.1..0.1, just to add some "noise" to the dependent variable value

I also fixed up a few minor syntax issues, which led to the attached.  Obviously your input data may be vastly different form what I have used, but the fitting process *ought* to be similar

dataFit.mw

If I solve the PDE in A(x,t) which is part of your "by-hand-real "solution, then the only real difference I see between this and the Maple solution is that the formerer is missing a simple additive arbitrary constant in the definition of the function nu().

pdetest() confirms that both the Maple solution and the "by-hand-real" solution are valid


 

eqs := {-(diff(tau(x, t, u), u)), -(diff(tau(x, t, u), x)), -(diff(xi(x, t, u), t)), -(diff(xi(x, t, u), u)), -(diff(tau(x, t, u), u, u)), -(diff(xi(x, t, u), u, u)), -(diff(tau(x, t, u), x, u))+diff(tau(x, t, u), x), -(diff(xi(x, t, u), u, t))+diff(xi(x, t, u), t), -(diff(eta(x, t, u), t))+diff(eta(x, t, u), u, t)+2*(diff(xi(x, t, u), u))-(diff(xi(x, t, u), x, t)), diff(xi(x, t, u), x)+diff(tau(x, t, u), t)-(diff(eta(x, t, u), u))+diff(eta(x, t, u), x, t), -(diff(tau(x, t, u), x, t))-(diff(eta(x, t, u), x))+2*(diff(tau(x, t, u), u))+diff(eta(x, t, u), x, u), -(diff(xi(x, t, u), x, u))-(diff(tau(x, t, u), u, t))-(diff(eta(x, t, u), u))+diff(eta(x, t, u), u, u)}
NULL

{-(diff(diff(tau(x, t, u), u), u)), -(diff(diff(xi(x, t, u), u), u)), -(diff(tau(x, t, u), u)), -(diff(tau(x, t, u), x)), -(diff(xi(x, t, u), t)), -(diff(xi(x, t, u), u)), -(diff(diff(tau(x, t, u), u), x))+diff(tau(x, t, u), x), -(diff(diff(xi(x, t, u), t), u))+diff(xi(x, t, u), t), -(diff(diff(tau(x, t, u), t), x))-(diff(eta(x, t, u), x))+2*(diff(tau(x, t, u), u))+diff(diff(eta(x, t, u), u), x), -(diff(diff(xi(x, t, u), u), x))-(diff(diff(tau(x, t, u), t), u))-(diff(eta(x, t, u), u))+diff(diff(eta(x, t, u), u), u), -(diff(eta(x, t, u), t))+diff(diff(eta(x, t, u), t), u)+2*(diff(xi(x, t, u), u))-(diff(diff(xi(x, t, u), t), x)), diff(xi(x, t, u), x)+diff(tau(x, t, u), t)-(diff(eta(x, t, u), u))+diff(diff(eta(x, t, u), t), x)}

(1)

 

maple_sol := pdsolve(eqs)

{eta(x, t, u) = _C4+_C5*exp(_c[1]*x)*_C6*exp(t/_c[1])*exp(u), tau(x, t, u) = -_C1*t+_C3, xi(x, t, u) = _C1*x+_C2}

(2)

real_sol := {diff(A(x, t), x, t)-A(x, t) = 0, eta(x, t, u) = A(x, t)*exp(u), tau(x, t, u) = -C2*t+C4, xi(x, t, u) = C2*x+C3}

{diff(diff(A(x, t), t), x)-A(x, t) = 0, eta(x, t, u) = A(x, t)*exp(u), tau(x, t, u) = -C2*t+C4, xi(x, t, u) = C2*x+C3}

(3)

NULL

#
# check the Maple solution
# NB - should return {0}
#
  pdetest( maple_sol, eqs);

{0}

(4)

#
# real sol is "incomplete" so solve the remaining
# pde which it contains
#
  real_sol2:=eval( real_sol, pdsolve(real_sol[1], build))[2..-1];
#
# Chech whether this is a solution
#
  pdetest( real_sol2, eqs);

{eta(x, t, u) = _C1*exp(_c[1]*x)*_C2*exp(t/_c[1])*exp(u), tau(x, t, u) = -C2*t+C4, xi(x, t, u) = C2*x+C3}

 

{0}

(5)

 

``


 

Download pdeSols.mw

that m is the length of the vectors, L0,u,v,w - then the attached wiil work

There really is no need to "subscript" all the variables in your loop using '||i' -  unless you wish to access these later. Doing this just makes your code 'cluttered' and difficult to read


 

restart;
A:=1:E:=0.1:m0g:=0.05:
L0:= <  5,    5,    5,    5 >:
u:=  <-5/4, -1/4, -3/4, -1/4>:
v:=  <-1/2, -3/2,  1/2,  1/2>:
w:=  <  0,    1,    0,   -1 >:
Force:=Matrix(op(1,L0),2):

for i from 1 to op(1,L0) do
    r1:=H*L0[i]/(A*E)+H/m0g*(arcsinh(V_A/H)+arcsinh((m0g*L0[i]-V_A)/H))-sqrt(u[i]^2+v[i]^2);
    r2:=V_A*L0[i]/(A*E)*(m0g*L0[i]/(2*V_A)-1)+H/m0g*(sqrt(1+((m0g*L0[i]-V_A)/H)^2)-sqrt(1+(V_A/H)^2))-abs(w[i]);
    Force[i,..]:=convert(rhs~(fsolve({r1, r2}, {H,V_A}, H=0..1)),Array);
end do:
Force;

_rtable[18446744074333389278]

(1)

 


 

Download eqSolve.mw

as in

f:=unapply(rsolve({s(1)=1, s(2)=2, s(n)=s(n-1)+6*s(n-2)}, s),n);
seq(f(n), n=1..30);

 

Not the greatest fit I've ever seen!


 

restart;
with(Statistics):
with(plots):

Xvals := [.53993447, .5599647, .57995479, .59995566, .61996118, .63994512, .65994136, .6799731, .69996782, .71997949, .73997422, .75995044, .77994976, .7999707, .81995146, .83995244, .85996729, .87994951, .8999458]:
Yvals := [-.79625455, -.75585259, -.67800183, -.47955884, -.25493698, -0.65747361e-1, .10114507, .26994542, .41484068, .50621122, .60363251, .65510417, .7251356, .75804148, .76002419, .8030069, .82774732, .8429559, .8692888]:

points := ScatterPlot(Xvals, Yvals, symbol=solidcircle, symbolsize=16):
fit:=Statistics:-Fit(a*y^3+b*x, Xvals, Yvals, y);
pfit:=plot(fit):
display([points, pfit], view=[0..1, -1..1]);

HFloat(2.996866026020011)*y^3-HFloat(0.9272751314432008)

 

 

 


 

Download fitdata.mw

is illustrated in the attached, which uses the standard example that the orthogonal trajectories of parabolas are ellipses.

orthoTraj.mw

I can't find your previous question either, so I am reposting the same code which I uploaded last time

anotherDSissue.mw

As I pointed out in my previous response, your problem is that for some values of the parameter 'uc', the fsolve() command fails (ie returns unevaluated), at which point all subsequent code breaks. You can circumvent this by using the SolveEquations() command from the DirectSearch package. This works for all 'uc' values which I tried.

The DirectSearch package can be downloaded from the Maple Application Centre

 

your previous question at

http://www.mapleprimes.com/questions/222319-How-Can-I-Draw-Odeplot-This-Numeric-Dsolve

So the answer is a slight variation of the one I posted there - see attached

ODEplot4.mw.

NB You will have to change the pathName in the ImportMatrix() command to somethng appropraite for your machine

First 161 162 163 164 165 166 167 Last Page 163 of 207