NTNUJAVA Virtual Physics LaboratoryEnjoy the fun of physics with simulations! Backup site http://enjoy.phy.ntnu.edu.tw/ntnujava/
October 22, 2018, 12:10:45 am

You cannot always have happiness but you can always give happiness. ..."Mother Teresa(1910-1997, Roman Catholic Missionary, 1979 Nobel Peace Prize)"

 Pages: [1]   Go Down
 Author Topic: I have a problem with Programming Simultaneous Coupled Complicated ODEs  (Read 11816 times) 0 Members and 1 Guest are viewing this topic. Click to toggle author information(expand message area).
Z-D
Newbie

Offline

Posts: 3

 « Embed this message on: February 17, 2009, 11:09:45 am » posted from:Kuala Lumpur,Kuala Lumpur,Malaysia

Hi all, hope you can help.

I have seen examples similar to my problem, for example the Double Pendulum.

For the Double Pendulum, in the Evolution Tab, the last two ODEs are given to be equal to a function (Attachment 1). And the functions are defined in the Custom tab (Attachment 2).

Due to the complicated non-linearity of the coupled simultaneous ODEs whichI am trying to programme, I am hoping to find a way to define all the derivatives in the Evolution Tab without defining a function (Attachment 3), and later writing out the equation of motion to define the system (Attachment 4). Attachment 4 is my attempt to do so, but I tried doing so in the Fixed Varaibles tab.
Right now it is not working, but it could be simply because the system is too complicated or they could be some mistakes in the equation.

But I would like to ask if what I am attempting, could be done. Perhaps there are example out there which I am no aware of which could help me.

Attachment 5 is EJS file.

... sorry, I actually have more to write right now, but I have to go...

Hope you can help. And thanks in advanced.
 Double Pendulum - Evolution.jpg (65.17 KB, 600x598 - viewed 371 times.)  Double Pendulum - Custom.jpg (81.66 KB, 1043x598 - viewed 381 times.)  Project Model - Evolution.jpg (82.38 KB, 618x739 - viewed 390 times.)  Project Model - Fixed Relations.jpg (87.45 KB, 618x739 - viewed 391 times.) *** There are 1 more attached files. You need to login to acces it! Logged
lookang
Hero Member

Offline

Posts: 1787

http://weelookang.blogspot.com

 « Embed this message Reply #1 on: February 17, 2009, 02:19:47 pm » posted from:SINGAPORE,SINGAPORE,SINGAPORE

could you upload your file that works first, it will help me to understand and appreciate your question better.

currently, it does not compile.

glad you are also using EJS. We are not alone! welcome to this community

variable double theta_A_double_dd

commented
//theta_A_dd + (y_dd+g)/l_p*Math.sin(theta_A) + l*phi_dd/l_p*Math.sin(phi+theta_A) + (l/l_p)*phi_d^2*Math.cos(phi+theta_A) = 0;
//theta_B_dd + (y_dd+g)/l_p*Math.sin(theta_B) - l*phi_dd/l_p*Math.sin(phi+theta_B) + (l/l_p)*phi_d^2*Math.cos(phi+theta_B) = 0;
// y_dd + 2*k/(M+2*m)*y + m*l_p/(M+2*m)*(theta_A_dd*Math.sin(theta_A) + theta_A_d^2*Math.cos(theta_A) + theta_B_dd*Math.sin(theta_B) + theta_B_d^2*Math.cos(theta_B)) + g = 0;
//phi_dd + g/l*Math.cos(phi) + k/m*Math.sin(phi)*Math.cos(phi) + (theta_A_dd*Math.sin(phi+theta_A) + l_p/(2*l)*theta_A_d^2*Math.cos(phi+theta_A) - theta_B_dd*Math.sin(phi+theta_B) - theta_B_d^2*Math.cos(phi+theta_B)) = 0;

i suspect you can't assign
theta_A_dd + (y_dd+g)/l_p*Math.sin(theta_A) + l*phi_dd/l_p*Math.sin(phi+theta_A) + (l/l_p)*phi_d^2*Math.cos(phi+theta_A) = 0;

the programming syntax says to assign a value of 0 to theta_A_dd + (y_dd+g)/l_p*Math.sin(theta_A) + l*phi_dd/l_p*Math.sin(phi+theta_A) + (l/l_p)*phi_d^2*Math.cos(phi+theta_A) .

