How to solve the second-oder differential equations with two variables using ODE45 or whatever

Hello,

I'm trying to solve the following equations using ODE45 but I'm not confident that my codes work correct. So please teach me whether the codes are correct or not, and if not, teach me how I can analyze the equation numerically.

----—the equation———

mx’’+Dx’+Kx=-f

Tf’ + f = Gx

%parameters

m=3

D=2

K=1

T=1

G=5

%variables

x, f

———my codes———

% I know ode45 only can solve one oder differential equations so I did a change of variables. However I would rather avoid using this change because it’s hard to see x and f in a plotted graph with the new variable appearing on it.

m=3

D=2

K=1

T=1

G=5

% I did the following change [X=x’] so the equation becomes[ X=x’ ][mX’+DX+Kx=-f][Tf’ + f = Gx]

%F(1) matches x, F(2) matches X, F(3) matches f

% x(0)=0.5

[t,F]=ode45(@(t,F)[F(2);(-F(3)-D*F(2)-K*F(1))/m;(G*F(1)+F(3))/T],[0,1],[0.5,0,0]);

plot(t,F);

Answer 1

[Bjorn Gustavsson]

To the best of my "reading-ability" you seem to have it right in ODE-defining function. From your equations I expect them to be something like this:

% mx’’+Dx’+Kx=-f
% Tf' + f = Gx
% From which we define our 3 first order ODEs as:
% dxdt = v;
% dvdt = -D/m*v-K/m*x-f/m
% dfdt = -f/T + G/T*x
% Which we can put into a dynamical function:
ODExvf = @(t,y,m,D,K,T,G) [y(2);
                   -D/m*y(2)-K/m*y(1)-y(3)/m;
                   -y(3)/T + G/T*y(1)];
m=3;
D=2;
K=1;
T=1;
G=5;
[t,F]=ode45(@(t,y) ODExvf(t,y,m,D,K,T,G),[0,1],[0.5,0,0]);
plot(t,F);

It might be preferable to define the ODE more "easily human-readable", and sometimes it's neat to have the opportunity to modify the parameters without having to redefine the function.

HTH

[嵩人 藤森]

Thanks for that!! It help me a lot.

The codes you posted result in the same graph that my codes answered so it seems that my codes are correct after all

I had two aim about this post:

1, To ensure my codes work correct

2, To avoid using the change of variable because it’s hard to see x and f in a plotted graph with the new variable appearing on it

aim.1 seems achived, thank you.

How about aim.2? Is there any way to avoid using the change of variable?

[Bjorn Gustavsson]

If you want to have all 3 unambiguously in one graph you can simply use the legend function something like this:

lh = legend(ph,'x','$\dot{x}$','f');
set(lh,'interpreter','latex','fontsize',12)

That should clarify what is what. Otherwise you'll just have to select which components of y to plot and keep track of what they corresponds to - for this clear code is important, preferably documented. (Many have spent years tracing what they've done and where they stored those results...)

HTH

[嵩人 藤森]

Thanks for that but I mean I wanna see only x and f in plotted graph because when X’ vibrates and diverges to infinity, the plotted line hides those of x and f. I just tried the following codes

plot(t, F(1)) plot(t,F(3))

But it didn’t work :(

How can plot x and f alone?

[Bjorn Gustavsson]

Oh, that's simple enough:

subplot(2,1,1)
ph = plot(t,F);         
lh = legend(ph,'x','$\dot{x}$','f');
set(lh,'interpreter','latex','fontsize',12)
subplot(2,1,2)                             
ph = plot(t,F(:,[1 3]));
legend(ph,'x','f','interpreter','latex')

You'll find all the information about the format of the output from ode45 in the help.

HTH

[嵩人 藤森]

It worked!! I'm so thankful for you.

[Bjorn Gustavsson]

Glad to have helped.