Carl Love

Carl Love

28020 Reputation

25 Badges

12 years, 300 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

Addressing the original question, how to collect all terms in a polynomial containing a certain variable without using select: Let p be the polynomial and v the variable. Then

p:= collect(p, v);
p - tcoeff(p,v);

See ?tcoeff (trailing coefficient). But I still don't understand why you don't want to use select.

I made up some fake data as an example, since I couldn't download your actual data. You have an interesting situation in that your model is given in implicit form, and you use zeros for the dependent values.

My plots are very boring because of the fake data, and my fit is very poor. I'd like to see it with the actual data.


(**)

restart:

Note that I use exact rational values in the model.

(**)

Model:= i0*(exp(1000*(v-i*rs)/n0/(259/100))-1)-i;

i0*(exp((100000/259)*(-i*rs+v)/n0)-1)-i

It turns out that Maple can solve the model for i in terms of v. This will make the plotting easier.

(**)

NewModel:= solve(Model, i);

(1/100000)*(259*n0*LambertW((100000/259)*i0*rs*exp((100000/259)*(i0*rs+v)/n0)/n0)-100000*i0*rs)/rs

(**)

dim:= 100:

(**)

V:= [.01*k $ k= 1..dim]:

(**)

idata:= map(exp, V):

(**)

Data:= < < V[] > | < idata[] > >:

I omitted the paramterranges because they wouldn't apply to my fake data, but should definitely include them

(**)

Fit:= Statistics:-NonlinearFit(
     Model, Data, Vector(dim), [v,i],
     #parameterranges= [ ],
     output= [residualsumofsquares, parametervalues]
);

[24.4439180403656, [i0 = HFloat(-1.7268875565927129), n0 = HFloat(1.0), rs = HFloat(1.0000000210837114)]]

(**)

Curve:= plot(eval(NewModel, Fit[2]), v= 0..1):

(**)

Points:= plot(Data, style= point, symbol= cross, color= green):

(**)

plots:-display([Curve, Points], axis[2]= [mode= log]);

(**)

 


Download Nonlinear_Fit.mw

 

 

Assuming that your goal is mainly to get the plots, as you say, and that the Matrix was only an intermediate step towards that, there is no need for the actual Matrix. This does the job:

 

 


(**)

restart:

(**)

de1a:= diff(y(eq), eq) = y(eq):

(**)

de1b:= diff(y(eq), eq) = 2*y(eq):

(**)

A:= dsolve(

     {de1a, y(0)=1}, {y(eq)}, numeric,

     events= [[y(eq)-2, halt]]

):

(**)

ic2:= A(1);

Warning, cannot evaluate the solution further right of .69314728, event #1 triggered a halt

[eq = HFloat(0.6931472897592507), y(eq) = HFloat(1.9999999999999998)]

(**)

B:= dsolve(

     {de1b, y(eval(eq, ic2))= eval(y(eq), ic2)}, {y(eq)},

     numeric

):

(**)

p1:= plots:-odeplot(

     A, 0..eval(eq, ic2), color= red, legend= "Elastic"

):

(**)

pt:= plot(

     [eval([eq, y(eq)], ic2)], style= point, symbolsize= 32

):

(**)

p2:= plots:-odeplot(

     B, eval(eq, ic2)..1, color= green, legend= "Plastic"

):

(**)

plots:-display([p1,pt,p2]);

(**)

 



Download splitting_event.mw



And if you still want to put the data into a Matrix, then the evaluation procedures A and B are still available for that.

I get the same sequence each time. For example, this code produces no output, indicating that the same sequence of 2^6 numbers was generated 2^8 times.

with(RandomTools):
randomize(1):
R:= 'round( Generate(distribution(DiscreteUniform(5,10))) )' $ 2^6:
for n to 2^8 do
     randomize(1);
     if R <> ('round( Generate(distribution(DiscreteUniform(5,10))) )' $ 2^6) then
          print(n)
     end if
end do:

So I guess that you should post your code.

Although undocumented, the code is visible by showstat, and is broken into easy-to-understand chunks. Start with

showstat(`convert/piecewise`);
showstat(PiecewiseTools:-Convert);
showstat(PiecewiseTools:-Import);
showstat(PiecewiseTools:-ImportImplementation:-StandardFunctions:-MAX);
etc.

Your code is quite bloated. It can be reduced to the following:

Chisqr1:= proc(w)
local
     j,
     x:= [seq(4*j, j= 0..12)],
     y:= [16759.00586, 19630.75977, 15190.14551, 11532.44141, 6447.006836,
          7827.111816, 3531.188965, 5167.17334, 10127.59277, 11626.91504,
          5197.712402, 6539.370605, 6164.602539
          ],
     A, B, C, D,
     S:= add((y[j]-A-x[j]*B-C*cos(w*x[j])-D*sin(w*x[j]))^2, j= 1..13),
     X:= VectorCalculus:-Gradient(S, [A, B, C, D]),
     M1:= < seq(`<|>`(map2(coeff, X[j], [A, B, C, D])[]), j= 1..4) >,
     M2:= LinearAlgebra:-LinearSolve(M1, map(tcoeff, X))
;
     eval(S, [A,B,C,D] =~ convert(M2, 'list'))    
end proc;

You asked:

