Applications, Examples and Libraries

Share your work here

Hi
It's been years since expressions like A %* B %+ C involving inert arithmetic operators used in infix form are correctly understood (parsed) when written on a 1D-Math input line. The idea is simple: have the operators %. %*, %+, %-, %^, %/ work on input as infix operators the same way their active forms: ., *, +, -, ^, / do. This useful functionality, however, remained elusive when using 2D-Math input notation, so one would have to resort to using the functional form of the operators. E.g., input the above as `%+`(`%*`(A, B) ,C), which for me is really ugly. Besides being a bit demoralizing: we do all this fuzz about how great computer algebra and 2D-Math input notation is, and then input things in that way …

So this is to mention that this elusive functionality of inert arithmetic operators used in infix form within a 2D-Math input line now works. The novelty is present in the latest Maplesoft Physics Updates for Maple 2023, which is version 1490. As usual, to install the Updates open a Maple worksheet and input Physics:-Version(latest).

Here is an image (worksheet at the end) showing the new thing


The implementation is pretty new; reports of anything related to these inert operators not displayed/working as you'd expect are much appreciated. 


Download Inert_arithmetic_operators_in_2D_Math.mw

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

Hi everyone,

I'd like to draw your attention to a package we recently uploaded to the Maple Cloud, here. You can download the package from the linked Cloud page, or directly from here as a workbook file: NaturalLanguage.maple. I'll include a lightly edited version of the "cover sheet" that introduces the package below. I have left the first four sections folded closed - you can see those in the linked Maple Cloud page or the workbook - because I think the last section is the most interesting.
 

Using natural language models in Maple


This package explores using large language models such as ChatGPT for processing natural language in Maple. Let's load the package from source, save it in this workbook, and load it.

restart

read "this:///lib/NaturalLanguage.mpl"

savelib('NaturalLanguage', "this://")

with(NaturalLanguage)

Warning, the NaturalLanguage package is an experimental package designed as an interface to publicly available large language models such as ChatGPT. Use of this package is entirely at the user's risk. Results may be inconsistent, misleading, or flat out incorrect. 

[Explain, GetCommand, GetMath, Query, RawQuery]

(1)

We note there is a warning: we will see output of large language models, which will often include inaccurate statements. Please keep this in mind - this is not (yet?) technology that you want to use to build a bridge!

We also note that there are five publically exposed commands. Let's look at them in turn.

RawQuery

   

Query

   

GetMath

   

GetCommand

   

Explain

 
• 

The Explain command asks the given model to explain the item in the query. After an explanation, you can ask for more detail by issuing the command Explain(more), Explain("go on"), Explain(elaborate), Explain(further), Explain("continue"), or minor variations. The item you ask about can be a mathematical expression, equation, list of equations, or a string.

• 

The default model is GPT-4; you can select the other model by using the model = ChatGPT option, or forcing use of the GPT-4 model by using the model = GPT4 option. Note that OpenAI may deprecate and disable models, so the set of models supported may change in the future.

display(Explain(x^2+y^2=1));display(Explain(x^2+y^2 = 1))

This equation represents a circle with a radius of 1 centered at the origin (0,0) on a coordinate plane.

 

display(Explain(more));display(Explain(more))

The equation x^2 + y^2 = 1 represents a circle in a two-dimensional plane (specifically, a Cartesian plane).

The center of the circle is at the origin of the coordinates (0,0) and its radius is 1. This is because any point (x, y) on the circle is a distance r (the radius) away from the center, and by the Pythagorean theorem, this distance is given by the square root of (x^2 + y^2). The fact that x^2 + y^2 = 1 implies that the radius r = sqrt(1) = 1.

This fundamental equation is also the unit circle in trigonometry, where angles are measured in radians. The coordinates (x, y) represent the cosine and sine of the angle respectively. The unit circle is a crucial concept in trigonometry, complex number theory, and calculus. It simplifies many mathematical concepts and provides a geometric interpretation of a variety of phenomena in physical sciences and engineering.

 

display(Explain(Re(exp(x*I + y*I)) = Re(exp(x*I)) * Re(exp(y*I)) - Im(exp(x*I)) * Im(exp(y*I))));display(Explain(Re(exp(I*x+I*y)) = Re(exp(I*x))*Re(exp(I*y))-Im(exp(I*x))*Im(exp(I*y))))

