Sign in to follow this  
Sirisian

Does someone have Matlab or Maple and mind crunching some numbers

Recommended Posts

Sirisian    2263
(sorry if this is in the wrong place, not sure if the math forum allows this. I'm desperate) Okay I have a parametric equation. X = Cos(polygon1RotRate * t + polygon1RotInit) * vertexX - Sin(polygon1RotRate * t + polygon1RotInit) * vertexY + (polygon1VelX - polygon2VelX) * t + polygon1X Y = Sin(polygon1RotRate * t + polygon1RotInit) * vertexX + Cos(polygon1RotRate * t + polygon1RotInit) * vertexY + (polygon1VelY - polygon2VelY) * t + polygon1Y integrate(sqrt(derivate(X,t)^2 + derivate(Y,t)^2),t); That is the arc length that I need. Only problem and anyone trying to do this will realize it. It will take a lot of time (a few hours or a day) to do and a lot of file paging (so increase file paging to the max if attempting). I would do it myself but when trying to use my companies matlab discs I realized they were scratched by whoever used them last. If someone that has a beefy computer and doesn't mind letting it run for a long time processing this I'd be very grateful. (it's for a small physics demo I'm coding and this is a key equation I need).

Share this post


Link to post
Share on other sites
mumpo    534
I'm taking a crack at it in Mathematica 5.2 right now. I have 4 GB of RAM and am running 64-bit Vista, so hopefully my system won't have any memory problems. My processor is a relatively tame 2.13 GHz Core 2 Duo, so who knows how long it will take. . . I'll let you know if it ever finishes (no guarantees).

Share this post


Link to post
Share on other sites
Bob Janova    769
I may be missing something but this looks like it can be done analytically, at least part of the way.

I'm going to restate variables to make me able to see what's going on so
polygonnVelX/Y => vX/Yn
polygon1X/Y => sX/Y
polygon1RotRate => r
polygon1RotInit => r0
vertexX => x
vertexY => y
So
X = x cos(rt+r0) - y sin(rt+r0) + t(vX1-vX2) + sX
Y = x sin(rt+r0) + y cos(rt+r0) + t(vY1-vY2) + sY

Assuming all the terms are not variant in t

dX/dt = xr sin (rt+r0) + yr cos(rt+r0) + (vX1-vX2)
dY/dt = xr cos (rt+r0) + yr sin(rt+r0) + (vY1-vY2)

Just setting R = rt+r0
and ∆vX/Y = (vX/Y1-vX/Y2) for ease of typing:

(dX/dt)² = x²r²sin²R + y²r²cos²R + ∆vX² + 2[xyr²cos R sin R+yr∆vX cos R+xr∆vXsin R]
and similarly for Y
(dY/dt)² = x²r²cos²R + y²r²sin²R + ∆vY² + 2[xyr²cos R sin R+yr∆vY sin R+xr∆vYcos R]

=> (dX/dt)²+(dY/dt)² =
x²r²(sin²R+cos²R) + y²r²(sin²R+cos²R) + ∆vY² + ∆vX² + 2[xyr²cos R sin R+yr∆vX cos R+xr∆vXsin R+xyr²cos R sin R+yr∆vY sin R+cos R]
= r²(x²+y²) + ∆vY² + ∆vX² + 2r cos R[y∆vX + x∆vY] + 2r sin R[x∆vX + y∆vY] + 2xyr²cos R sin R
As R there is the only thing that depends on t, you can park [r²(x²+y²) + ∆vY² + ∆vX²] in a constant (A) and you get
= A + 2r(cos R[y∆vX + x∆vY] + sin R[x∆vX + y∆vY] + xyr cos R sin R)

Hmm, okay, it looks like that as far as you go analytically, unless I'm missing things. However it shouldn't take that long to integrate that numerically?

Of course if other parts depend on T then I just wasted 10 minutes <.<

Share this post


Link to post
Share on other sites
Sirisian    2263
Thanks for trying it mumpo. I had someone else that started doing it yesterday but had to stop it after 50 minutes because he had something to do. That's the only reason I know that it will take a long time.

And yes, Bob Janova, I too thought I could do it by hand but a computer has a much higher chance of telling if it is possible.

Share this post


Link to post
Share on other sites
alvaro    21246
I entered (-Sin[w x+p]*w*p1-Cos[w x+p]*w*p2+v1)^2+(Cos[w x+p]*w*p1-Sin[w x+p]*w*p2+v2)^2 into http://integrals.wolfram.com/index.jsp and out came a formula. You should be able to use that.

EDIT: Ooops! Nevermind. I forgot the Sqrt[]. It doesn't work once you add it.

Share this post


Link to post
Share on other sites
mumpo    534
Just a quick update -- my computer is still at it. Mathematica is using about 800 MB now. I wonder if the computer can actually solve it, or if it is just going ever deeper into some infinitely recursive algorithm? Either way, I'll try to give it at least a couple more days before I give up on it.

Share this post


Link to post
Share on other sites
Sirisian    2263
Quote:
Original post by mumpo
Just a quick update -- my computer is still at it. Mathematica is using about 800 MB now. I wonder if the computer can actually solve it, or if it is just going ever deeper into some infinitely recursive algorithm? Either way, I'll try to give it at least a couple more days before I give up on it.

Thanks :) I can't tell you how much this is going to help me if this works. If not then substituting sine and cosine with a taylor series or something will be my next guess and just solve on a range. If it's the sin and cos that actually cause the problem.

