## Posts Tagged ‘maple’

### Convincing maple to solve an ODE with Neumann conditions at a symbolic valued location

Friday, November 17th, 2017

I can use maple to solve a 1D second-order ODE with Dirichlet boundary conditions at symbolic-valued locations:

# Z'' = 0, Z(a)=0, Z(b) = 1
dsolve({diff(Z(r),r,r) = 0,Z(a)=0,Z(b)=1});

This correctly returns

r       a
Z(r) = - ----- + -----
a - b   a - b

I can also easily convince maple to solve this ODE with some Neumann (normal derivative) boundary conditions at at fixed-value, numeric location:

# Z'' = 0, Z(a) = 1, Z'(0) = 0
dsolve({diff(Z(r),r,r) = 0,Z(a)=1,eval(diff(Z(r),r),r=0)=0});

produces

Z(r) = 1

But if I try naively to use a Neumann condition at a symbolic value location

# Z'' = 0, Z(a) = 1, Z'(b) = 0
dsolve({diff(Z(r),r,r) = 0,Z(a)=1,eval(diff(Z(r),r),r=b)=0});

then I get an error:

Error, (in dsolve) found differentiated functions with same name but depending on different arguments in the given DE system: {Z(b), Z(r)}

After a long hunt, I found the solution. dsolve takes an optional second argument that can tell it what the dependent variable actually is. So the correct call is:

# Z'' = 0, Z(a) = 1, Z'(b) = 0
dsolve({diff(Z(r),r,r) = 0,Z(a)=1,eval(diff(Z(r),r),r=b)=0});

and this gives the correct answer

Z(r) = 1

### Energy optimization, calculus of variations, Euler Lagrange equations in Maple

Tuesday, August 16th, 2016

Here’s a simple demonstration of how to solve an energy functional optimization symbolically using Maple.

Suppose we’d like to minimize the 1D Dirichlet energy over the unit line segment:

min  1/2 * f'(t)^2
f
subject to: f(0) = 0, f(1) = 1

we know that the solution is given by solving the differential equation:

f''(t) = 0, f(0) = 0, f(1) = 1

and we know that solution to be

f(t) = t

How do we go about verifying this in Maple:

with(VariationalCalculus):
E := diff(f(t),t)^2:
L := EulerLagrange(E,t,f(t)):

so far this will output:

L := {-2*diff(diff(x(t),t),t), -diff(x(t),t)^2 = K[2], 2*diff(x(t),t) = K[1]}

Finally solve with the boundary conditions using:

dsolve({L[1],f(0)=0,f(1)=1});

which will output

t(t) = t

### Inverse of common ease curves

Wednesday, July 13th, 2016

Awkwardly I ended up needing the inverse of an ease curve (or S-curve) used in tweening animations. For trigonometric ease curves, this is easy. For example, if your ease filter is:

f = 0.5-cos(x*pi)*0.5;

then the inverse is:

x = acos(-2*(f-0.5))/pi

But if you’re using the famous cubic ease curve,

f = 3.*x.^2 - 2.*x.^3

then the relevant inverse is a complex function (involving the i = sqrt(-1)) that produces real values for f in [0,1]:

x = -1./4.*(1-2.*f+2.*(f.^2-f).^(1./2)).^(1./3)-1./4./(1-2.*f+2.*(f.^2-f).^(1./2)).^(1./3)+1./2-1./2.*i.*3.^(1./2).*(1./2.*(1-2.*f+2.*(f.^2-f).^(1./2)).^(1./3)-1./2./(1-2.*f+2.*(f.^2-f).^(1./2)).^(1./3));

I haven’t bothered to see if this can be simplified into something tidy. It’d be great to get rid of the i.

### Uninstalling Maple toolbox breaks matlab app

Tuesday, November 12th, 2013

I tried to uninstall the maple matlab toolbox via Uninstall Maple 17.app. This seem to execute OK, but then when I tried to launch Matlab the app immediately exited. Investigating the console I saw that it was producing an error:

2013-11-12 19:23:22.290 StartMATLAB[494:707] execl("/Applications/MATLAB_R2012a.app/Contents/MacOS/../../bin/matlab", "/Applications/MATLAB_R2012a.app/Contents/MacOS/../../bin/matlab", "-desktop", "-arch=maci64", nil) failed!  errno: 13

Upon closer inspection it seems that the matlab file existed but it’s permissions were not setup for executing. I fixed this by issuing:

sudo chmod a+x /Applications/MATLAB_R2012a.app/Contents/MacOS/../../bin/matlab

My guess is that maple’s script tried to replace matlab with a backup but forgot to recent the permissions.

### atan2 is harmonic

Thursday, August 22nd, 2013

Don’t you forget it. For the lazy among us, here’s a maple proof:

simplify(diff(arctan(y,x),x,x)+diff(arctan(y,x),y,y),size);

produces:

0

### Get values of solve variables as list in maple

Wednesday, August 21st, 2013