The equation seems to be using Euler's formula, which states that e^(ix) = cos(x) + i*sin(x). The 'Re' refers to the real part of a complex number, and 'Im' refers to the imaginary part. So, if you put Euler's formula in there your equation is equivalent to cos(x + y) = cos(x)cos(y) - sin(x)sin(y), which is basically the well-known trigonometric formula for the cosine of a sum. In a nutshell, it's dealing with the trigonometric identities and properties of exponential functions with complex numbers.

 

display(Explain("go on", model = ChatGPT));display(Explain("go on", model = ChatGPT))

The given equation appears to be dealing with the real part of a complex exponential expression. Let's break it down step by step to understand its meaning:

1. The equation involves the expression e^(i*x+i*y), where e represents Euler's number (~2.718), i is the imaginary unit (√(-1)), and x and y are real numbers.
2. The expression e^(i*x+i*y) can be rewritten using Euler's formula: e^(i*x+i*y) = cos(x+y) + i*sin(x+y). This is a complex number in the form a + bi, where 'a' represents the real part (Re) and 'b' represents the imaginary part (Im).
3. The real part of e^(i*x+i*y) is, therefore, Re(e^(i*x+i*y)) = cos(x+y).
4. The right side of the equation consists of the product of two complex exponential expressions: Re(e^(i*x))*Re(e^(i*y))-Im(e^(i*x))*Im(e^(i*y)).
5. The expression Re(e^(i*x)) represents the real part of e^(i*x), which is cos(x), and Im(e^(i*x)) represents the imaginary part of e^(i*x), which is sin(x).
6. Similarly, Re(e^(i*y)) represents the real part of e^(i*y), which is cos(y), and Im(e^(i*y)) represents the imaginary part of e^(i*y), which is sin(y).
7. Therefore, the right side becomes: (cos(x)*cos(y)) - (sin(x)*sin(y)).

So, the equation Re(e^(i*x+i*y)) = Re(e^(i*x))*Re(e^(i*y))-Im(e^(i*x))*Im(e^(i*y)) essentially states that the real part of e^(i*x+i*y) (i.e., cos(x+y)) is equal to the real part of e^(i*x) (i.e., cos(x)) multiplied by the real part of e^(i*y) (i.e., cos(y)), minus the imaginary part of e^(i*x) (i.e., sin(x)) multiplied by the imaginary part of e^(i*y) (i.e., sin(y)).

In simpler terms, the equation is connecting the real part of the sum of two complex exponential expressions with the product of the real parts and the difference of the imaginary parts of those individual expressions.

 

display(Explain(Re(sigma) < 0));display(Explain(Re(sigma) < 0))

The inequality "Re(sigma) < 0" is a statement in complex number theory. It is saying that the real part (Re) of a complex number sigma is less than zero. This means that the complex number is to the left of the imaginary axis in the complex number plane. This kind of inequality is often used in the study of complex functions, in areas such as control systems or signal processing in electrical engineering, and also in number theory and the study of dynamical systems in mathematics.

 

display(Explain(continue));display(Explain(continue))

The inequality "Re(sigma) < 0" is typically related to the field of complex numbers. Here, "Re" stands for the real part of a complex number and "sigma" is usually used to represent a complex number. So, "Re(sigma) < 0" is saying that the real part of the complex number sigma is negative.

Complex numbers are usually written in the form a + bi, where a and b are real numbers and i is the imaginary unit. The "Re" in "Re(sigma)" refers to the a part of this form, i.e. the real part of the complex number.

This inequality is often found in fields like control systems engineering or signal processing. For instance, in control systems, all the poles of a stable system must have a negative real part. Another context you might come across this inequality is in complex analysis, one might examine the region where the real part of a complex function is negative. So, the inequality "Re(sigma) < 0" expresses that a complex number is in the left half of the complex plane, which is often crucial for the stability of many systems in engineering.

 

display(Explain("The Maple command simplify"));display(Explain("The Maple command simplify"))

The Maple command "simplify" is used to simplify mathematical expressions. In Maple, this command applies various rules to reduce an equation or expression to its simplest or most compact form.

For example, if you have an equation like `(x^3 - x^2)/x^2`, using the simplify command would reduce it to `x - 1`. It simplifies algebraic expressions, performs factorization, combines like terms, cancel out common factors, etc.

The general syntax is `simplify(expr)`, where `expr` is the expression that you want to simplify. The simplified expression is returned as output.

Please note that the simplifications made by this command are based on general rules and may not always give the expected results for some specific expressions due to the complexity of mathematical simplifications.

 