I'm gonna find out next Tuesday where my university's cluster server is. I heard they use it for astronomy, maybe I can get someone to plug in equations.

I might not have made this clear. After finding the arc length (if it exists) I need to solve it for t. So that given an arc length it returns the t value associated with that length.

Share this post


Link to post
Share on other sites
daftasbrush    252
From My Understanding...I don't think there's definitive Analytical solution.

MatLab (or whatever) may kick out a formula which will initially appear to be "analytical" but I reckon it's going to contain a "numerical approximation" component....

Here's how I arrived at this conclusion...
Assuming the same as Bob ... that most everything is constant wrt t
Your problem boils down to a point rotating about a centre, and that centre is moving.

Ultimately you are trying to compute s = integral(V(t) dt)

It seems to me that you can find V(t) much more simply if you don't "decompose" everything into 2 independant cartesian axes, but consider the problem "as a whole".
Consider the 2 motion components (Rotational and Linear) as vectors, and simply add them and determine the length of the result.

The point is rotating at a constant (angular) rate...
it's distance from the "centre" r = sqrt(VertexX^2 + VertexY^2)
so it's rotational speed Vr = r * polygon1RotRate
(assuming polygon1RotRate in radians per "Unit Time")

It's Linear Velocity Vl = sqrt((polygon1VelX - polygon2VelX)^2 + (polygon1VelY - polygon2VelY)^2))

Now all you need is "Cosine rule" to determine the actual speed.

V(Theta)^2 = Vl^2 + Vr^2 + 2.Vl.Vr.cos(Theta)
so V(Theta) = sqrt(Vl^2 + Vr^2 + 2.Vl.Vr.cos(Theta))

Theta = Angle between rotational and linear velocity
= polygon1RotRate * t + polygon1RotInit - "Angle of Linear Velocity" (ALV)

Simplifying for clarity :
a = Vl^2 + Vr^2
b = 2.Vl.Vr
c = polygon1RotInit - ALV
f = polygon1RotRate

we get V(t) = sqrt(a + b*cos(c + f*t))

Now if you drop that into the Wolfram Web Integrator you get.... a result.
I imagine Matlab will give a similar result (eventually...)

(2 Sqrt[a + b Cos[c + f t]] *

c + f t 2 b
EllipticE[-------, -----]) /
2 a + b

a + b Cos[c + f t]
(f Sqrt[------------------])
a + b

But the (one?) problem is.... The Elliptic Integral
Quote:
From Wiki - "elliptic integrals originally arose in connection with the problem of giving the arc length of an ellipse."

Not (I think) a coincidence it's come up for this problem...

