# NaN values when using trapz and second order coupled ode

21 views (last 30 days)
Meva on 16 Sep 2016
Edited: Meva on 20 Sep 2016
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. %

Walter Roberson on 18 Sep 2016
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.

Meva on 20 Sep 2016
Yes Walter, In the previous comment, I meant that in order to use dsolve(), I need to use int() as the integration should be symbolic so that the inputs of dsolve() should be symbolic as well. My physical system uses just initial conditions. Boundary conditions are not known in advance. So, could you please inform me that is there any alternative way to solve this initial value problem?
Thanks so much.
Walter Roberson on 20 Sep 2016
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 on 20 Sep 2016
Thanks so much Walter. Looking forward to hearing from you. I would be very grateful if you come across anything.