21 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

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

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

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.