display(Explain(continue, model=ChatGPT));display(Explain(continue, model = ChatGPT))

The Maple command "simplify" is used to simplify mathematical expressions. It takes an expression as input and attempts to find a simplified form by applying various algebraic rules and simplification techniques.

When you use the "simplify" command, Maple applies a series of simplification steps, which can include:

1. Expanding and factoring polynomials: This step expands products of terms and factors out common terms.
2. Simplifying fractions: The command simplifies fractions by canceling common factors in the numerator and denominator.
3. Simplifying trigonometric expressions: Maple applies trigonometric identities to simplify expressions involving trigonometric functions.
4. Simplifying logarithmic and exponential expressions: The command applies rules of logarithms and exponents to simplify expressions with these functions.
5. Simplifying square roots and radicals: It simplifies square roots by expressing them in simplest radical form.
6. Simplifying expressions involving absolute values: Maple simplifies absolute value expressions by considering different cases based on the sign of the argument.
7. Applying mathematical identities and rules: The command uses various mathematical identities and rules to simplify expressions further.

The "simplify" command in Maple is very powerful and can handle a wide range of mathematical expressions. However, it may not always produce the most simplified form, as simplification can sometimes be a subjective task. In such cases, you can try using other commands or options available in Maple to refine the simplification process.

To use the "simplify" command in Maple, you can simply type "simplify(expression)" where "expression" represents the mathematical expression you want to simplify.

 

display(Explain("Weight of a W10X54 beam?"))

The weight of a W10X54 beam refers to the weight specification of a beam in the W Beam series, commonly used in construction and engineering projects. The nomenclature "W10X54" is decoded as follows:

- W: it stands for 'wide flange', a type of structural steel beam where the flanges are parallel to each other.
- 10: it indicates the nominal depth of the wide flange, which is about 10 inches.
- 54: it indicates the weight of the beam per linear foot, which is about 54 pounds.

So, a W10X54 beam weighs 54 pounds per foot.

 

display(Explain(further))

"Weight of a W10X54 beam" refers to the weight of a specific kind of structural steel beam used in construction and architecture. The "W" stands for "wide flange", which is a particular type of I-beam.

The "10" means the nominal depth of the beam in inches (also known as its height when standing upright).

The "54" represents the beam's weight per foot in pounds. In other words, for every foot of length the beam has, it weighs 54 pounds.

 So if you have a W10X54 beam that is 10 feet long, its total weight would be 540 pounds (10 feet times 54 pounds per foot).

Please note that these are nominal dimensions and may not exactly match the true dimensions of the beam due to manufacturing variables.

 

NULL

 

What do you think? Will this be an important part of mathematical software in the near future, or are we still far away from that point? We'd love to hear your opinions.

Maple's fsolve command can quickly solve expressions involving large floating point numbers where the (symbolic) solve command can take much longer attempting to solve the equivalent rational expression. For example, consider the following worksheet:

restart;

sys := {a + b^0.2784982189 = c+1, b + c^0.9575068354 = a+2, c + a^0.1576130817 = b+3};

{a+b^.2784982189 = c+1, b+c^.9575068354 = a+2, c+a^.1576130817 = b+3}

(1)

fsolve_start:=time[real]():

fsolve(sys);

{a = 3.561242843, b = 1.994950678, c = 3.773320855}

(2)

fsolve_elapsed_seconds:= time[real]()-fsolve_start;

0.50e-1

(3)

solve_start:=time[real]():

###warning, the following command may crash and/or execute indefinitely###

solve(sys);

solve_elapsed_hours:=(time[real]()-solve_start)/3600;


 

Download solve-fsolve-primes.mw

If you've seen Paulina's announcement then you know that we are once again holding a virtual Maple Conference this year.  As well, we are once again going to have a virtual gallery featuring artwork and creative projects submitted by the Maple community!

Last year we had a number of great submissions to our Maple Art Gallery and our Maple Learn Creative Showcase.  These were our excellent prize winners.

From left to right we have A visualization of all the primitive roots of 10037 created by Simon Plouffe, winner of the Judge’s Choice, Mother’s Day Rose created with Maple plots by Greg Wheaton, winner of the People’s Choice, and Mona Lisa in Maple Learn created by Paul DeMarco (with help from Leonardo DaVinci), the winner of the People’s Choice for the Maple Learn Showcase.