As I understand them (admittedly in limited fashion, so I'm not 100% on this bit), there is no exact analytical way to evaluate an elliptic integral, it has to be done numerically or by approximation (there are a few different methods)....
But either way... it doesn't look like you can get an exact answer.

So the question is....
How accurate vs complex do you want it to be....

As I see it, you could...
i) Do the original integration numerically. (using whatever method)
or
ii)Use the "analytical" integral above.. that includes a numerical integration or an approximation.

I can't say which would be more accurate....

On a happy Note:
MatLab has a function for 2nd form elliptic integrals under GNU license... might make things a bit easier!?

HTH

[Edited by - daftasbrush on March 21, 2008 1:05:49 PM]

Share this post


Link to post
Share on other sites
mumpo    534
Quote:
Original post by Sirisian
Quote:
Original post by mumpo
Just a quick update -- my computer is still at it. Mathematica is using about 800 MB now. I wonder if the computer can actually solve it, or if it is just going ever deeper into some infinitely recursive algorithm? Either way, I'll try to give it at least a couple more days before I give up on it.

Thanks :) I can't tell you how much this is going to help me if this works. If not then substituting sine and cosine with a taylor series or something will be my next guess and just solve on a range. If it's the sin and cos that actually cause the problem.

You're welcome :) . I'm not very optimistic about the current calculation ever finishing, though; I can't imagine why it would need so much memory to solve a relatively concise symbolic problem unless the answer has an infinite number of terms or it just can't do it, but won't admit it. Let me know if you want me to try something else, instead.
Quote:
I'm gonna find out next Tuesday where my university's cluster server is. I heard they use it for astronomy, maybe I can get someone to plug in equations.

I might not have made this clear. After finding the arc length (if it exists) I need to solve it for t. So that given an arc length it returns the t value associated with that length.

My calculus is quite rusty, and I never took all that much of it to begin with, but it seems to me that if my computer ever comes up with a sensible result for that integral, solving it for t shouldn't take too long.

Share this post


Link to post
Share on other sites
Sirisian    2263
Quote:
Let me know if you want me to try something else, instead.

Okay try a taylor series using power of degree 9 for sine and degree 8 for cosine

cos(x) = 1 - x^2/2 + x^4/24 - x^6/720 + x^8/40320
sin(x) = x - x^3/6 + x^5/120 - x^7/5040 + x^9/362880

I'm not sure what the hangup is. If it's the sin or cosine the sqrt might be able to work with those substitutions.

daftasbrush, I'm looking for an approximation or anything that I can use to find the arclength of that moving vertex. A vertex will never rotate more than 180 degrees in a timestep so I picked those two equations above.