that's is all i can help i don't know how to solve your problem yet!
 *** There are 1 more attached files. You need to login to acces it! « Last Edit: February 17, 2009, 02:33:14 pm by lookang » Logged
Fu-Kwun Hwang
Hero Member

Offline

Posts: 3080

 « Embed this message Reply #2 on: February 17, 2009, 03:57:51 pm »

Quote
I am hoping to find a way to define all the derivatives in the Evolution Tab without defining a function

The error will become larger and larger if you try to put function in the "Fixed relations" page instead of define a function.

You might need to have some background about numerical method in order to fully understand it.
I will try to explain it starts with how EJS works:
1. EJS will read "Variables" page: it will ask computer to allocate memory for each variables defined and assing value to the memory. for example x=5.; it will put 0.05 to memory assigned for x.
2. EJS will read "Initialization" page: This is an optional page. EJS will do whatever java command
3. EJS will read "Evolution" page: This is the page to move time step forward and cacculate variable value at the next time step. EJS will use user selected "solver" method to solve the ODE.
4. EJS will read "Fixed relation" page: The system has been moved from time t to t+dt. And java code in this page will be executed. But all the command execute here is done after step 3.
5. EJS will show all the GUI elements according to value assigned to element's properties.
6. EJS will go back to step 3 if it in play mode and continue process 3-4-5 cycles, otherwise it will stop.

Let me use a simple hormonic motion as examples to explain why relations for derivatives can not be defined in the "Fixed relation page".
If the equestion we try to solve is
F=m*d2x/dt2= -k*x-b*vx;
we will change the above second order differtial equation into two first order differential equations:
dx/dt=vx;
dvx/dt=-k/m*x-b/m*vx;

You need to select Runge-Kutta (4th) order to reduce the error.
If the Euler's method were selected , the error will accumuate will soon.
The middle method will reduced error to 10% of the Euler's method .
The error for Runge-Kutta (4th) order is 0.00% of the Euler's method.
And this is the best method for most of the cases.

If you change the above equestions to
dx/dt=vx;(STEP 3 )
dvx/dt=ax; (STEP 3)
and define ax=-k/m*x-b/m*vx; at fixed relation page (STEP 4)
The error will be the same as Euler's method even if Runge-Kutta (4th) order was selectd.

Because EJS use the above equation to move time from t to t+dt (STEP3 to STEP4)
But vx(t) is not the same as vx(t+dt), and x(t) is not the same as x(t+dt)

For Euler's method :
x(t+dt)=x(t)+vx*dt;
vx(t+dt)=vx(t)+k/m*x(t)-b/m*x(t);

For midpoint method:
Two time steps(t+dt/2 and t+dt) were calculated to move from t-> t+dt.

For Runge-kutta 4th order method:
Four time steps were calculated to reduced the error to 1/10000 as Euler's method.
The right hand side need to be evaluate 4 times. If the correct formula is only available at STEP4.
The right hand side will always be the same , the Runge-kutta 4th order method is not calculated with correct value. The error move back to the same as Euler's method.

You need to read document about Runge-kutta 4th order method to fully understand how the numerical method works and why it can reduced the error.

You need to define function or write the equations into the right hand side of the evolution page, if the right hand side contain variables which is also in the left hand side (like x and vx in the above examples.).
If the right hand side did not contain any variables at the left hand side, like k,b,m are all constant . Then, you can defined those at fixed relation page.
 Logged
lookang
Hero Member

Offline

Posts: 1787

http://weelookang.blogspot.com

 « Embed this message Reply #3 on: February 17, 2009, 04:53:35 pm » posted from:SINGAPORE,SINGAPORE,SINGAPORE

Quote
The error will become larger and larger if you try to put function in the "Fixed relations" page instead of define a function.

You might need to have some background about numerical method in order to fully understand it.
I will try to explain it starts with how EJS works:
1. EJS will read "Variables" page: it will ask computer to allocate memory for each variables defined and assigning value to the memory. for example x=5.; it will put 0.05 to memory assigned for x.
2. EJS will read "Initialization" page: This is an optional page. EJS will do whatever java command
3. EJS will read "Evolution" page: This is the page to move time step forward and caculate variable value at the next time step. EJS will use user selected "solver" method to solve the ODE.
4. EJS will read "Fixed relation" page: The system has been moved from time t to t+dt. And java code in this page will be executed. But all the command execute here is done after step 3.
5. EJS will show all the GUI elements according to value assigned to element's properties.
6. EJS will go back to step 3 if it in play mode and continue process 3-4-5 cycles, otherwise it will stop.

