Hi, if you inject the value of x2'' from the second equation in the first one, the first equation becomes :

x1x2' + x2x1'-tx2-t^2x1 = 0

Thus it seems that your system reduces to only one ODE of the first order in which appears x1' and x'2.

Perhaps a second equation is missing that would lead to a system of 2 ODE possibly solved with ODE45?

Hi Vikash.

Further to Richard's reply, if you gather together all the second derivative terms onto the left hand sides of these equations, then those two left hand sides will be identical. Therefore you cannot isolate $x_1''$ from $x_2''$. If any of your second derivatives were to have a different coefficient then it would be possible to rewrite your equations in the form x1'' = f(x1,x2,x1',x2'), x2"=g(x1,x2,x1,x2').

Please see below and confirm if I am correct or wrong. Please check my approach and my code properly. I have always used FORTRAN to solve such equations, trying Matlab first time for such a problem.

PLEASE check the expression for y-dots carefully. Is my understanding correct?

I am just checking my understanding of ode45 in this scaled-down version of the actual problem. I wish to make sure, I understand ode45 correctly before I go to the actual problem.

I am not doubting your coding of the numerical method (and Matlab should give exactly the same as Fortran with identical methods, as will C++ etc.), but rather I am dubious of your equations. They are inconsistent unless Richard's expression is satisfied. If one of the equations had x1" + c*x2" on the left hand side then x1" and x2" could be decoupled. What is the application, and how were these derived?

What worries me is that if I change the initial conditions to

initial_conditions = [2 1 0 4 3 0]; which means acceleration is zero for both x1 and z2 at t= 0. Then I get the following plots. x1,2 and x1,2-dot plots change, but x1,2-ddot plots do not show any sign of change at all. How is that possible? Velocity is changing but not the acceleration? There must be some bug in the code.

