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

4 次查看（过去 30 天）

Vikash Pandey 2018-11-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/429363-how-to-implement-tightly-coupled-nonlinear-odes-using-ode45-in-matlab

评论： Huy Nguyen 2022-12-7

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.

8 个评论
显示 6更早的评论隐藏 6更早的评论

Vikash Pandey 2018-11-12

编辑：Vikash Pandey 2018-11-12

在 MATLAB Online 中打开

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.

Torsten 2018-11-12

在 MATLAB Online 中打开

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

Vikash Pandey 2018-11-12

编辑：Vikash Pandey 2018-11-12

@ 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

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten 2018-11-12

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/429363-how-to-implement-tightly-coupled-nonlinear-odes-using-ode45-in-matlab#answer_346494

在 MATLAB Online 中打开

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

类别

MATLAB Programming

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by

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

8 个评论
显示 6更早的评论隐藏 6更早的评论

采纳的回答

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

更多回答（1 个）

12 个评论
显示 10更早的评论隐藏 10更早的评论

另请参阅

类别

标签

Community Treasure Hunt

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

8 个评论 显示 6更早的评论隐藏 6更早的评论

采纳的回答

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

更多回答（1 个）

12 个评论 显示 10更早的评论隐藏 10更早的评论

另请参阅

类别

标签

Community Treasure Hunt

8 个评论
显示 6更早的评论隐藏 6更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

12 个评论
显示 10更早的评论隐藏 10更早的评论