NaN values when using trapz and second order coupled ode

Question

Meva 2016-9-16

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/303285-nan-values-when-using-trapz-and-second-order-coupled-ode

编辑： Meva 2016-9-20

采纳的回答： Walter Roberson

在 MATLAB Online 中打开

Hi, Firstly, apologise for a long question. I could not make it shorten. Let me explain the problem.

The image above shows that I have two copupled nonlinear equations for h and theta. These h and theta are unknown functions of time T. The target is to find these two functions h and theta. Initial conditions are given in the code.

I used trapz to compute integrations. Then I used modified Euler method, Runga-Kutta 4th order and Matlab ode45 to solve this system and compare solutions. When I debug my code, integral_func has NaN values then it takes numerical values. It's size is 1*101. The code is below:

function euler  
    % A second order differential equation
    Tsim = 1;                                            % simulate over Tsim seconds
    dt = 0.01;                                           % step size
    N= Tsim/dt;                                          % number of steps to simulate
    x=zeros(4,N);
    t = zeros(1,N);
%----------------------------------------------------------------
b = 0.4;                                                       %|
kappa = 1;                                                     %|
M = .1;                                                        %|
g = .1;                                                        %|
Gamma = 0.2;                                                   %|                                                             
h0 = kappa*sin(pi*b);        % h0 = F(b);                      %|
theta0 = kappa*pi*cos(pi*b);                                      
%-------------------------------------------------------
x(1,1)= h0;        % h(t = 0) ;
x(2,1)= 0;         % dh/dt (t = 0) = h1; 
x(3,1)= theta0;    % theta(t = 0);
x(4,1)= 0;         % dtheta/dt(t = 0) = theta1; 
for k=1:N-1
    t(k+1)= t(k) + dt;
    xnew=x(:,k)+ dt *MassBalance(t(k),x(:,k));
    x(:,k+1) = x(:,k) + dt/2* (MassBalance(t(k),x(:,k)) + MassBalance(t(k),xnew)); %Modified Euler algorithm
end
%Compare with ode45
x0 = [h0; 0; theta0; 0];
[T,Y] = ode45(@MassBalance,[0 Tsim],x0);   % where Y is the    solution vector      
end
function dx= MassBalance(t,w)
  %----------------------------------------------------------------
  b = 0.4;                                                       %|
  kappa = 1;                                                     %|
  M = .1;                                                        %|
  g = .1;                                                        %|
  Gamma = 0.2;                                                   %|                                                             
  h0 = kappa*sin(pi*b);        % h0 = F(b);                      %|
  theta0 = kappa*pi*cos(pi*b);                                      
  %--------------------------------------------------------------
  % Equations of motion for the body
  dx = zeros(4,1);
  xx = (0:0.01:1);
  integral_func = 1/2*(1 - ((w(1)+ (1-b)*w(3))./(w(1)+ (xx-b).*w(3)-kappa.*sin(pi.*xx))).^2)
dx(1)=w(2);
dx(2)=trapz(xx, integral_func) - M*g
dx(3)=w(4);
dx(4)= 1/Gamma.* trapz(xx, (xx-b).*integral_func); 
end

When I plot h and theta versus time, I got different results from these three solvers: modified Euler, Runga-Kutta, ode45. I think I have more than one mistake. Please, any help will be appreciated ..

Thanks.

%

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

请先登录，再进行评论。

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

Answer 1

Walter Roberson 2016-9-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/303285-nan-values-when-using-trapz-and-second-order-coupled-ode#answer_235131

在 MATLAB Online 中打开

In your line

integral_func = 1/2*(1 - ((x+ (1-b)*y)./(x+(xx-b).*y-kappa.*sin(pi.*xx))).^2)

x and y are undefined. We might suspect that x should be replaced by xx, but there is no obvious substitution for y.

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

Walter Roberson 2016-9-20

It is not generally a problem to have an unresolved int() inside a dsolve(), or at least it is not an error (but possibly the solver is not powerful enough to deal with the situation.)

Where I wrote "boundary condition", that is the same thing as an initial condition at the initial time: your system looks to me like it needs more initial conditions. However I need to think about this more. (I am not sure I will come up with anything.)

Meva 2016-9-20