Good overview of computing in EJS.

Quote
Let me use a simple harmonic motion as examples to explain why relations for derivatives can not be defined in the "Fixed relation page".
If the equation we try to solve is
F=m*d2x/dt2= -k*x-b*vx;
we will change the above second order differential equation into two first order differential equations:
dx/dt=vx;
dvx/dt=-k/m*x-b/m*vx;

Evolution page can do first order differential equation, that's cool.

Quote
You need to select Runge-Kutta (4th) order to reduce the error.
If the Euler's method were selected , the error will accumuate will soon.
The middle method will reduced error to 10% of the Euler's method .
The error for Runge-Kutta (4th) order is 0.00% of the Euler's method.
And this is the best method for most of the cases.

If you change the above equestions to
dx/dt=vx;(STEP 3 )
dvx/dt=ax; (STEP 3)
and define ax=-k/m*x-b/m*vx; at fixed relation page (STEP 4)
The error will be the same as Euler's method even if Runge-Kutta (4th) order was selectd.

Because EJS use the above equation to move time from t to t+dt (STEP3 to STEP4)
But vx(t) is not the same as vx(t+dt), and x(t) is not the same as x(t+dt)

For Euler's method :
x(t+dt)=x(t)+vx*dt;
vx(t+dt)=vx(t)+k/m*x(t)-b/m*x(t);

For midpoint method:
Two time steps(t+dt/2 and t+dt) were calculated to move from t-> t+dt.

For Runge-kutta 4th order method:
Four time steps were calculated to reduced the error to 1/10000 as Euler's method.
The right hand side need to be evaluate 4 times. If the correct formula is only available at STEP4.
The right hand side will always be the same , the Runge-kutta 4th order method is not calculated with correct value. The error move back to the same as Euler's method.

You need to read document about Runge-kutta 4th order method to fully understand how the numerical method works and why it can reduced the error.

So i get the idea to always use Runge-kutta 4th order method if i want to reduce computing error.
So is the draw back, responsiveness and speed of using the applet?
Cool, i will remember to use Runge-kutta 4th order method in the future

Quote
You need to define function or write the equations into the right hand side of the evolution page, if the right hand side contain variables which is also in the left hand side (like x and vx in the above examples.).
Do you mean another tab on the right of evolution page?

Quote
If the right hand side did not contain any variables at the left hand side, like k,b,m are all constant . Then, you can defined those at fixed relation page.

I am a bit lost here.

Thanks for the tips.
But it that enough for Z-D ?
 « Last Edit: February 17, 2009, 04:59:13 pm by lookang » Logged
Fu-Kwun Hwang
Hero Member

Offline

Posts: 3080

 « Embed this message Reply #4 on: February 17, 2009, 08:38:58 pm »

Sorry for the confusion: for dvx/dt=-k/m*x-b/m*vx; at evolution page
The left hand side means dx/dt , the variable at the left hand side means x.
And the variables at right hand side means: x,vx, k,b,m

What I means was that: you can not write the above equation as
dvx/dt=ax; //at evolution page
and define ax=-k/m*x-b/m*vx; //at fixed relation page

The following information are copied from http://calculuslab.deltacollege.edu/ODE/7-C-3/7-C-3-h.html

### Theoretical Introduction

In the last lab you learned to use Heuns's Method to generate a numerical solution to an initial value problem of the form:

y ' = f(x, y)
y(xo) = yo

We discussed the fact that Heun's Method was an improvement over the rather simple Euler Method, and that though it uses Euler's method as a basis, it goes beyond it, attempting to compensate for the Euler Method's failure to take the curvature of the solution curve into account.  Heun's Method is one of the simplest of a class of methods called predictor-corrector algorithms.  In this lab we will address one of the most powerful predictor-corrector algorithms of all—one which is so accurate, that most computer packages designed to find numerical solutions for differential equations will use it by default—the fourth order Runge-Kutta Method.  (For simplicity of language we will refer to the method as simply the Runge-Kutta Method in this lab, but you should be aware that Runge-Kutta methods are actually a general class of algorithms, the fourth order method being the most popular.)

