NaN values when using trapz and second order coupled ode

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.

%

Answer 1

[Walter Roberson]

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.

[Meva]

Thanks so much for your answer Walter. Yes, you are right, it was a copy paste mistake. I have edited the question. w(1) and w(3) are solution vector components . w(1)is denoted as h, w(3) is denoted as theta.

[Walter Roberson]

In the situation that w(1) == kappa*sin(pi*b) then you will integral_func will have a denominator of 0 at the point that x == b . Your b happens to be a subset of your xx points so this will happen if w(1) ever happens to equal kappa*sin(pi*b) . And for whatever reason, it does come to exactly equal it.

Your integral_func has this potential singularity at x == b, going to negative infinity when it happens. I suspect that the singularity acts to "tune" the values until they match, creating a "well of potential"

[Meva]

Yes, you are absolutely right. There is a singularity. I see that the integrand has a singularity at x==b, I have decreased the increment in this line: xx = (0:0.003:1); which was before xx = (0:0.01:1). So, the when the index 41, my xx was equal to b which gives the -infinity. Now, It jumps that b.Thanks so much Walter..

[Meva]

Thanks @Walter again. Now, the problem appears as why I have a huge difference in outputs from modified Euler, ode45 and Runga-Kutta for h and theta that are the solution vector components. I have just inserted an image in the question as I cannot insert here. However, if it is not suitable for this context, I can post another question.

[Walter Roberson]

The different algorithms make predictions of additional points in different ways, that are going to end up including lesser or greater parts of that singularity. The singularity has a large but narrow negative contribution; if you end up skipping over it due to your predictions being different then you might not get that negative contribution so your output might be much higher.

[Meva]

Oww I see, thanks @Walter, do you mind informing me that what do you suggest me to make this integration possible then ? (without skipping over x==b)

[Walter Roberson]

Did you try running it through the Symbolic Toolbox dsolve() function?

[Meva]

Thanks @Walter. However, I suspected that the system cannot be analytically solvable because of the integrations in the RHS. I used

int

before numerically solving the integrations. However, it seems to me it is not analytically possible (as far as I know) since the command window returned

int

again (did not give me the analytical answer of the integration).

[Walter Roberson]

You should be using dsolve() not just int()

However, dsolve() will probably not be able to solve it.

The symbolic toolbox is able to do numeric solutions of ode, but in order to do that you would need to define boundary conditions for h'' and theta'' . You need a boundary condition for each level of differentiation, for each function, so with two functions and three levels of differentiation (not differentiated, prime, double-prime) you are going to need 6 boundary conditions.

[Meva]

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]

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]

Thanks so much Walter. Looking forward to hearing from you. I would be very grateful if you come across anything.