JacquesCarette / Drasil

Generate all the things (focusing on research software)
https://jacquescarette.github.io/Drasil
BSD 2-Clause "Simplified" License
142 stars 26 forks source link

Feedback on DblPendulum #2275

Closed smiths closed 1 year ago

smiths commented 3 years ago

@oluowoj, I've reviewed the physics behind the DblPendulum example and I have the following feedback:

smiths commented 3 years ago

I've been thinking about the pendulum problem further. I think the work will be much easier if you switch to a polar coordinate system. If you continue on the path of Cartesian coordinates, you will get two ODEs (one for ax and one for ay). To get one ODE, we would have to convert to a polar coordinate system at this point, so why not just start there.

The web-page that you used for the figure of the physical system shows how to derive the equation of motion for a pendulum using a polar coordinate system. This switch will greatly simplify your document.

oluowoj commented 3 years ago

@smiths, for the derivation of the equation of pendulum motion, the equation is a 2nd order derivative image

I don't think we have any functions in Drasil to write a 2nd derivative equation.

@JacquesCarette, any thoughts about a workaround?

smiths commented 3 years ago

@oluowoj, we discussed this problem in a previous meeting. There isn't currently a way to solve second order ODEs in Drasil. However, Drasil is close to being able to solve problems like this. Naveen (CAS 741 student) has said that he will look into it.

My guess is that you don't want to tackle adding second order ODEs to Drasil, and I know you are anxious to have a complete example. Why don't we "cheat" a little, just temporarily, and use the first derivative (as you show in your SRS above), even though we know it is really the second derivative? We can go back and fix is later, but for now, you want to be able to complete your example.

I looked at solving the second order ODE as two first order ODEs (a standard solution technique), but this doesn't work for us because the two first order ODEs are coupled. We would need to extend Drasil to solve a system of first order ODEs. Drasil is close to being able to this as well, Brooks has the infrastructure in there, but it hasn't been completed.

oluowoj commented 3 years ago

Yes, we did discuss this. I would just keep the equation as is. I was trying some hacks initially. Could you please point me to Brooks's instructions to generate code in Drasil?

smiths commented 3 years ago

Brooks didn't write instructions, as far as I know, he wrote code for interfacing to an external library, like an ODE solver library. He had the input ODE as a list of ODEs (as opposed to a single ODE) to facilitate future expansion to a system of ODEs. This is what I meant when I said infrastructure.

JacquesCarette commented 3 years ago

You can simply do deriv (deriv (sy XXX) time) time where XXX is the symbol you want to get a second derivative. It won't necessarily display very well right now, but it is still correct. We can fix things later, but the knowledge should be properly encoded.

smiths commented 3 years ago

@oluowoj, I thought the question was about generating code to solve second order equations. I didn't realize the question was about displaying the correct equation. Please use the suggestion from @JacquesCarette for displaying the equation. However, we still need to solve the problem of generating code to solve second order (and higher) ODEs.

oluowoj commented 3 years ago

@JacquesCarette, I am getting this error when I build the pendulum project. Any thoughts?

image

I just pushed the code to pendulumDerv branch

JacquesCarette commented 3 years ago

Oh my. <<loop>> is one of the hardest things to debug in Haskell. It means it has detected that something has gone into an infinite loop.

This usually means that you've got something defined in terms of itself, in a way that you really didn't mean. Luckily, it should be the case that the problem was introduced in your latest changes. So I would back them out, and re-introduce them as slowly as possible, to know which line is the cause.

oluowoj commented 3 years ago

Right, it was a case of input and input constraint mismatch.

initialPendAngleCons = constrained' initialPendAngleCons [gtZeroConstr] (dbl 2.1)

instead of

initialPendAngleCons = constrained' initialPendAngle [gtZeroConstr] (dbl 2.1).

JacquesCarette commented 3 years ago

Glad it wasn't too hard to spot and fix!

balacij commented 1 year ago

@cd155 did some work on DblPendulum, including, for example, switching it to an actua double pendulum. We'll just need to verify which feedback is still relevant after @cd155's changes.

smiths commented 1 year ago

I've checked through this issue and the feedback has either been addressed, or is no longer relevant because we have changed the pendulum and double pendulum examples. We can close this issue, although I will create another for a review of the current double pendulum example.