NaN values when using trapz and second order coupled ode

2 views (last 30 days)
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
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
%Compare with ode45
x0 = [h0; 0; theta0; 0];
[T,Y] = ode45(@MassBalance,[0 Tsim],x0); % where Y is the solution vector
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(2)=trapz(xx, integral_func) - M*g
dx(4)= 1/Gamma.* trapz(xx, (xx-b).*integral_func);
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 ..

Accepted Answer

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

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!