ODE for system of non linear, second order differential equations

I need to solve a system of 3 equations in the variable x1,x2,x3. Normaly I solve differential equations with ode solvers but in this system I have some problem with non linearity. I do not know how write the ode function that takes into account a term of a second order derivative of x2 in equation 1 or similar situations. I have a system of differential equation like that:

ddx1=(F1(t)-b11*dx1-a13*ddx3-b13*dx3-a15*ddx5-b15*dx5)/(m+a11)      %%Eq. 1
ddx3=(F2(t)-b33*dx3-c33*dx3-a35*ddx5-b35*dx5-c35*dx5-a31*ddx1-b31*dx1)/(m+a33)                                                    %%Eq. 2
ddx5=(F3(t)-b55*dx5-c55*x5-a51*ddx1-b51*dx1-a53*ddx3-b53*dx3-c53*x3)/(I22+a55)      
                                           %%Eq. 3

All the cofficients are known.

I do not know how write in the ode function for this system. Can somebody please explain or write an example of the ode function required to solve a non linear system like that? I would be grateful. I found a similar post where I wrote but I can not get the meaning.

Best regards Alessandro Antonini

Answer 1

[Jonathan Epperl]

I'm not entirely sure I understand your notation, but assuming by dx5 you mean the first derivative of x5, and by ddx5 you mean the 2nd derivative, there are 2 possibilities:

• Either you can rewrite above system so that the left hand sides are ddx1, ddx3, ddx5 and the right hand sides only contain up to first derivatives. If so, then define new variables Y, as Y1=x1, Y2=dx1, Y3=x3, Y4=dx4 and so on. In those new variables you now have a system of 1st order ODEs that you can solve with ode23 or whatever. If you don't understand, http://www.math.uiowa.edu/ftp/atkinson/ENA_Materials/Overheads/sec_8-7.pdf seems to be a decent explanation.
• If you are not able to do the above, then you might be dealing with what is called a descriptor system (there might be other names I am not aware of). That is more complicated, but a quick google search brought up a Matlab toolbox for those kinds of things, maybe that'll help you some more: http://elib.dlr.de/11629/1/varga_cacsd2000p2.pdf

[Alessandro Antonini]

Thanks for the answer Jonathan,

Yes, you are right about the notation. I tried to write a function for this system,doesn't give me any error but I'm not sure about the results, the function would be:

    function dxdt=Newton(t,x)

    global m FeS alphaS b11 a13 b13 a15 b15 a11 FeH alphaH b33 c33 a35 b35 c35 a31 b31 a33...
        MeP alphaP b55 c55 a51 b51 a53 b53 c53 I22 a55 ome 
    dxdt=zeros(6,1);
    dxdt=[x(2);

        ((F1(t))-(b11*x(2))-(b13*x(4))-(a13*(dxdt(4)))-(a15*(dxdt(6)))-(b15*...(x(5)))-(500000*x(1)))/(m+a11);

        x(4);

        ((F2(t))-(b33*x(4))-(c33*(x(3)))-(a35*(dxdt(6)))-(a31*(dxdt(2)))-(b35*(x(6)))-(c35*...(x(5)))-(b31*x(2)))/(m+a33);

    x(6);

        ((F3(t))-(b55*(x(6)))-(c55*(x(5)))-(a51*(dxdt(2)))-(a53*...(dxdt(4)))-(b51*x(2))-(b53*x(4))-(c53*(x(3))))/(I22+a55)];

    return

Could you please have a look on what I did and comment on it? If you think the way it was done is fine? I don't exactly get the first possibility that you described above, can you give me an example for just one variable?It would be really helpfull.

Thanks again Alessandro

[Jonathan Epperl]

So what you did is wrong, you are using elements of dxdt when you are computing dxdt, e.g. in the 2nd row you use dxdt(4) and dxdt(6). At this time though, dxdt is still a vector of all zeros, so dxdt(4) and so on all evaluate to zero.

Here is a bad example (bad because it is not an example of anything physical):

(1)

dx1  = -x1  + ddx2
ddx2 = ddx1 - x2

Now try and rewrite that so that the 2nd derivatives are on the left, and on the left only. Check it, but the above should be the same as

(2)

ddx1 = x1 + x2 + dx1
ddx2 = x1 + dx1

Now introduce a new variable, say y as y = [x1; dx1; x2; dx2] and you can write the original system of odes as a 1st order system:

(3)

dy = [ y(2);
       y(1) + y(3) + y(2);
       y(4);
       y(1) + y(2)];

In more complicated cases, the step from (1) to (2) should be done by rewriting the system as one equation, using matrices, (1) can be written as

[1 0;0 1]*x + [1 0; 0 0]*dx + [0 -1; -1 1]*ddx = 0

Then you invert the matrix in front of ddx and you arrive at (2). If in your case the matrix in front of the 2nd derivatives is not invertible, then you're in trouble and will have to look at the descriptor systems stuff.

Lastly: I don't think your system is nonlinear, for that to be the case it would have to involve nonlinear functions of x, dx, ddx and so on.

[Alessandro Antonini]

Thank you very much, I've really appreciated!