Carl Love

Carl Love

26488 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

There are several issues.

  1. Remove the inner brackets. If brackets were needed for mathematical grouping (they're not in this case), they would have to be round parentheses. The outer brackets are fine---they indicate that you're building a list.
  2. You need to give the command evalf(%). This is the command that reduces things, as much as possible, to decimal values.
  3. I can't tell how you entered Pi. If you typed pi, you need to change it to Pi.

If you check these things and it still doesn't work, please upload your worksheet so that I can see what you actually typed.

 

 

Here is a procedure to generate the Matrix A such that A.variables = polynomials. I had to take off the (t)s to do this. That is, I changed a[1](t) to a[1], etc. They can be added back easily at any later time. If you prefer, I can change the procedure so that it takes them off and puts them back.




restart:

eq1:= a[1]*a[2]^2+a[2]*a[3]*a[1]+a[3]^3;

a[1]*a[2]^2+a[1]*a[2]*a[3]+a[3]^3

eq2:= a[1]^3*a[2]^2+a[2]*a[3]+a[2]^3;

a[1]^3*a[2]^2+a[2]^3+a[2]*a[3]

eq3:= a[1]*a[3]^2+a[3]^2*a[1]^3+a[2];

a[1]^3*a[3]^2+a[1]*a[3]^2+a[2]

(**)

MVfactor:= (eqns::{list,Vector}(polynom), vars::{list,Vector}(name))->
# Returns a Matrix A such that A.vars = eqns
     Matrix([seq([seq(quo(r,v,v,'r'), v= vars)], r= eqns)])
;                 

proc (eqns::(({Vector, list})(polynom)), vars::(({Vector, list})(name))) options operator, arrow; Matrix([seq([seq(quo(r, v, v, 'r'), v = vars)], r = eqns)]) end proc

(**)

eqns:= < eq||(1..3) >:
vars:= < seq(a[k], k= 1..3) >:
A:= MVfactor(eqns, vars);

A := Matrix(3, 3, {(1, 1) = a[2]^2+a[2]*a[3], (1, 2) = 0, (1, 3) = a[3]^2, (2, 1) = a[2]^2*a[1]^2, (2, 2) = a[2]^2+a[3], (2, 3) = 0, (3, 1) = a[1]^2*a[3]^2+a[3]^2, (3, 2) = 1, (3, 3) = 0})

 

Test (should get zero vector)

(**)

simplify(A.vars - eqns);

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

 


Download problem_4A.mw


Download problem_4A.mws

 

Change g:-`.` to g:-`*`.

Also, highly recommended: Change array to Array and get rid of evalm.

Recommended: Check out ifactors to replace ifactor. The output format is easier to work with in programs (do op manipulations, etc.); ifactor is mainly for interactive display.

No integer whatsoever has the stated property since 2014 is divisible by primes greater than 9, namely 19 and 53.

ifactor(2014);
                         (2) (19) (53)

A more interesting year for this question will be 2016 = 8*7*4*3*3.

The technique for solving this problem is essentially the same as for your previous problem.


(**)

restart:

readdata("C:/Users/Carl/desktop/t.txt"):

T := Vector(%, datatype = float[8]):

readdata("C:/Users/Carl/desktop/i2.txt"):

idata := Vector(%, datatype = float[8]):

(**)

Model:= i[s]*(1-exp(-(t/tau)^d));

i[s]*(1-exp(-(t/tau)^d))

FIT := Statistics:-NonlinearFit(Model, T, idata, [t], output = [residualsumofsquares, parametervalues]);

[0.212347302318356e-8, [d = HFloat(0.3048276105378924), tau = HFloat(1.361653413160606), i[s] = HFloat(1.1766215946973502e-4)]]

(**)

Curve:= plot(eval(Model, FIT[2]), t= min(T)..max(T)):

(**)

Points:= plot(< T | idata >, style= point, color= green):

(**)

plots:-display([Curve, Points]);

(**)

Model:= unapply(eval(Model, FIT[2]), t);

proc (t) options operator, arrow; HFloat(1.1766215946973502e-4)-HFloat(1.1766215946973502e-4)*exp(-HFloat(0.9101915620796802)*t^HFloat(0.3048276105378924)) end proc

(**)

writedata("C:/Users/Carl/desktop/out_i2.txt", convert(map(Model, T), list));

(**)

 

``


Download I-T.mw

It can be done as

N:= Sum(Sum(Sum('nn'(i, j, sigma), sigma= 1..2), j= 1 ..3), i= 1 ..f);

 

Please put questions in the Questions section rather than the Posts section. And please try to avoid screen shots -- they are not very readable, and they are not downloadable. You can upload worksheets instead.

Sorry about the 259/100 instead of 259/10---that was just a dumb mistake on my part.

Answering your question: It can be done with the writedata command:

NewModel:= unapply(eval(NewModel, FIT[2]), v):
writedata("C:/Users/Carl/desktop/out_i.txt", convert(map(NewModel, V), list));

It is not hard. It's a standard calculus II textbook exercise. Maple can do it directly, but here I show some steps of l'Hopital's rule.

 

restart:

(**)

y:= simplify(x/(x-1) - 1/ln(x));

(x*ln(x)-x+1)/((x-1)*ln(x))

(**)

eval(numer(y), x= 1), eval(denom(y), x= 1);

0, 0

Apply l'Hopital's rule:

(**)

y1:= simplify(diff(numer(y), x) / diff(denom(y), x));

x*ln(x)/(x*ln(x)+x-1)

(**)

eval(numer(y1), x= 1), eval(denom(y1), x= 1);

0, 0

Apply it again:

(**)

y2:= simplify(diff(numer(y1), x) / diff(denom(y1), x));

(ln(x)+1)/(ln(x)+2)

(**)

eval(numer(y2), x= 1)/eval(denom(y2), x= 1);

1/2

(**)

 

 

Download lHopitals.mw

The pattern repeats modulo 12 since:

ilcm(3,4);

                                        12

There are 167 periods:

iquo(2014, 12);
                                     167

with 2 occurrences each.

The partial period at the end is

irem(2014, 12);
                                      10

so that contributes another 2 occurrences.

2*167 + 2;
                                    336

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;

First 338 339 340 341 342 343 344 Last Page 340 of 382