The Runge-Kutta algorithm may be very crudely described as "Heun's Method on steroids."  It takes to extremes the idea of correcting the predicted value of the next solution point in the numerical solution.  (It should be noted here that the actual, formal derivation of the Runge-Kutta Method will not be covered in this course.  The calculations involved are complicated, and rightly belong in a more advanced course in differential equations, or numerical methods.)  Without further ado, using the same notation as in the previous two labs, here is a summary of the method:

xn+1 = xn + h

yn+1 = yn + (1/6)(k1 + 2k2 + 2k3 + k4)

wherek1 = h f(xnyn)
k2 = h f(xn + h/2, yn + k1/2)
k3 = h f(xn + h/2, yn + k2/2)
k4 = h f(xn + hyn + k3)

Even though we do not plan on deriving these formulas formally, it is valuable to understand the geometric reasoning that supports them.  Let's briefly discuss the components in the algorithm above.

First we note that, just as with the previous two methods, the Runge-Kutta method iterates the x-values by simply adding a fixed step-size of h at each iteration.

The y-iteration formula is far more interesting.  It is a weighted average of four values—k1, k2, k3, and k4.  Visualize distributing the factor of 1/6 from the front of the sum.  Doing this we see that k1 and k4 are given a weight of 1/6 in the weighted average, whereas k2 and k3 are weighted 1/3, or twice as heavily as k1 and k4.  (As usual with a weighted average, the sum of the weights 1/6, 1/3, 1/3 and 1/6 is 1.)  So what are these ki values that are being used in the weighted average?

k1 you may recognize, as we've used this same quantity on both of the previous labs.  This quantity, h f(xnyn), is simply Euler's prediction for what we've previously called Δy—the vertical jump from the current point to the next Euler-predicted point along the numerical solution.

k2 we have never seen before.  Notice the x-value at which it is evaluating the function fxn + h/2 lies halfway across the prediction interval.  What about the y-value that is coupled with it?  yn + k1/2 is the current y-value plus half of the Euler-predicted Δy that we just discussed as being the meaning of k1.  So this too is a halfway value, this time vertically halfway up from the current point to the Euler-predicted next point.  To summarize, then, the function f is being evaluated at a point that lies halfway between the current point and the Euler-predicted next point.  Recalling that the function f gives us the slope of the solution curve, we can see that evaluating it at the halfway point just described, i.e. f(xn + h/2, yn + k1/2), gives us an estimate of the slope of the solution curve at this halfway point.  Multiplying this slope by h, just as with the Euler Method before, produces a prediction of the y-jump made by the actual solution across the whole width of the interval, only this time the predicted jump is not based on the slope of the solution at the left end of the interval, but on the estimated slope halfway to the Euler-predicted next point.  Whew!  Maybe that could use a second reading for it to sink in!

k3 has a formula which is quite similar to that of k2, except that where k1 used to be, there is now a k2.  Essentially, the f-value here is yet another estimate of the slope of the solution at the "midpoint" of the prediction interval.  This time, however, the y-value of the midpoint is not based on Euler's prediction, but on the y-jump predicted already with k2.  Once again, this slope-estimate is multiplied by h, giving us yet another estimate of the y-jump made by the actual solution across the whole width of the interval.

k4 evaluates f at xn + h, which is at the extreme right of the prediction interval.  The y-value coupled with this, yn + k3, is an estimate of the y-value at the right end of the interval, based on the y-jump just predicted by k3.  The f-value thus found is once again multiplied by h, just as with the three previous ki, giving us a final estimate of the y-jump made by the actual solution across the whole width of the interval.

In summary, then, each of the ki gives us an estimate of the size of the y-jump made by the actual solution across the whole width of the interval.  The first one uses Euler's Method, the next two use estimates of the slope of the solution at the midpoint, and the last one uses an estimate of the slope at the right end-point.  Each ki uses the earlier ki as a basis for its prediction of the y-jump.

This means that the Runge-Kutta formula for yn+1, namely:

yn+1 = yn + (1/6)(k1 + 2k2 + 2k3 + k4)