Does anyone know why, if you plot a discontinuous function of one variable, Maple's default behavior is to try to make the function continuous by adding a vertical line segment?

It is not trying to make it continuous; it is assuming that it is continuous. Detecting discontinuities involves significant symbolic processing. The default is to plug values into the expression and evaluate numerically, which is much faster.

I don't know why anyone would anyone would ever want to graph a function this way.

It is not assumed that anyone would want it. The behaviour is a side effect of the numerical technique; it is not intended behaviour.

I couldn't find a place on Maple's Web site to submit a complaint (this is not really a "support question").

From the row of tabs near the top of this (MaplePrimes) page, pull down "More" and select "Submit Software Change Request".

If Maple is installed on a computer, can I change some settings so I don't have to do this?

Yes, it is easy. Add the following line to your initialization file maple.ini:

local plot:= proc({discont:= true}) :-plot(args, :-discont= discont) end proc:

Then the default for discont will be true, and you'll need to specify discont= false on individual plots for which you want to turn off discontinuity checking. (Also see the other suboptions for discont at ?plot,discont .)

Related question: if you plot a discontinuous function of two variables, such as Arg(x,y) (the argument of x + yi; I apologize if this is not the exact syntax)

I am assuming that you mean plot3d(argument(x+I*y), x= -1..1, y= -1..1).

you may see a vertical "cliff" which does not belong in the graph.  I do not blame the Maple people for this - it seems like a difficult issue.  Does anyone know a good way to fix such plots?

Yes, make it a parametric plot in cylindrical coordinates.

plot3d([r, theta, theta], theta= -Pi..Pi, r= 0..1, coords= cylindrical);

That is surprising.

Here's a workaround: Create a calendar by Joining the desired calendar with the Null calendar:

c:= JoinHolidays(Null, Toronto);

Your system of two difference equations can be solved by mapping ztrans and invztrans over the system:


(**)

Eqns:= {y(n+1) + x(n) = 1+2^(n+1), y(n) + 2*x(n+1) = 2 + 2^n}:

(**)

ICs:= {x(0)=1, y(0)=1}:

(**)

eval(map(ztrans, Eqns, n, z), ICs);

{z*ztrans(y(n), n, z)-z+ztrans(x(n), n, z) = z/(z-1)+z/((1/2)*z-1), ztrans(y(n), n, z)+2*z*ztrans(x(n), n, z)-2*z = 2*z/(z-1)+(1/2)*z/((1/2)*z-1)}

(**)

solve(%, map(ztrans, {y(n), x(n)}, n, z));

{ztrans(x(n), n, z) = z/(z-1), ztrans(y(n), n, z) = z/(z-2)}

(**)

map(invztrans, %, z, n);

{x(n) = 1, y(n) = 2^n}

Or it can all be done with a single call to rsolve:

(**)

rsolve({y(n+1) + x(n) = 1+2^(n+1), y(n)+2*x(n+1)=2+2^n, x(0)=1, y(0)=1}, {y(n), x(n)});

{x(n) = 1, y(n) = 2^n}

(**)

 


Download ztrans.mw

From reading your attached file, I think what you want is

for k from 0 to 1 do
     for m from 0 to 4 do
          C[m,k]:= coeff(Eq, Y[k], m) = 0
     end do
end do;

But also see my Reply to your previous version of this Question.

Just take the ln of both sides of each equation, and you then have a simple system of linear equations.


sys:= {2^x*3^y*5^z = 240, 3^x*2^y*5^z = 360,  5^x*3^y*2^z = 1350}, {x,y,z}:

(**)

solve(sys):

(**)

simplify(%);

{x = 2*(ln(5)^2+2*ln(5)*ln(3)-2*ln(2)^2-ln(2)*ln(3))/(ln(5)^2+ln(5)*ln(3)-ln(2)^2-ln(2)*ln(3)), y = (ln(5)^2+3*ln(5)*ln(3)-3*ln(2)^2-ln(2)*ln(3))/(ln(5)^2+ln(5)*ln(3)-ln(2)^2-ln(2)*ln(3)), z = (ln(5)^2+2*ln(5)*ln(2)+ln(5)*ln(3)-ln(2)^2-ln(2)*ln(3)-2*ln(3)^2)/(ln(5)^2+ln(5)*ln(3)-ln(2)^2-ln(2)*ln(3))}

(**)

evalf(%);

{x = 2.82637650486254, y = 1.82637650486254, z = .941362406163020}

(**)

 


Download lin_sys.mw

It seems like the seventh and eighth rows of your Matrix are being converted to lists by CodeGeneration:-Matlab. I have no idea why this is (there are no lists in your input Matrix), but here is a workaround: Change your last command to

Matlab(
     subsop([-1,1]= J, eval([codegen:-optimize](tmp, tryhard), pow= `^`)),
     output = string, defaulttype = numeric
);

I don't have Matlab to test this, so please let me know how this works for you.

y(seq(x[k], k= 1..n))

If x is actually a list, you can do y(x[]).

If you make assumptions such as assume(m>0, m<1), assume(m>2, m<3), etc., that place m between two consecutive integers, then MultiSeries:-multiseries will work.

You need to use one section per slide. So, if you have a section that is too long to fit on one slide, then you need to divide it into multiple sections.

First 351 352 353 354 355 356 357 Last Page 353 of 395