I issue a solve command in maple like this one:

s:= solve({2*x1-x2=0,x1+x2=1},{x1,y2});

and I see as output:

s := {x1 = 1/3, x2 = 2/3}

I’d like to get the results as an ordered list. I’ve already carefully chosen my variable names so that they’re ordered lexicographically, which is respected by the solve output. I tried to use convert, table, entries, indices to no avail. What I came up with is:

convert(map(rhs,s),list);

which produces:

[1/3, 2/3]

This can be easily dumped into matlab for processing.

Note: If you use square brackets in your solve:

s:= solve({2*x1-x2=0,x1+x2=1},[x1,x2]);

Then the map above will give you an error like this one:

Error, invalid input: rhs received [x1 = 1/3, x2 = 2/3], which is not valid for its 1st argument, expr

You need to get the first index of the solve output:

map(rhs,s[1]);

### Minimize a quadratic energy with linear equality constraints symbolically in maple

Wednesday, August 21st, 2013

I’m working with Maple to fit a polynomial to specific values and derivatives so that it minimized a quadratic energy. I tried the Optimization toolbox/package and came up with this:

with(Optimization):
p := (t) -> c3*t^3 + c2*t^2;
C1 := eval(p(t),t=1)-1;
C2 := eval(diff(p(t),t),t=1);
m := Minimize(int(diff(p(t),t)^2,t=0..1),{C1=0,C2=0});

but this seems to give me a numerically optimized solution:

m := [1.19999999999997, [c2 = 3., c3 = -2.]]

Rather I want a symbolically optimized solution. So far, I can do this by explicitly computing the Lagrangian and finding the saddle point:

with(VectorCalculus):
Lambda := int(diff(p(t),t)^2,t=0..1) + lambda_1*(C1) + lambda_2*(C2);

This gives me the exact solution:

[[c3 = -2, c2 = 3, lambda_1 = -12/5, lambda_2 = 1/5]]

### Trilinear interpolation in Maple

Friday, July 19th, 2013

Here’s a maple function for trilinear interpolation in a cube. I tried to match the notation in the wikipedia entry. I define C = (x,y,z) where x grows to the right, y grows up and z grows into the screen

tri := (x,y,z,C000,C001,C010,C011,C100,C101,C110,C111) ->
y*(z*(x*C000+(1-x)*C100)+(1-z)*(x*C010+(1-x)*C110)) +
(1-y)*(z*(x*C001+(1-x)*C101)+(1-z)*(x*C011+(1-x)*C111));

### Areas of quadrilaterals within a triangle determined by its circumcenter (or perpendicular bisectors), using matlab

Thursday, February 11th, 2010

First, I need to locate the circumcenter of the triangle. The circumcenter is the center of the circle circumscribing the three points of the triangle. Thus, it is equidistance from all three points.

I do not need the actual world coordinates of the circumcenter, rather it is okay for me to just know the relative distance to each corner. The law of cosines tells me that the relative distance of the circumcenter from the line opposite point A (called a in the picture) is cos A (that is cosine of the angle at A). Likewise for B and C. Knowing these I can convert to barycentric coordinates by multiplying by the opposite edge lengths.