(I don't need to solve for t, not sure what I was thinking. Been writing up the algorithms for a unit test. Non-rotating seems to be fine.)

Share this post


Link to post
Share on other sites
daftasbrush    252
I guess it depends on how "approximate" you want your approximation to be,

The "Analytical" Integral (in my post above), will likely yield results accurate to the precision of whatever float type you are using.. even though it's "strictly" an approximation.

But providing you remain within your stated limits : d(Theta) < pi
Numerical Integration of the equation I gave
Quote:
a = Vl^2 + Vr^2
b = 2.Vl.Vr
c = polygon1RotInit - ALV
f = polygon1RotRate
V(t) = sqrt(a + b*cos(c + f*t))

even using something as simple as Simpsons rule will likely give less than 0.5% error, in all but the most "extreme" cases.

Those "extreme" cases being.. Vr ~ Vl and d(theta) > pi/2

In those cases you could test whether (c + ft) spans pi.
If it does, then breaking the numerical integration into 2..
t=starttime to t(c+ft = pi) and t(c+ft= pi) to t = endtime
should improve things greatly.

As I said before, which method you choose depends how accurate you want it to be.

Share this post


Link to post
Share on other sites
Sirisian    2263
Well I want accuracy over speed with this physics simulation. The idea behind the arclength is that I'd be using it to determine a per vertex accuracy value. So basically it has to be fairly accurate in order for the simulation to run well.

I'm not sure if I understand what your saying 100%. When you say using simpson's rule are you suggesting that I iterate and calculate the distance at intervals? I'm looking for one equation that when given the values will give an approximation for t. I'd rather not do any numerical approximation to get this value. But if it comes down to it I'll have to.

It seems like these equations are very difficult to work with. I mean I tried last week to find out if it was possible to get an approximation of the time t when a rotating and moving vertex hit a line segment but the equations were unsolvable (even with approximations). I'm gonna talk to my math professor to see if he has any other ideas though.

Share this post


Link to post
Share on other sites
daftasbrush    252
Quote:
I'm looking for one equation that when given the values will give an approximation for t.

That's what the Analytical Integral is!

If you calculate the following values... (per my previous post)
a = Vl^2 + Vr^2
b = 2.Vl.Vr
c = polygon1RotInit - ALV
f = polygon1RotRate

Then, like any normal integral, evaluate this function for t = starttime and t = endtime.

c + f t 2 b
(2 Sqrt[a + b Cos[c + f t]] * EllipticE[-------, -----])
2 a + b
/
a + b Cos[c + f t]
(f Sqrt[------------------])
a + b


Subtract the 2 results and you've got the arc length.

You'll just need to get hold of a function to calculate EllipticE.
Like I said there's a MatLab one available for download under GNU License, there are probably others out there.
Once you've got it, just treat it like any other "black box" function like sin, cos, sqrt etc.

That'll get you the most accurate answers, and shouldn't take much time to compute.

Share this post


Link to post
Share on other sites
Sirisian    2263
I seem to be having a problem finding this EllipticE matlab program online. Is it on their site? I thought maybe it was the symbolic math toolbox but that requires a license. Hmm.

//edit, yeah I can't seem to find the equation. And this is the second kind right?
I found this link talking about them, but not much else. http://reference.wolfram.com/mathematica/ref/EllipticE.html

[Edited by - Sirisian on March 22, 2008 11:05:01 PM]

Share this post


Link to post
Share on other sites
daftasbrush    252
The MatLab one i'd found turns out to be written in (for) MatLab, so probably won't be very useful to you. but here's the link anyway.
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8805&objectType=file
Quote:
And this is the second kind right?
Yup, strictly the Incomplete Elliptic Integral of the Second Kind.

Your best bet is probably Stephen Moshier's Cephes Library.
http://www.moshier.net/
The whole package is available for download, cephes28.zip
but you pretty much just need the ellie.c file.

There's serveral versions of the file in the package, for different precisions
Single, double, Ldouble etc... take your pick.

His "license" statement is a bit wooly... (see the readme in the zip file)
but it includes the line...
Quote:
What you see here may be used freely but it comes with no support or guarantee.


I've seen the library billed as "Free" .. somewhere..
And it's used by several other packages (Grace, Labplot) which are themselves under GPL.

I you feel you need futher clarification on useage / licensing..
shoot the guy an email.

[Edited by - daftasbrush on March 23, 2008 8:12:31 AM]

Share this post


Link to post
Share on other sites
mumpo    534
Quote:
Original post by Sirisian
Quote:
Let me know if you want me to try something else, instead.

Okay try a taylor series using power of degree 9 for sine and degree 8 for cosine

cos(x) = 1 - x^2/2 + x^4/24 - x^6/720 + x^8/40320
sin(x) = x - x^3/6 + x^5/120 - x^7/5040 + x^9/362880

I'm not sure what the hangup is. If it's the sin or cosine the sqrt might be able to work with those substitutions.

daftasbrush, I'm looking for an approximation or anything that I can use to find the arclength of that moving vertex. A vertex will never rotate more than 180 degrees in a timestep so I picked those two equations above.

(I don't need to solve for t, not sure what I was thinking. Been writing up the algorithms for a unit test. Non-rotating seems to be fine.)

I ran it with the taylor series approximations. Mathematica finished thinking after a while with those, but it was unable to actually solve the integral. I suggest you focus on the method daftasbrush was kind enough to work out for you. I'll keep an eye on this thread for a bit longer, though, if you have any more ideas.

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this