This year we are expanding the Gallery into two collections to encourage more people to submit.  They are

  • The Art Gallery - A small gallery to highlight high effort, mathematically interesting works (with stricter criteria)

  • The Creative Works Showcase - A larger showcase for nearly any interesting visual works created with Maplesoft products like Maple Learn and Maple

Feel free to submit nearly anything cool for the Creative Works Showcase, if we find it particularly impressive we might even ask you to let us consider it for the gallery.  Also, do not be intimidated by the title "Art Gallery" we're looking for anything that has taken some artistic effort and tells a mathematical story.

For more information on critera and how to submit, please visit our Call for Creative Works.  The important deadline to know is the September 14th deadline for submission of works with virtual gallery reception and awards ceremony durring the conference October 26-27.

I look forward to seeing all the submissions for the Maple community again this year!

We've just launched Maple Flow 2023!

The new release offers many enhancements that help you calculate and write reports faster, resulting in polished technical documents. Let me describe a few of my favorite new features below.

You can now change the units of results inline in the canvas, without taking your hands off the keyboard. You can still use the Context Panel, but the new method is faster and enhances the fluid workflow that Flow exemplifies.

You can also enter a partial unit inline; Flow will automatically insert more units to dimensionally balance the system.

This is useful when results are returned in base dimensions (like time, length and mass) but you want to rescale to higher-level derived units. For an energy analysis, for example, you might guess that the result should contain units of Joules, plus some other units, but you don't know what those other units are; now, you can request that the result contains Joules, and Flow fills the rest in automatically.

The new Variables Palette lists all the user-defined variables and functions known to Flow at the point of the cursor. If you move your grid cursor up or down, the variables palette intelligently removes or adds entries.

You can now import an image by simply dragging it from a file explorer into the canvas.

This is one of those small quality-of-life enhancement that makes Flow a pleasure to use.

You can now quickly align containers to create ordered, uncluttered groups.

We've packed a lot more into the new release - head on over here for a complete rundown. And if you're tempted, you can get a trial here.

We have a lot more in the pipeline - the next 12 months will be very exciting. Let me know what you think!

On this day 181 years ago, Christian Doppler first presented the effect that would later become known as the Doppler effect. In his paper “On the coloured light of the binary stars and some other stars of the heavens”, he proposed (with a great deal of confidence and remarkably little evidence) that the observed frequency of a wave changes if either the source or observer is moving. Luckily for Doppler, he did turn out to be right! Or at least, right about the effect, not right about supernovas actually being binary stars that are moving really fast. The effect was experimentally confirmed a few years later, and it’s now used in a whole variety of interesting applications.

To learn more about how the Doppler effect works, take a look at this Maple MathApp. You can adjust the speed of the jet to see how the frequency of the sound changes, and add an observer to see what they perceive the sound as. You can even break the sound barrier, although the poor observer might not like that so much!

 

A screenshot of a Maple MathApp, showing a visual representation of sound waves coming off a moving jet, with sliders to adjust the speed.

 

For Maple users, you can also check out the MathApp on the relativistic Doppler effect. You’ll find it in the Natural Sciences section, under Astronomy and Earth Sciences. Settle in to watch those colours come to life!

A screenshot of a Maple MathApp showing a spectrum of colours, with sliders to control the initial wavelength of the light and a dial showing the current velocity of the viewer

 

But wait, I mentioned interesting applications, didn’t I? And don’t worry, I’m not just here to talk about sirens moving past you or figuring out the speed of stars (although admittedly, that one is pretty interesting too). No, I’m talking about robots. Some robots make use of the Doppler effect to help monitor their own speed, by bouncing sound waves off the floor and measuring the frequency of the reflected wave. A large change in frequency means that robot is zooming!

The Doppler effect is also used in the medical field—Doppler ultrasonography uses the Doppler effect to determine and visualize the movement of tissues and body fluids like blood. It works by bouncing sound waves off of moving objects (like red blood cells) and measuring the result. The difference in frequency tells you the speed and direction of the blood flow, in accordance with the Doppler effect! Pretty neat, if you ask me.

And like any good scientific phenomena, the Doppler effect can be used for both work and pleasure. The Leslie speaker is a type of speaker invented in the 1940s that modifies the sound by rotating a baffle chamber, or drum, in front of the loudspeakers. The change in frequency dictated by the Doppler effect causes the pitch to fluctuate, creating a distinct sound that I can only describe as “woobly”. The speaker can be set to either “chorus” or “tremolo”, depending on how much woobliness the user wants. It was typically used with the Hammond organ, and you can hear it in action here!

