MATLAB: NaN values when using trapz and second order coupled ode

modified eulerode45trapz

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.
%

Best Answer

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.
Related Question