Hi, according to our remarks (Andrew's one and mine) your problem formulation is not correct. In MATLAB, ode45 solves an ODE system of the form y'=f(t,y) so, in the right and side only y must appear and not its derivatives. Or, if on the right hand side a variable y appears, an additional ODE given y' must be defined. Best

The actual problem that I am trying to solve is attached. Note that the equation is for R2 and double derivative of R2 depends on the double derivative of R1. One can get equation for R1 by interchanging the indices 1 and 2.

And I am just trying to gain insight on ode45 on how to take care of the updated values of the double derivatives as code proceeds in its loop.

@ Andrew and Epenoy

I get your comment of the dubiousness of the equations. May be I should not use the term distance, speed and accelerations. Lets regard the variables from a pure mathematical setting; i.e., dimensionless. So in that case, how is that possible that the variable and its first derivative is changing but not its second derivative.

I assume that the equation for R_1 has a T_2 defined in the way that you said!

Now the coefficients of R_1" and R_2" are different from one another, so I think that it is possible to solve these equations. That said, I don't that you can do it immediately with the equations in their current form. I think that you will still need to decouple the equations and write down an equation for R_1 which doesn't contain R_2" and vice versa. You have

aR1" + bR2" = function of R1, R2, R1', R2'

cR1" + dR2" = function of R1, R2, R1', R2'

where a, b, c, and d are functions of first-derivative terms at worst. Therefore R2" can be removed from the first equation and R1" from the second by appropriate multiplications and additions. This will result in a somewhat untidy pair of equations, but all the computer will see is some RHSs which are longer than usual.

May be I should have formulated the problem little different as in this case, the double derivatives can be brought into one side. What if the equation was like attached.

Even in this case my double derivative plots has zero values all the way if the ICs for the double derivatives are zero.

But the situation changes if I have non-zero values of double derivatives as the ICs.

Perhaps this is important, but the equations for R1 and R2 which you quote are very similar indeed to the ones which you have highlighted in a PhysReE article in a different thread. The double-derivative term on the present right hand side is only a single-derivative term in the paper by Mettin et al., or did Mettin et al write down the derivative of R1*R1*R1' incorrectly?

@ Andrew. yes you are right. Mettin had a typo in his equations. I have developed a theory on the top of that. but before I solve my equations I wish to reproduce Mettin's result and so this mess is all about.

So, Andrew are you suggesting that I should substitute the value of R1 doubledot in the equation for R2-double dot? And so the expression for R1 doubledot will only have terms R1, R2, R1doubledot(obvious), R1dot, R2dot but NOT R2doubledot? Same goes for the expression for R2 doubledot which will have all terms but NOT R1doubledot?

Yes, but given that your equations are of the form,

aR1" + bR2" = function of R1, R2, R1', R2'

cR1" + dR2" = function of R1, R2, R1', R2'

one can solve these simultaneous equations for R1" and R2". It'll be a mess, but it is doable.

I do not wanna take approach. Do you think ode45 should be an okay approach to solve following what I mentioned above about the substitution?

ode45 can't be implemented, nor any other IVP solver. Somewhere in the implementation it will be necessary to decouple the equations. If you had the system,

R1" = f(R1,R2,R1',R2')

R2" = g(R1,R2,R1',R2') + something*R1"

then you could find R1" first and then use that to update R2". But in your case R1" and R2" have to be obtained simultaneously. What I suggested earlier works nicely: let the equations be

aR1" + bR2" = F(R1, R2, R1', R2')

cR1" + dR2" = G(R1, R2, R1', R2')

then

R1" = (dF-bG)/(ad-bc)

R2" = (-cF+aG)/(ad-bc)

and hence all these RHS terms involve only R1 and R2 and their first derivatives.

I can always start with R1 and R2 at rest t1. At t2, I first calculate, R1(t2) that depends on the value of R1(t1), R1'(t1), R2(t1), and R2'(t1). Once R1(t2) is obtained, I can obtained R1' and R1'' in a similar manner. The obtained values of R1 and its derivatives at t2 are then used to calculate R2 that will depend on R1(t2), R1'(t2), R2(t1), and R2'(t1). And the process repeats. Do you agree?

Also, I do not think, R2" = g(R1,R2,R1',R2') + something*R1" will be the case. I guess it will be just:

R2" = g(R1,R2,R1',R2') after the substitutions.

yes i agree. Decoupling will be done by substitution of the double derivatives in the "other" equation. By the way, most papers have solved this or variants of this problem using ode45. And then I can use just four standard ICs against the present six ICs issues that I have been facing.

Yes, it's a fourth order system and so four ICs are all that is needed.

A standard initial value problem solver will require ICs for R1, R2, R1' and R2' at t=t1. For each of these four variables, their value at the next timestep, t=t1+dt, will depend on the values of their slopes at t=t1, i.e. on the values of R1', R2', R1" and R2", respectively, and the questions is then: from where does one get the values of R1" and R2".

You mentioned that others have used ode45, but do you know if the authors who have done so have given details about how they did it? It is sometimes the case authors wish to hide certain little secrets and don't mention the one thing that is the most difficult about performing the simulation.

Personally, I would do the simultaneous equations thing that I mentioned and not worry about it (until ad-bc becomes close to zero!). The coding doesn't need to be messy even if the full expressions for R1"=... and R2"=... are very long. I would write something like,

a=...

b=...

c=...

d=...

use these in my Fortran subroutine which defines the right hand sides of the four equations.

@ Andrew Thanks for the input. Well, all authors just mentioned the Rk family they used for solving; Dormand-Prince, Burlisch-Stoer etc.

I do not think that the values of double derivatives is required at all. The expressions finally reduce to the variables and their first derivatives after the substitution. I am trying. I will get back.

Ah, if we have the single second order equation, y"=f(t,y), and we reduce it to first order form using v1=y, v2=y', then the equation becomes,

v1' = v2

v2'= f(t,v1).

For the sake of argument let us use Euler's method, and designate v1_n to be the numerical value of v1 at t=t_n, then Euler stepping for this problem is,

v1_{n+1} = v1_n + dt*v2_n

v2_{n+1} = v2_n + dt*f(t_n,v1_n)

where dt is the timestep. The quantities multiplying dt are the values of the derivatives of v1 and v2 at t=t_n, i.e. the values of y' and y''. So yes, one needs to know the values of the second derivatives. Most people will interpret these solely as the right hand sides of the equations for v1 and v2, but they are the derivatives of the v-variables, and one of them is, in this context, y''. Maybe it is a matter of interpretation rather than of mathematics!

I have solved the problem. Thank you all so much for helping me fix this problem. :)

I decoupled the ODEs before implementing ODE45. It was a hell of a task. But, now I can reproduce the results from published paper quite well, though very minor difference is still there. Little relative shift in the position of the curve is there, otherwise they are almost identical. I guess that could be due to the relative normalization of the axes, or may be the ode solver ode45 has evolved in last 20 years. :)

Hi Vikash.

Congratulations on getting the code to work. Sounds as if it was quite tricky!

Re the very slight differences, they could possibly be due to the number of significant figures used. Or it could be that the maximal Lyapunov exponent could occasionally be positive - this could be tested by varying very slightly the initial conditions which you used.

Interesting

A very relevant discussion!

@ All Finally Solved and published. https://journals.aps.org/pre/abstract/10.1103/PhysRevE.99.042209

Dear Vikash

Can I have your final matlab code?

Regards

Reza