You know who else uses the Doppler effect? Bats. Since they rely on echolocation to get around, they need some way to account for the fact that the returning sound waves won’t be at the same pitch that they were sent out at. This fantastic video explains it far better than I ever could, and involves putting bats on a swing, which I think should be enough of a recommendation all on its own.

That’s it for our little foray into the Doppler effect, although there’s still a lot more that could be said about it. Try checking out those Maple MathApps for inspiration—who knows, maybe you’ll find a whole new use for this fascinating effect!

Ever wonder how to show progress updates from your executing code without printing new lines each time?

One way to do this is to use a TextArea component and the DocumentTools package. The TextArea could be inserted from the Components Palette in Maple, or programmatically like so:

restart;

with(DocumentTools):

with(DocumentTools:-Components):

with(DocumentTools:-Layout):

s := "0": #initial text value

T := TextArea(s, identity = "TextArea0"):
xml := Worksheet(Group(Input(Textfield(T)))):

insertedname:=InsertContent(xml)[1,1]: #find the inserted component name in case changed

for i to 10 do #start the demonstration procedure
   Threads:-Sleep(1);
   SetProperty(insertedname,value,sprintf("%d",i),refresh=true);
end do:

Maplets:-Examples:-Message("Done");


Download text-area-update-progress.mw

This post discusses a solution for modeling a traveling load on Maplesim's Flexible Beam component and provides an example of a bouncing load.

The idea for the above example came from an attempt to reproduce a model of a mass sliding on a beam from MapleSim's model gallery. However, reproducing it using contact components in combination with the Flexible Beam component turned out to be not straightforward, and this will be discussed in the following.

To simulate a traveling load on the Flexible Beam component, one could apply forces at discrete locations for a certain duration. However, the fidelity of this approach is limited by the number of discrete locations, which must be defined using the Flexible Beam Frame component, as well as the way in which the forces are activated.

One potential solution to address the issue of temporal activation of forces is to attach contact elements (such as Rectangle components) at distinct locations along the beam, which are defined by Flexible Beam Frame components, and make contact using a spherical or toroidal contact element. However, this approach also introduces two new problems:

  • An additional bending moment is generated when the load is not applied at the center of the contact element's attachment point, the Flexible Beam Frame component. Depending on the length of the contact element, deformations caused by this moment can be greater than the deformation caused by the force itself when the force is applied at the ends of the contact element. Overall, this unwanted moment makes the simulation unrealistic and must be avoided.
  • When the beam bends, a gap (see below) or an overlap is created between adjacent Rectangle components. If there is a gap, the object exerting a force on the beam can fall through it. Overlaps can create differences in dynamic behavior when the radius of curvature of the beam is on the opposite side of the point of contact.

To avoid these problems, the solution presented here uses an intermediate kinematic chain (encircled in yellow below) that redistributes the contact force on the Rectangle component on two support points (ports to attach Flexible Beam Frame components) in a linear fashion.

 

 

To address gaps, the contact element (Rectangle) attached to the kinematic chain has the same width as the chain and connects to the adjacent contact elements via multibody frames. In the image below, 10 contact elements are laid on top of a single Flexible Beam component, like a belt made out of tiles. The belt has to be pinned to the flexible beam at one location (highlighted in yellow). The location of this fixed point determines how the flexible beam is loaded by tangential contact forces (friction forces) and should be selected carefully.

 

 

Some observations on the attached model:

  • Low damping and high initial potential energy of the ball can result in a failed simulation (due to constraint projection failure). Increasing the number of elastic coordinates has a similar effect. Constraint projection can be turned off in the simulation settings to continue simulation.
  • The bouncing ball excites several eigenmodes at once, causing the beam to wiggle chaotically in combination with the varying bouncing frequency of the ball. A similar looking effect can also be achieved with special initial conditions, as demonstrated with Maple in this excellent post on Euler beams and partial differential equations.
  • Repeated simulations with low damping lead to different results (an indication of chaotic behavior; see three successive simulations below (gold) compared to a saved solution(red)). The moment in the animation when the ball travels backward represents a metastable equilibrium point of the simulation. This makes predictions beyond this point difficult, as the behavior of the system is highly dependent on the model parameters. Whether the reversal is a simulation artifact or can happen in reality remains to be seen. Overall, this example could evolve into a nice experimental fun project for students.

  • Setting the gravitational constant for Mars, everything is different. I could not reproduce the fun factor on Earth. A reason more to stay ,-)