is simply the y-value of the current point plus a weighted average of four different y-jump estimates for the interval, with the estimates based on the slope at the midpoint being weighted twice as heavily as the those using the slope at the end-points.

As we have just seen, the Runge-Kutta algorithm is a little hard to follow even when one only considers it from a geometric point of view.  In reality the formula was not originally derived in this fashion, but with a purely analytical approach.  After all, among other things, our geometric "explanation" doesn't even account for the weights that were used.  If you're feeling ambitious, a little research through a decent mathematics library should yield a detailed analysis of the derivation of the method.

 Logged
lookang
Hero Member

Offline

Posts: 1787

http://weelookang.blogspot.com

 « Embed this message Reply #5 on: February 17, 2009, 08:55:42 pm »

Sorry for the confusion: for dvx/dt=-k/m*x-b/m*vx; at evolution page
The left hand side means dx/dt , the variable at the left hand side means x.
And the variables at right hand side means: x,vx, k,b,m

What I means was that: you can not write the above equation as
dvx/dt=ax; //at evolution page
and define ax=-k/m*x-b/m*vx; //at fixed relation page

oic.... i understand this. use dvx/dt=-k/m*x-b/m*vx; at evolution page.

give Z-D time to understand this.

cheers!
 Logged
Z-D
Newbie

Offline

Posts: 3

 « Embed this message Reply #6 on: February 18, 2009, 09:17:25 pm »

WOW!!! Thank you!

About Runge-Kutta, I have studied it and I have done C-Programming on Runge-Kutta as a coursework, so I think I should be able to understand (if I can recall what I've learnt) what's going on.

But the problem here is that the equations of motion which I am trying to work out (Attachment), includes more than one second order derivatives in a single equation. (Btw, I just realised that the equations I used which I uploaded in the first post has some mistakes.)

How do I go about this?... Do I have no choice but to substitute and rearrange the equations of motion so that there will be only one double derivative term in each?

Thanks again!!
 Equation of Motion.jpg (24.86 KB, 725x281 - viewed 369 times.) Logged
Fu-Kwun Hwang
Hero Member

Offline

Posts: 3080

 « Embed this message Reply #7 on: February 18, 2009, 10:48:42 pm »

You can transform the above 4 equations into another set of (4) equations which only contain one second order derivate in each equation.
Treat those four second order derivates as indepentent variables, all the other variables including first order derivates as constant.
You have 4 indepent variables and 4 equations, so you can find solution for each  indepent variables.
(solve linear equations)
 Logged
Z-D
Newbie

Offline

Posts: 3

 « Embed this message Reply #8 on: February 23, 2009, 09:11:54 pm »

^
Thank you very much. I very much appreciate your help.
 Logged
HUMSHAPHY
watchlist

Offline

Posts: -10

 « Embed this message Reply #9 on: November 16, 2009, 10:51:18 pm »

I have a small problem with 'Elecard Converter Studio AVC HD Edition' that I couldn't quite solve myself.
Im trying to Split a .mkv file in two different parts but anytime I specify the part to convert by the gray area which it is supposed to convert only...
it converts straight to the end of the video file without noticing the end part of the encoding, Though it only starts the convert by the specified spot.
How do I make it stop at the specified spot?
 Logged

Электронные каталоги автозапчастей
 Pages: [1]   Go Up
You cannot always have happiness but you can always give happiness. ..."Mother Teresa(1910-1997, Roman Catholic Missionary, 1979 Nobel Peace Prize)"

 Related Topics Subject Started by Replies Views Last post Coupled springs Dynamics Fu-Kwun Hwang 0 24139 June 25, 2005, 10:36:53 am by Fu-Kwun Hwang I need help on java programming specifically on applet devel Information and Download barffour 2 11156 June 06, 2006, 11:19:21 pm by rhipple Coupled springs dynamics ahmedelshfie 0 3312 May 28, 2010, 07:07:09 pm by ahmedelshfie simultaneous combination of lenses and mirrors (with images) Request for physics Simulations sen-h 9 10500 July 24, 2018, 04:04:19 pm by sophiahustler Coupled Pendulum Dynamics Fu-Kwun Hwang 0 6433 April 11, 2013, 09:00:07 am by Fu-Kwun Hwang
Page created in 0.094 seconds with 22 queries.since 2011/06/15