Carl Love

Carl Love

28055 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

Use command march (short for Maple archive). Use march(create, ...) to create an archive. Use march(add, ...to add the .m file to the archive. Use march(list, ...to list the contents of the archive.

You use a global variable t symbolically; it's never given a value. Thus, the value returned from f is an expression depending on t[1], t[2], and t[3]. The value returned by f needs to be numeric if it's to be used by NLPSolve. If you do trace(f) before calling NLPSolve, you can see this.

Don't set warnlevel to 0 before your code is debugged. That's just asking for trouble. If the warning tells you that the variables are local, then declare them local.

You are using some very old coding styles. Perhaps you are emulating a textbook. It's giving you bad coding habits. There's no need for with(linalg). Change array to Array.

Don't rely on the with command to allow the use of so-called "short form" names inside a procedure. Instead, use a uses clause in the procedure header, or just use the "long form" name.

The following worksheet solves your original problem using Statistics:-NonlinearFit. Then it solves the same problem using NonlinearFit with procedural input. You can rewrite the procedure for your more-complicated model.

An example of the two ways to use Statistics:-NonlinearFit: with algebraic input (for simple models) and with procedural input (for more-complicated models)

 

Simple problem from Bronstejn-Semengyajev math book.

 

 

restart:

X:= ExcelTools:-Import("C:/users/owner/desktop/BroSzem_Data.xlsx", "XY", "A2:A13"):

X^+;

Matrix([[.1, .2, .3, .4, .5, .6, .7, .8, .9, 1.0, 1.1, 1.2]])

Y:= ExcelTools:-Import("C:/users/owner/desktop/BroSzem_Data.xlsx", "XY", "B2:B13"):

Y^+;

Matrix([[1.78, 3.18, 3.19, 2.54, 1.77, 1.14, .69, .4, .23, .13, 0.7e-1, 0.4e-1]])

XY:= <X|Y>:

fig1:= plot(XY, style= point, symbolsize= 16, view= [0..1.5, 0..3.5], gridlines= false);

Solution method 1: Algebraic input

Model:= a*(x^b)*exp(c*x);

a*x^b*exp(c*x)

MinSol:= Statistics:-NonlinearFit(
     Model, XY, x, initialvalues= [a=375, b=3, c=-7.5],
     output= [residualsumofsquares, parametervalues]
);

[0.916398873982076e-4, [a = HFloat(396.6019850867548), b = HFloat(1.9980985400641966), c = HFloat(-8.05357304777323)]]

fig2:= plot(eval(Model, MinSol[2]), x= 0..1.5):

plots:-display([fig||(1..2)], gridlines= false);

Solution method 2: Procedural input

Model:= proc(x, a, b, c)
     a*(x^b)*exp(c*x)
end proc:
   

MinSol:= Statistics:-NonlinearFit(
     Model, XY, initialvalues= [375, 3, -7.5],
     output= [residualsumofsquares, parametervector]
);

MinSol := [0.916398873982076e-4, Vector(3, {(1) = 396.6019850867548, (2) = 1.9980985400641966, (3) = -8.05357304777323}, datatype = float[8])]

params:= convert(MinSol[2], list)[];

HFloat(396.6019850867548), HFloat(1.9980985400641966), HFloat(-8.05357304777323)

ModelSubs:= subs(_P= params, x-> Model(x, _P));
     

proc (x) options operator, arrow; Model(x, HFloat(396.6019850867548), HFloat(1.9980985400641966), HFloat(-8.05357304777323)) end proc

fig3:= plot(ModelSubs, 0..1.5):

plots:-display([fig1,fig3], gridlines= false);

``

NULL

``

``

``

``

``

``


Download NonlinearFit_proc_input.mw

The default symbol for the imaginary unit in Maple is uppercase I, not the lowercase i that is more usual in mathematics and which you are using above. That default can be changed to any other symbol:

interface(imaginaryunit= i);

But beware that whatever symbol you use, you won't be able to use it as a variable. The use of i as a variable is so common in most programming languages that Maple decided to use the less commomly used symbol for the imaginary unit.

Form-conversion commands often need to be composed. The evalc command will take most expressions and express them as the sum of their real and imaginary parts. That is, the expression is put into cartesian or rectangular form. When polar is passed a single expression in this form, it converts it to polar form. Then evalf converts to decimal form.

polar(evalc(z));

     polar(sqrt(2), Pi/4)

evalf[4](%);

     polar(1.414, 0.7854)

The command randpoly might be a place to start:

randpoly(x);
randpoly(x, degree= 2);

etc.

Use literally the word undefined. Maple knows what this means.

A simpler way to solve the problem is to pass solve an inequality that forces the solution to be positive. This avoids the need for select, is, or another command line.

f:= gamma*sqrt(4) * sqrt(omega^2*C2^2*R4^2 /
     (C2^4*R4^4*beta^2*gamma^2*omega^4+C2^2*(1+gamma^2*(beta+1)^2-2*gamma) *
      R4^2*omega^2+1)):

solve({diff(f,omega), omega > 0}, omega) assuming positive;

     {omega = 1/(C2*R4*sqrt(beta*gamma))}

I have a hunch that the above is more robust than the flaky and ancient is because it uses the more modern and more sophisticated package SolveTools:-Inequality. However, I also know that it's still easy to confuse solve with inequalities, especially when square roots are involved.

An equivalent command is

solve(diff(f,omega), omega, useassumptions) assuming positive;

It's important for you to understand that the roots that you originally computed with RootFinding:-Isolate (or, equivalently, fsolve) are accurate to 10 digits. The large errors are introduced by your attempt to compute the residuals. Kitonum's Answer may give you the false impression that you need to increase Digits to compute the roots, although that wasn't his intention. No, you only need to increase Digits for the computation of the roots if you also want to compute the residuals, which isn't a necessary part of the computation of the roots.

restart:

p:=
     -12116320194738194778134937600000000*t^26 +
     167589596741213731838990745600000000*t^24 +
     1058345691529498270472972795904000000*t^22 -
     4276605572538658673086219419648000000*t^20 -
     23240154739806540070988490473472000000*t^18 -
     5442849111209103187871341215744000000*t^16 +
     49009931453396028716875310432256000000*t^14 +
     74247033158233643322704589225984000000*t^12 -
     2762178990802317464801412907008000000*t^10 -
     25947900993773120244883450232832000000*t^8 -
     7468990043547273070742668836864000000*t^6 -
     567730116675454293925108383744000000*t^4 +
     3703566799705707258760396800000000*t^2 -
     4742330812072533924249600000000
:

Some exact computations can be used to simplify your polynomial. These aren't necessary.

p1:= p/icontent(p):   #Remove common integer factor.

p2:= algsubs(t^2= T, p1):   #The polynomial is even.

p3:= unapply(convert(p2, horner), T)@(t-> t^2):   #Reduce amount of computation for residuals.

The rest of the computations are floating point.

Digits:= 10:

display:= proc(e, digits::posint) print(evalf[digits](e)); e end proc:

Compute the positive roots.

R10:= display([fsolve(p2, T= 0..infinity)^~(1/2)], 10):

[0.4212510777e-1, 0.6596008501e-1, .8048430445, 1.300314688, 2.295186769, 4.162501845]

Compute the residuals.

display(p3~(R10), 2):

[0.2e8, 0.1e8, -0.67e13, -0.32e17, 0.14e22, 0.11e27]

Now we do the exact same computations at 50 digits.

Digits:= 50:

R50:= display([fsolve(p2, T= 0..infinity)^~(1/2)], 10):

[0.4212510777e-1, 0.6596008501e-1, .8048430445, 1.300314688, 2.295186769, 4.162501845]

Note that the above roots are exactly the same (to 10 digits) as the previously computed roots.

display(p3~(R50), 2):

[0., 0.3e-32, -0.12e-26, 0.32e-23, -0.92e-19, 0.42e-13]

...but the residuals computed at 50 digits are much smaller in magnitude. It's also interesting that as we increased Digits by 40, the magnitudes were reduced by a factor of about "10^(40)." That's typical, but not guaranteed.



Download accuracy_of_roots.mw

To say that a matrix is "augmented", in this context, means that the right-side vector b has been appended as the rightmost column of the matrix. To achieve this, you simply need to make the command

IterativeApproximate(<A|b>, ...);

So, the b is incorporated into the first argument. In this case, don't use b as the second argument.

I get this solution:

     Vector[column](4, [3.98912386882810, 2.99292198093478, 2.99317583909596, 2.99337057407024])

which has residual

LinearAlgebra:-Norm(A.% - b);

     0.103910640752635786e-2

which is very close to your specified tolerance of 1e-3.

The above is a response to your originally posted Question, using b = <1,1,1,1>.

Responding to your reformulated Question, using b = <1,0,0,0>: The solution that you got seems correct to me. The residual is 0.235228835640644007e-3, which is well less than your specified tolerance.

 

It's important for a beginner to note that print doesn't "give" anything that can be stored for later use. It's often difficult for a beginner to understand this. The print command only displays something on the screen; it's rarely used in formal Maple programming.

Here's how to create A in one line of code and then create in a second line of code:

Download Matrices.mw
A:= Matrix(

     [[0$7], [0$6], [a,-d,a,-d,a], [-a,d,d,a], -[d$3], [d,d]],
     scan= band[0,5],
     shape= antisymmetric
);

C:= [seq(seq(eval(A, {a= i, d= j}), i= [-1,1]), j= [-1,1])];

Now C is a list of four matrices. The matrices can be accessed as C[1]C[2], etc. Note that for loops can often be replaced by seq expressions, and it's usually beneficial to do so.

Example:

A:= Array([4,3,7]);

If your objective function is a procedure, then the constraints must also be procedures. Your current constraints are expressions. So your command should be like this:

sol:= Optimization:-NLPSolve(
     f, {}, {(x1,x2,x3,x4)-> x1+x2+x3+x4-65}, 0..20, 0..15, 0..20, 0..15, minimize
);

The Units subpackages (Standard in this case) overload the arithmetic operators to behave like this. You can temporarily revert any operator to its default meaning by giving it the :- prefix and using it in prefix form. So, for this case, you could do

eq:= :-`+`(2.*Unit(m), d) = 5.*Unit(m);

     (the equation will display normally: the + will be infix.)

Now you can pass this equation to solve:

solve(eq, d);

e:= w^2/(w^4+2*w^2+1):
expand(numer(e)/w^2)/expand(denom(e)/w^2);

The first expand isn't really needed in this case.

First 218 219 220 221 222 223 224 Last Page 220 of 395