Ball_bouncing_on_a_flexible_beam.msim

 

 

2-dimensional motion and displacement are some of the first topics that high school students learn in their physics class. In my physics classes, I loved solving 2-dimensional displacement problems because they require the use of so many different math concepts: trigonometry, coordinate conversions, and vector operations are all necessary to solve these problems. Though displacement problems can seem complicated, they are easy to visualize.
For example, below is a visualization of the displacement of someone who walked 10m in the direction 30o North of East, then walked 15m in the direction 45o South of East:

From just looking at the diagram, most people could identify that the final position is some angle Southeast of the initial position and perhaps estimate the distance between these two positions. However, finding an exact solution requires various computations, which are all outlined in the Directional Displacement Example Problem document on Maple Learn.

Solving a problem like this is a great way to practice solving triangles, adding vectors, computing vector norms, and converting points to and from polar form. If you want to practice these math skills, try out Maple Learn’s Directional Displacement Quiz; this document randomly generates displacement questions for you to solve. Have fun practicing!

 

TODAY I GOT AN INSPIRATION TO CREATE 3D GRAPH EQUATION OF WALKING ROBOT (ED-209) IN CARTESIAN SPACE USING ONLY WITH SINGLE IMPLICIT EQUATION.

ENJOY...

 

How to Create Graph Equation of Wankel Engine on Cartesian Plane using Single Implicit Function run by Maple Software

Enjoy...

 

It seems once in a while someone asks about converting from degrees minutes seconds format to decimal or the other way around. 

I created a little procedure for just that purpose. 

restart; gc()

NULL

A little procedure to convert the form of degrees, minutes, seconds to decimal and vice versa.  

NULL

dms := proc (d, m := 0, s := 0) local con, d1, d2, d3; if 0 < frac(d) then d1 := floor(d); d2 := floor(60*frac(d)); d3 := 60*frac(60*frac(d)); con := cat(d1, `° `, d2, `' `, d3, `"`) else con := d+m/60.+s/3600. end if; print(con) end proc

NULL

Examples of use.

 

dms(45.2365)

`45° 14' 11.4000"`

(1)

dms(45, 14, 11.4)

45.23650000

(2)

NULL

Download dms.mw

edit - a quick realization is to remove the decimals that changes the fractions to floating decimals and change print(con) to print(evlaf(con)) to avoid rounding issues.

This post is motivated by a question asked by @vs140580  ( The program is making intercept zero even though There is a intercept in regression Fit (A toy code showing the error attached) ).

