Solve a second order differential equation with ODE45

2 次查看（过去 30 天）
Peter Grant
Peter Grant 2014-1-24
评论： Peter Grant 2014-1-25
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 个评论
Peter Grant
Peter Grant 2014-1-24
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.

请先登录，再进行评论。

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

采纳的回答

Mischa Kim
Mischa Kim 2014-1-25
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 个评论
Peter Grant
Peter Grant 2014-1-25
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.

请先登录，再进行评论。

更多回答（0 个）

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by