cosines = [ ...
(edge_norms(:,3).^2+edge_norms(:,2).^2-edge_norms(:,1).^2)./(2*edge_norms(:,2).*edge_norms(:,3)), ...
(edge_norms(:,1).^2+edge_norms(:,3).^2-edge_norms(:,2).^2)./(2*edge_norms(:,1).*edge_norms(:,3)), ...
(edge_norms(:,1).^2+edge_norms(:,2).^2-edge_norms(:,3).^2)./(2*edge_norms(:,1).*edge_norms(:,2))];
barycentric = cosines.*edge_norms;
normalized_barycentric = barycentric./[sum(barycentric')' sum(barycentric')' sum(barycentric')'];

A quality of the barycentric trilinear coordinates of a point in a triangle is that they specify the relative area of inner triangles form by connecting the corners to that point. In our case, the triangles made by connecting the circumcenter to each corner. Here, the aquamarine, red, and purple triangles corresponding to points A, B, and C.

Now I find the area of the triangle and then multiple against the normalized barycentric coordinates of the circumcenter to get the areas of each of these inner triangles:

areas = 0.25*sqrt( ...
(edge_norms(:,1) + edge_norms(:,2) - edge_norms(:,3)).* ...
(edge_norms(:,1) - edge_norms(:,2) + edge_norms(:,3)).* ...
(-edge_norms(:,1) + edge_norms(:,2) + edge_norms(:,3)).* ...
(edge_norms(:,1) + edge_norms(:,2) + edge_norms(:,3)));
partial_triangle_areas = normalized_barycentric.*[areas areas areas];

Finally the quadrilateral at each point is built from half of the neighboring inner triangles: the green, blue and yellow quads corresponding to A, B, and C.

Note, that the fact that these inner triangles are cut perfectly in half follows immediately from the fact the the circumcenter is equidistance from each corner. So the last step is just to add up the halves:

quads = [ (partial_triangle_areas(:,2)+ partial_triangle_areas(:,3))*0.5 ...
(partial_triangle_areas(:,1)+ partial_triangle_areas(:,3))*0.5 ...
(partial_triangle_areas(:,1)+ partial_triangle_areas(:,2))*0.5];

Note: I’ve suffered trying to use mathematic or maple to reduce these steps to a simple algebra formula (I must be mistyping something). As far as I can tell the above steps check out (of course, it only makes sense for cases where the circumcenter is within or on the triangle: i.e. the triangle is not obtuse). I have checked out my equations on maple with the following functions:

mycos:=(a,b,c)->(b^2+c^2-a^2)/(2*b*c);
barycentric:=(a,b,c)->a*mycos(a,b,c);
normalized_barycentric:=(a,b,c)->barycentric(a,b,c)/(barycentric(a,b,c)+barycentric(b,c,a)+barycentric(c,a,b));
area:=(a,b,c)->sqrt((a+b-c)*(a-b+c)*(-a+b+c)*(a+b+c))/4;
inner_triangle:=(a,b,c)->normalized_barycentric(a,b,c)*area(a,b,c);

Which then you can run:

To see a beautiful long algebra formula for the quadrilateral at point A.
Which in matlab corresponds to:

[-sqrt(-(edge_norms(:,1) + edge_norms(:,2) + edge_norms(:,3)).* ...
(edge_norms(:,1) - edge_norms(:,2) - edge_norms(:,3)).* ...
(edge_norms(:,1) - edge_norms(:,2) + edge_norms(:,3)).* ...
(edge_norms(:,1) + edge_norms(:,2) - edge_norms(:,3))).* ...
(2.*edge_norms(:,2).^2.*edge_norms(:,3).^2  + edge_norms(:,1).^2.* ...
edge_norms(:,2).^2  - edge_norms(:,2).^4  + edge_norms(:,1).^2.* ...
edge_norms(:,3).^2  - edge_norms(:,3).^4)./ ...
(8.*(edge_norms(:,1) + edge_norms(:,2) + edge_norms(:,3)).* ...
(edge_norms(:,1) - edge_norms(:,2) - edge_norms(:,3)).* ...
(edge_norms(:,1) - edge_norms(:,2) + edge_norms(:,3)).* ...
(edge_norms(:,1) + edge_norms(:,2) - edge_norms(:,3))), ...
-sqrt(-(edge_norms(:,2) + edge_norms(:,3) + edge_norms(:,1)).* ...
(edge_norms(:,2) - edge_norms(:,3) - edge_norms(:,1)).* ...
(edge_norms(:,2) - edge_norms(:,3) + edge_norms(:,1)).* ...
(edge_norms(:,2) + edge_norms(:,3) - edge_norms(:,1))).* ...
(2.*edge_norms(:,3).^2.*edge_norms(:,1).^2  + edge_norms(:,2).^2.* ...
edge_norms(:,3).^2  - edge_norms(:,3).^4  + edge_norms(:,2).^2.* ...
edge_norms(:,1).^2  - edge_norms(:,1).^4)./ ...
(8.*(edge_norms(:,2) + edge_norms(:,3) + edge_norms(:,1)).* ...
(edge_norms(:,2) - edge_norms(:,3) - edge_norms(:,1)).* ...
(edge_norms(:,2) - edge_norms(:,3) + edge_norms(:,1)).* ...
(edge_norms(:,2) + edge_norms(:,3) - edge_norms(:,1))), ...
-sqrt(-(edge_norms(:,3) + edge_norms(:,1) + edge_norms(:,2)).* ...
(edge_norms(:,3) - edge_norms(:,1) - edge_norms(:,2)).* ...
(edge_norms(:,3) - edge_norms(:,1) + edge_norms(:,2)).* ...
(edge_norms(:,3) + edge_norms(:,1) - edge_norms(:,2))).* ...
(2.*edge_norms(:,1).^2.*edge_norms(:,2).^2  + edge_norms(:,3).^2.* ...
edge_norms(:,1).^2  - edge_norms(:,1).^4  + edge_norms(:,3).^2.* ...
edge_norms(:,2).^2  - edge_norms(:,2).^4)./ ...
(8.*(edge_norms(:,3) + edge_norms(:,1) + edge_norms(:,2)).* ...
(edge_norms(:,3) - edge_norms(:,1) - edge_norms(:,2)).* ...
(edge_norms(:,3) - edge_norms(:,1) + edge_norms(:,2)).* ...
(edge_norms(:,3) + edge_norms(:,1) - edge_norms(:,2)))];

Note: It seems also that the wikipedia page for the circumscribed circle has a direct formula for the barycentric coordinates of the circumcenter. Though, wolfram’s mathworld has an ostensibly different formula. Upon a quick review Wolfram’s is at least easily verifiable to be the same as mine (recall these only need to be relative to each other).