NaN values when using trapz and second order coupled ode

6 次查看（过去 30 天）
Meva
Meva 2016-9-16
编辑： Meva 2016-9-20
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
Walter Roberson 2016-9-18
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更早的评论
Walter Roberson
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
Meva 2016-9-20
Thanks so much Walter. Looking forward to hearing from you. I would be very grateful if you come across anything.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by