Solve a second order differential equation with ODE45

1 vue (au cours des 30 derniers jours)
Peter Grant
Peter Grant le 24 Jan 2014
Commenté : Peter Grant le 25 Jan 2014
Hello,
I had managed to make this file work but recently it's very slow and I'm not able to run it properly. I think the issue might come form the vectors. Here is the script and function file:
function xdot = OWC (t,x)
global g h a
xdot = [((g*h-((x(2).^2)/2)-g*x(1))/(x(1)+a)); x(1)];
and the script:
global g h a;
g=9.81;
h=0;
for a=0;
[t,x] = ode45('OWC',[0 10],[1.0 1.0]);
plot (t,x)
hold on
end
for a=1;
[t,x] = ode45('OWC',[0 10],[1.0 1.0]);
plot (t,x)
hold on
end
The original equation is: y"=(g*h-(y'^2/2)-g*z)/(z+a)
I haven't used Matlab much so I wouldn't be surprised if I'm using the wrong technique but I did manage to make this script work before and somehow I've ruined it.
Thank you.
  1 commentaire
Peter Grant
Peter Grant le 24 Jan 2014
Whenever I change the y0 in
[t,x] = ode45('OWC',[0 10],[1.0 1.0]);
From [1.0 1.0] to [0 0], the script finishes but without giving an answer. That's why I think the slow process of the script might be because of these initial conditions.

Connectez-vous pour commenter.

Réponse acceptée

Mischa Kim
Mischa Kim le 25 Jan 2014
Hello Peter, there are a couple of issues:
  • It looks like you want to use if rather than for statements to execute either of the blocks depending on the value of a , correct?
  • When you set the initial conditions to [0; 0] and with h = 0 the derivatives are all zero at all time so you get a flat line (provided that a != 0).
  • When you set the initial conditions to [0; 0] and with a = 0, h = 0 you have a devide by zero the differential equations resulting in NaN.
  1 commentaire
Peter Grant
Peter Grant le 25 Jan 2014
Thank you! This really helped, I had realized that by changing 'h' my script seemed to run better but I didn't really try to figure out why. It was stupid of me not to look into the mathematical aspect of the script rather than technical aspect.

Connectez-vous pour commenter.

Plus de réponses (0)

Community Treasure Hunt

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

Start Hunting!

Translated by