The problem met by @vs140580 comes from the large magnitudes of the (two) regressors and the failure to Fit/LinearFit to find the correct solution unless an undecent value of Digits is used.
This problem has been answerd by @dharr after scaling the data (which is always, when possible, a good practice) and by 
myself while using explicitely the method called "Normal Equations" (see https://en.wikipedia.org/wiki/Least_squares).

The method of "Normal Equations" relies upon the inversion of a symmetric square matrix H whose dimension is equal to the number of coefficients of the model to fit.
It's well known that this method can potentially lead to matrices H extremely ill-conditionned, and thus face severe numerical problems (the most common situation being the fit of a high degree polynomial).
 

About these alternative methods:

  • In English: http://www.math.kent.edu/~reichel/courses/intr.num.comp.1/fall09/lecture4/lecture4.pdf
  • In French: https://moodle.utc.fr/pluginfile.php/24407/mod_resource/content/5/MT09-ch3_ecran.pdfI


The attached file illustrates how the QR decomposition method works.
The test case is @vs140580's.

Maybe the development team could enhance Fit/LinearFit in future versions by adding an option which specifies what method is to be used?

 

restart:

with(Statistics):

interface(version)

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

(1)

Data := Matrix([[4.74593554708566, 11385427.62, 2735660038000], [4.58252830679671, 25469809.77, 12833885700000], [4.29311160501838, 1079325200, 11411813200000000], [4.24176959154225, 1428647556, 18918585950000000], [5.17263072694618, 1428647556, 18918585950000000], [4.39351114955735, 1877950416, 30746202150000000], [4.39599006758777, 1428647556, 18918585950000000], [5.79317412396815, 2448320309, 49065217290000000], [4.48293612651735, 2448320309, 49065217290000000], [4.19990181982522, 2448320309, 49065217290000000], [5.73518217699046, 1856333905, 30648714900000000], [4.67943831980476, 3071210420, 75995866910000000], [4.215240105336, 2390089264, 48670072110000000], [4.41566877563247, 3049877383, 75854074610000000], [4.77780395369828, 2910469403, 74061327950000000], [4.96617430604669, 1416936352, 18891734280000000], [4.36131111330988, 1416936352, 18891734280000000], [5.17783192063198, 1079325200, 11411813200000000], [4.998266287191, 1067513353, 11402362980000000], [4.23366152474871, 2389517120, 48661380410000000], [4.58252830679671, 758079709.3, 5636151969000000], [6.82390874094432, 1304393838, 14240754750000000], [4.24176959154225, 912963601.2, 8621914602000000], [4.52432881167557, 573965555.4, 3535351888000000], [4.84133601918601, 573965555.4, 3535351888000000], [6.88605664769316, 732571773.2, 5558875538000000], [4.35575841415627, 1203944381, 13430693320000000], [4.42527441640593, 955277678, 8795128298000000], [6.82390874094432, 997591754.9, 8968341995000000], [4.35144484433733, 143039477.1, 305355143300000]]):

# Direct use of LinearFit.
#
# As far as I know LinearFit is based on the resolution of the "Normal Equations"
# (see further down), a system of equations that is known to be ill-conditioned
# when regressors have large values (in particular when polynomial regression
# is used).

X := Data[.., [2, 3]]:
Y := Data[.., 1]:


LinearFit(C1+C2*v+C3*w, X, Y, [v, w]);

Warning, model is not of full rank

 

HFloat(6.830889923844425e-9)*v-HFloat(2.275143726335622e-16)*w

(2)

# For roundoff issues the 3-by-3 matrix involved in the "Normal Equations" (NE)
# appears to of rank < 3.
# The rank of this matrix is rqual to 1+rank(X) and one can easily verify that
# the 2 columns of X are linearly independent:

LinearAlgebra:-LinearSolve(X, Vector(numelems(Y), 0));
LinearAlgebra:-Rank(X);

 

Vector[column]([[0.], [-0.]])

 

2

(3)

# Solve the least squares problem by using explicitely the NE.
#
# To account for an intercept we augment X by a vector column of "1"
# traditionally put in column one.
Z := `<|>`(Vector(numelems(Y), 1), X):  
A := (Z^+ . Z)^(-1) . Z^+ . Y;          # Normal Equations

A := Vector(3, {(1) = 4.659353816079307, (2) = 0.5985084089529947e-9, (3) = -0.27350964718426345e-16}, datatype = float[8])

(4)

# What is the rank of Z?
# Due to the scale of compared to "1", Rank fails to return the good value
# of rank(Z), which is obviously equal to rank(X)+1.

LinearAlgebra:-LinearSolve(Z, Vector(numelems(Y), 0));
LinearAlgebra:-Rank(Z);

Vector[column]([[0.], [0.], [-0.]])

 

2

(5)


A WORKAROUND : SCALING THE DATA

model := unapply( LinearFit(C1+C2*v+C3*w, Scale(X), Scale(Y), [v, w]), [v, w] );

proc (v, w) options operator, arrow; -HFloat(1.264691577813453e-15)+HFloat(0.6607154853418553)*v-HFloat(0.8095150669884322)*w end proc

(6)

mX, sX := (Mean, StandardDeviation)(X);
mY, sY := (Mean, StandardDeviation)(Y);

mX, sX := Vector[row](2, {(1) = 1447634550.7963333, (2) = 24441399854567932.}, datatype = float[8]), Vector[row](2, {(1) = 871086770.7242773, (2) = 23354440973344224.}, datatype = float[8])

 

HFloat(4.857279402730572), HFloat(0.789073010656694)

(7)

MODEL := model((x1-mX[1])/sX[1], (x2-mX[2])/sX[2]) * sY + mY

HFloat(4.659353816079309)+HFloat(5.985084089530032e-10)*x1-HFloat(2.7350964718426736e-17)*x2

(8)

# Check that the vector of regression coefficients is almost equal to A found above
# relative error lesst than 10^(-14)

A_from_scaling       := < coeffs(MODEL) >:
Relative_Discrepancy := (A_from_scaling - A) /~ A

Relative_Discrepancy := Vector(3, {(1) = 0.5718679809000842e-15, (2) = 0.14166219140514066e-13, (3) = 0.14308415396983588e-13}, datatype = float[8])

(9)


THE QR DECOMPOSITION  (applied on raw data)

The QR decomposition, as well as Given's rotation method, are two alternatives to the the NE method
to find the vector of regression coefficients.
Both of them are known to be less sensitive to the magnitudes of the regressors and do nt require (not
always) a scaling of the data (which can be quite complex with polynomial regression or when some
transformation is used to liearize the statistical model, for instanc Y=a*exp(b*X) --> log(Y)=log(a)+b*X).

N := numelems(Y);
P := numelems(Z[1]);

30

 

3

(10)

# Perform the QR decomposition of Z.

Q, R := LinearAlgebra:-QRDecomposition(Z, fullspan);

Q, R := Vector(4, {(1) = ` 30 x 30 `*Matrix, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order}), Vector(4, {(1) = ` 30 x 3 `*Matrix, (2) = `Data Type: `*float[8], (3) = `Storage: `*triangular[upper], (4) = `Order: `*Fortran_order})

(11)

# Let C the column vector of length P defined by:

C := (Q^+ . Y)[1..P];

C := Vector(3, {(1) = -26.6044149698075, (2) = -.517558353158176, (3) = -.881007519371895})

(12)

# Then the vector of regression coefficients is given by:

A_QR                 := (R[1..P, 1..P])^(-1) . C;
Relative_Discrepancy := (A_QR - A) /~ A

A_QR := Vector(3, {(1) = 4.65935381607931, (2) = 0.5985084090e-9, (3) = -0.2735096472e-16})

 

Relative_Discrepancy := Vector(3, {(1) = 0.3812453206e-15, (2) = 0.1105656128e-13, (3) = 0.1216778632e-13})

(13)

# The matrix H = Z^+ . Z writes

H                    := Z^+ . Z:
H_QR                 := R^+ . Q^+ . Q . R:
Relative_Discrepancy := (H_QR - H) /~ H

Relative_Discrepancy := Matrix(3, 3, {(1, 1) = -0.1184237893e-15, (1, 2) = 0., (1, 3) = -0.3491343943e-15, (2, 1) = 0., (2, 2) = 0.1930383052e-15, (2, 3) = 0.3369103254e-15, (3, 1) = -0.1745671971e-15, (3, 2) = -0.5053654881e-15, (3, 3) = -0.1366873891e-15})

(14)

# H_QR expression is required to obtain the covariance matrix of the regression coefficients.


 

Download LeastSquares_and_QRdecomposition.mw


 

 

CREATING GRAPH EQUATION OF "DNA" IN CARTESIAN SPACE USING PARAMETRIC SURFACE EQUATION RUN ON MAPLE SOFTWARE

ENJOY...

 

Happy Springtime to all in the MaplePrimes Community! Though some in our community may not live in the northern hemisphere where flowers are beginning to bloom, many will be celebrating April holidays like Ramadan, Passover, and Easter.

One of my favorite springtime activities is decorating eggs. Today, the practice is typically associated with the Christian holiday of Easter. However, painted eggs have roots in many cultures.

For over 3,000 years, painting eggs has been a custom associated with the holiday of Nowruz, or Persian New Year, during the spring equinox. Furthermore, in the Bronze Age, decorated ostrich eggs were traded as luxury items across the Mediterranean and Northern Africa. Dipped eggs have also played an important role in the Jewish holiday of Passover since the 16th century.

To celebrate this tradition, I would like to invite all of the Maplesoft community to create a decorated egg of their own with the Easter Egg Art Maple Learn document. In this document, an ovoid egg equation is used to define the shape of an egg. 



The ovoid egg equation mimics the shape of a typical hen’s egg. Each bird species lays differently shaped eggs. For example, an ostrich’s egg is more oblong than an owl’s, and an owl’s egg is rounder than a goose’s. Surprisingly, every egg can be described by a single equation with four parameters:



Learn more about this equation and others like it with John May’s Egg Formulas Maple Learn document.

The Easter Egg Art document includes 9 different decorative elements; users can change the color, position, and size of each in order to create their own personal egg! The egg starts out looking like this:



In just a couple of minutes, you can create a unique egg. Have fun exploring this document and share a screenshot of your egg in the comments below!  Here’s one I made:


How to Create Graph Equation of Water Drop Wave in Cartesian Space using single Implicit Function only run by Maple Software

The Equation is:   z = - cos( (x2+y2)0.5 - a)  with paramater a is moving form 0 to 2pi 

Enjoy...

Plese Click link below to see full equation in Maple software

Water_Drop_Wave.mw

1 2 3 4 5 6 7 Last Page 3 of 71