Maximum Calculus with Maxima
We looked at Maxima in the February 2011 issue to do algebra and rearrange some equations. But those aren't the only tricks up Maxima's sleeve. This month, I describe how Maxima can help with differential equations, but I'm going to leave out some of the intermediate results to save some space.
A lot of science involves figuring out how systems change over time and what causes those changes. When you start looking at changes, and especially rates of change, that is essentially calculus. Calculus and rates of change also are linked to slopes of lines on graphs. When you plot data and find an equation that describes it, you can find the slope of the line by taking the derivative of the equation. Let's look at a falling object and see what theory has to say about it.
You should start by looking at how you get a derivative. Let's say you have the equation:
(%i1) f(x):= 2 + x^2; 2 (%o1) f(x) := 2 + x
You would find the first derivative by calling the function diff, giving it the equation to differentiate along with the variable to differentiate by. So, you would write:
(%i2) answer:diff(f(x),x); (%o2) 2x
Maxima can do differentiation of expressions too. If you have a couple equations, you can derive their ratio with:
(%i3) g(x):= x^(1/2); (%i4) ratio_diff:diff(g(x)/f(x),x); 3/2 1 2 x (%o4) ----------------- - ----------- 2 2 2 2 sqrt(x) (x + 1) (x + 1)
This might be a bit messy to work with, so you might want to refactor it to a more concise form:
(%i5) factor(ratio_diff); 2 3 x - 1 (%o5) - ------------------ 2 2 2 sqrt(x) (x + 1)
Maxima also can handle trigonometric functions, but there are lots of identities you can use to help simplify equations with trig functions in them. By default, Maxima does not try to apply these unless you specifically say so, using special functions. As an example, let's say you have the following equation:
(%i6) diff(sin(x)/(1 + cos(x)),x); 2 sin (x) cos(x) (%o6) ------------- + ---------- 2 cos(x) + 1 (cos(x) + 1) (%i7) factor(%); 2 2 sin (x) + cos (x) + cos(x) (%o7) -------------------------- 2 (cos(x) + 1)
That's still not very simple. If you then apply the function trigsimp, you can force Maxima to apply trigonometric simplification rules to the equation and see what you get:
(%i8) trigsimp(%); 1 (%o8) ---------- cos(x) + 1
You should be aware of some important caveats regarding how Maxima treats trig
functions. The first is that
sin(x)^(-1) is the reciprocal of sine, not
arcsine. To get the arcsine, you would use
asin(x). The other is another
trig simplification function, trigreduce. This function is used to reduce
the powers of trig functions by using the multiple angle formulas. For
(%i9) trigsimp(cos(x)^2 + 2*sin(x)^2); 2 (%o9) sin (x) + 1 (%i10) trigreduce(cos(x)^2 + 2*sin(x)^2); cos(2 x) + 1 1 cos(2 x) (%o10) ------------ + 2 (- - --------) 2 2 2
That may not look simpler than what you would get from trigsimp, but it is an easier form of the equation to use with other functions, like integration.
Maxima can apply the chain rule when doing a derivative. Say you have the equation:
(%i11) f(x):= x^3); 3 (%o11) f(x) := x (%i12) depends(x,u)$ (%i13) diff(f(x),u); 2 dx (%o13) 3 x -- du
The line at
%i12 uses a new function,
depends. This is a way of telling
Maxima that x is a function of u, without explicitly defining a function
describing this relationship. If you decide later that you want to define an
actual equation for this relation, you can use:
(%i14) remove([x,u],dependency); (%o14) done (%i15) x:sin(u); (%o15) sin(u) (%i16) diff(f(x),u); 2 (%o16) 3 cos(u) sin (u)
Along the same lines, Maxima can handle implicit differentiation. Say you have the equation x^2 + y^2 = 25, and you want to find dy/dx. You need to use the depends function I just mentioned to handle this:
(%i17) eqn := x^2 + y^2 = 25; 2 2 (%o17) y + x = 25 (%i18) depends(y,x); (%o18) [y(x)] (%i19) deriv_of_eqn:diff(eqn,x); dy (%o19) 2 y -- + 2 x = 0 dx (%i20) solve(deriv_of_eqn,'diff(y,x)); dy x (%o20) [-- = - -] dx y
The other side of calculus is integration. The basic function to do that in
Maxima is called
integrate. This function can do both definite and
indefinite integrals. Indefinite integrals are the symbolic form of
integration you likely learned in school. For example:
(%i21) integrate(x^2,x); 3 x (%o21) - 3
A definite integral actually is evaluated over an interval. This form of an integral can be visualized as the area under the curve defined by the equation you are integrating. To do definite integrals, simply add two arguments giving the start and end points of the interval:
(%i22) integrate(x^2,x,0,1); 1 (%o22) - 3
Putting all these techniques together, you can solve a differential
equation for a given variable—for example, solve dy/dx = f(x) for y. You can do
this by doing all the required algebra and calculus, but you don't
really need to. Maxima has the very powerful function,
ode2, which can do
it in one step. Start with your equation:
(%i23) eq: 'diff(y,x) = sqrt(1/x^2 - 1/x^3); dy 1 1 (%o23) -- = sqrt(-- - --) dx 2 3 x x (%i24) ode2(eq,y,x); 2 2 2 sqrt(x - x) (%o24) y = log(2 sqrt(x - x) + 2 x - 1) - ------------- + %c x
This one function call does the integration and the solve steps and gives you a final answer to the differential equation.
Let's say you're doing an experiment dropping a coin and timing how long it takes to fall. How do you know whether the times you are measuring actually make sense? Let's start with the most basic law: force = mass * acceleration.
The mass of the coin is a constant, so ignore that for now. The force is the force due to gravity, pulling the coin down to the ground, and the acceleration describes the coin's motion due to this force. The force due to gravity is a constant, at least here on Earth, and it depends linearly on the mass, so you can define the force as:
(%i1) force: mass * g; (%o1) g mass
The acceleration also is a constant, because both the mass and the force are constants. Acceleration is simply the rate of change of the velocity, and the velocity is the rate of change of the position, so you can set that up as:
(%i2) depends(y,t); (%o2) [y(t)] (%i3) acceleration: 'diff('diff(y,t),t); 2 d y (%o3) --- 2 dt
Putting it all together, you get:
(%i4) eq_of_motion: force = mass * acceleration; 2 d y (%o4) g mass = mass --- 2 dt (%i5) solve(eq_of_motion, y); 2 d y (%o5) [--- = g] 2 dt
You can see right away that how fast an object falls doesn't depend on the mass at all. Galileo was right! The next step is to do some integrating and see what you end up with:
(%i6) integrate(%,t); dy (%o6) [-- = g t + %c1] dt
At this step, you would be able to find out the velocity (dy/dt) at time t.
The additional term
%c1 is a constant of
integration. In this case, you can
see that it represents the initial velocity of your penny. One more round of
integration gives this:
(%i7) integrate(%,t); / 2 [ dy g t (%o7) [I -- dt = ---- + %c1 t + %c2] ] dt 2 /
Now you can find the position, y, of your coin at any time, t. Again, a new
constant of integration is introduced,
%c2. In this
case, you can see that
this represents the starting height of your coin. But that's not what you
were measuring. You were measuring how long it took the coin to drop a given
distance. So you need to do a bit of rearranging. Because you are dropping
coin, you know that the start velocity is 0 (that is,
%c1=0). You can rewrite
things a little to make it a bit clearer:
(%i8) eqn: y = (g * t^2)/2 + y0; 2 g t (%o8) y = y0 + ---- 2 (%i9) solve(eqn,t); y y0 y y0 (%o9) [t = - sqrt(2) sqrt(- - --), t = sqrt(2) sqrt(- - --)] g g g g
There you go. You now have an equation for the time, given a height that your coin is dropping. With this theoretical relation under your belt, you can check to see whether gravity is working correctly in your local lab. If not, you should contact the Nobel committee straightaway.
This only scratches the surface of Maxima's capabilities in dealing with calculus and differential equations, but hopefully, this article gives you a starting point. Happy integrating.
Joey Bernard has a background in both physics and computer science. This serves him well in his day job as a computational research consultant at the University of New Brunswick. He also teaches computational physics and parallel programming.
Practical Task Scheduling Deployment
July 20, 2016 12:00 pm CDT
One of the best things about the UNIX environment (aside from being stable and efficient) is the vast array of software tools available to help you do your job. Traditionally, a UNIX tool does only one thing, but does that one thing very well. For example, grep is very easy to use and can search vast amounts of data quickly. The find tool can find a particular file or files based on all kinds of criteria. It's pretty easy to string these tools together to build even more powerful tools, such as a tool that finds all of the .log files in the /home directory and searches each one for a particular entry. This erector-set mentality allows UNIX system administrators to seem to always have the right tool for the job.
Cron traditionally has been considered another such a tool for job scheduling, but is it enough? This webinar considers that very question. The first part builds on a previous Geek Guide, Beyond Cron, and briefly describes how to know when it might be time to consider upgrading your job scheduling infrastructure. The second part presents an actual planning and implementation framework.
Join Linux Journal's Mike Diehl and Pat Cameron of Help Systems.
Free to Linux Journal readers.Register Now!
- SUSE LLC's SUSE Manager
- My +1 Sword of Productivity
- Murat Yener and Onur Dundar's Expert Android Studio (Wrox)
- Non-Linux FOSS: Caffeine!
- Managing Linux Using Puppet
- Doing for User Space What We Did for Kernel Space
- Tech Tip: Really Simple HTTP Server with Python
- Parsing an RSS News Feed with a Bash Script
- SuperTuxKart 0.9.2 Released
- Google's SwiftShader Released
With all the industry talk about the benefits of Linux on Power and all the performance advantages offered by its open architecture, you may be considering a move in that direction. If you are thinking about analytics, big data and cloud computing, you would be right to evaluate Power. The idea of using commodity x86 hardware and replacing it every three years is an outdated cost model. It doesn’t consider the total cost of ownership, and it doesn’t consider the advantage of real processing power, high-availability and multithreading like a demon.
This ebook takes a look at some of the practical applications of the Linux on Power platform and ways you might bring all the performance power of this open architecture to bear for your organization. There are no smoke and mirrors here—just hard, cold, empirical evidence provided by independent sources. I also consider some innovative ways Linux on Power will be used in the future.Get the Guide