How to implement tightly coupled nonlinear odes using ode45 in matlab?

I am solving a problem from fluid dynamics; in particular tightly coupled nonlinear ordinary differential equations. The following is a scaled-down version of my actual problem. I have solved system of coupled odes many times in the past but this case is different since double derivatives of one variable depends on the double derivative of another variable. How do I implement it in ode45? I need 3 x 2 = 6 plots of x, x-dot and x-ddot versus time for t, 0 to 2. All required initial conditions have zero values at t = 0 How do I store the updated value of the double derivatives as the ode45 code runs? The way ode45 works, I get x and x-dot as output but not the double derivatives. Any help will be highly appreciated.

[John D'Errico]

Did you read the documentation for ODE45? In there as I recall, it shows you how to convert higher order differentials into a system of first order problems.

[Vikash Pandey]

@ John, I know how to do that. I have solved system of odes multiple times before. But for me the problem is how to take control of the updated values of the double derivatives as the code runs. I can put up the code within 15 minutes. Wait.

[madhan ravi]

what are the intial conditions??

[Vikash Pandey]

Please see below and confirm if I am correct or wrong.

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

>

function ydot = f(t, y);

ydot(1) = y(2);
ydot(2) = t*y(4) - y(1)*y(5) - y(6);
ydot(3) = y(3);

ydot(4) = y(5);
ydot(5) = (t^2)*y(1) + y(4)*y(2) - y(3);
ydot(6) = y(6);

ydot = ydot';

------

clc;
clear all;
close all;
time_range = [0 2];
initial_conditions = [2 1 0.5 4 3 1]; % x1(0) y1(0) y1-dot(0) x2(0) y2(0) y2-dot(0)
[t, y] = ode45('bubble', time_range, initial_conditions);
x1=y(:,1);
x1dot=y(:,2);
x1ddot = y(:,3);
x2=y(:,4);
x2dot=y(:,5);
x2ddot = y(:,6);
figure(201);
h1 = plot(t, x1, 'b-', t, x2, 'k:');
set(h1,'linewidth',4);
legend(['x1'], ['x2'])
legend('Location','Northwest')
set(gca,'xscale','linear','FontSize',26)
set(gca,'yscale','linear','FontSize',26)
set(gca,'XMinorTick','on')
set(gca,'YMinorTick','on')

grid off
xlabel('Time, $t\mbox{ } (s)$','interpreter','latex','FontSize',34);
ylabel('Distance, $x_{1,2 \mbox{ }} (m)$','interpreter','latex','FontSize',34)
hold off
figure(202);
h2 = plot(t, x1dot, 'b-', t, x2dot, 'k:');
set(h2,'linewidth',4);
legend(['x1-dot'], ['x2-dot'])
legend('Location','Northwest')
set(gca,'xscale','linear','FontSize',26)
set(gca,'yscale','linear','FontSize',26)
set(gca,'XMinorTick','on')
set(gca,'YMinorTick','on')
grid off
xlabel('Time, $t\mbox{ } (s)$','interpreter','latex','FontSize',34);
ylabel('Speed, $\dot{x}_{1,2} \mbox{ } (ms^{-1})$','interpreter','latex','FontSize',34)
hold off

figure(203);
h3 = plot(t, x1ddot, 'b-', t, x2ddot, 'k:');
set(h3,'linewidth',4);
legend(['x1-ddot'], ['x2-ddot'])
legend('Location','Northwest')
set(gca,'xscale','linear','FontSize',26)
set(gca,'yscale','linear','FontSize',26)
set(gca,'XMinorTick','on')
set(gca,'YMinorTick','on')
grid off
xlabel('Time, $t\mbox{ } (s)$','interpreter','latex','FontSize',34);
ylabel('Acceleration, $\ddot{x}_{1,2} \mbox{ } (ms^{-2})$','interpreter','latex','FontSize',34)
hold off

-------

And my plots are:

<<

<<

>>

>>

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.

[madhan ravi]

I have always used FORTRAN to solve such equations,

so does the result coincide with fortran?

[Vikash Pandey]

@ Madhan, I have not gone to FORTRAN for this problem. That is a bit more time consuming there. And the actual problem is much more complex. 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.

[Torsten]

 ydot(3)=y(3) and ydot(6)=y(6) ? Looks wrong to me.

[Vikash Pandey]

@ Torsten, you may be correct. How do I fix that? 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 x2 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.

I cannot post more images, so I am putting up the links below.

https://postimg.cc/xJJknsmz

https://postimg.cc/DJ9SZqQY

Answer 1

[Torsten]

Solve the equations

t*x2 - x1*x2' = t^2*x1 + x2*x1'
x1'' + x2'' = t*x2 - x1*x2'

Setting

y1 = x1
y2 = x1'
y3 = x2
y4 = x2'

you arrive at the system

y1' = y2
y2' + y4' = t*y3 - y1*y4
y3' = y4
y1'*y3 + y1*y3' = t*y3 - t^2*y1

Now either solve for y1', y2, y3' and y4' in each time step or use the mass matrix facility of the ODE solvers by writing your system as

M(t,y)*y' = F(t,y)

with

M = [1 0 0 0; 0 1 0 1; 0 0 1 0; y3 0 y1 0]
F = [y2;t*y3-y1*y4;y4;t*y3-t^2*y1]

Best wishes

Torsten.

Answer 2

[Vikash Pandey]

@ Torsent, Thanks. 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 this: Modified problem

[Vikash Pandey]

Even in this case my double derivative plots has zero values all the way if the ICs for the double derivatives are zero. 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?

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

[Torsten]

You don't need initial conditions for the second derivatives.

You have four equations, and initial values have to be prescribed for x1, x1', x2 and x2'.

[Torsten]

If your equations are like in the modified problem, set

y1=x1
y2=x1'
y3=x1''
y4=x2
y5=x2'
y6=x2''

and solve the DAE system

y1'=y2
y2'=y3
y4'=y5
y5'=y6
y3+t*y6^2-t*y4+y1*y5=0
t*y3+y6-t^2*y1-y4*y2=0

As for the former case, initial conditions are only necessary for x1, x1', x2 and x2'.

[Vikash Pandey]

I get it Torsten. So essentially we transforming this ode problem to system of equations, correct? I just got another idea: I should substitute the value of x1 doubledot in the equation for x2-double dot? And so the expression for x1 doubledot will only have terms x1, x2, x1doubledot(obvious), x1dot, x2dot but NOT x2doubledot? Same goes for the expression for x2 doubledot which will have all terms but NOT x1doubledot? And then implement ode45? What do you think?

[Torsten]

I don't understand how you want to formulate the problem for ODE45 in this case.

Usually, you can explicitly solve for x1'' and x2'' as expressions of x1,x1',x2,x2' and t (two equations in the two unknowns x1'' and x2''). This is in principle possible for your modified problem (with ugly right-hand side F) , but not for your original problem (why ?).

If it's possible, you can write your system as a first-order system consisting of 4 equations.

If it's not possible, you can write your system as a system of 4 ODEs and 2 AEs as I did for the modified problem.

[Vikash Pandey]

The actual problem that I wish to solve is this Real problem The equation is for R2 but it depends on R1'' (look at the last term). Once can obtain equation for R1 by interchanging the indices 1 and 2. And most papers have solved this or variants of this problem using ode45. What do you say?

[Vikash Pandey]

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.

[Torsten]

I see no difference to the modified problem from above. So you can solve it the way I suggested.

[Vikash Pandey]

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

[Huy Nguyen]

hi, I have the same problem. can you help me or share the code please?

[Sam Chak]

Hi @Huy Nguyen

Would advise you to post your specific problem in a New Question. Also to clarify whether it's a math-related problem or a technical problem in MATLAB code.

[Huy Nguyen]

hi, thanks for your reply. I fix the problem.