# Solving non linear ode using ode45

3 views (last 30 days)

Show older comments

I want to solve a problem of die casting where a liquid alloy contained in a crucible is pressurised by argon into a tube against gravity. The ODE i have derived is P.A- pho.a.y.g= pho.a.y.(y'').

where pho=2375kg/m^3

P=7x10^5 Pa.

A=0.015362 m^2

A/a=581.45

g=9.8m/s^2

y''= 2nd order time derivative of y;

at t=0,y=0,

at t=tf,y=L=1m;

The simulation ends when y=L is reached.

tf=time to fill the tube.

The problem I am facing is during discretization where the variable y is coming in the denominator. so for the first iteration y=0 yeilds an infinitely large no. Please help

##### 1 Comment

Walter Roberson
on 16 Dec 2022

### Answers (2)

Walter Roberson
on 16 Dec 2022

P.A- pho.a.y.g= pho.a.y.(y'')

P and A are nonzero finite constants as are everything else except y and y''

y(0)=0

P*A - pho*a*y(0) = pho*a*y(0)*y''(0)

P*A - 0*constant = constant*0*y''(0)

Assuming y''(0) is finite this gives us

P*A - 0 = 0

and the system is inconsistent unless P or A are 0.

However that logic does not work if in fact y''(0) is infinite, in which case the right hand side has 0*infinity in which case it could potentially be the case that even though you have 0*infinity, that potentially the limit of the product is P*A. But for that to be the case, y(0) = 0 would have to be discontinuous with the rest of y(t), which is a situation that the ode routines cannot deal with.

So... you cannot do that. The equations are inconsistent in a way that cannot be worked with.

##### 0 Comments

Sam Chak
on 16 Dec 2022

Edited: Sam Chak
on 16 Dec 2022

I don't know how you derived this, but after rearranging the equation

you can see the state variable is at the denominator.

Thus, the initial value cannot be absolutely zero, or else it triggers a singularity caused by the Division-by-zero event.

If you select a sufficiently small value for , then you can still run the simulation.

rho = 2375; % liquid aluminum density

P = 7e5; % pressure

A = 0.015362; % area 153.62 cm^2

Aoa = 581.45; % ratio A/a

a = A/Aoa;

g = 9.8;

odefcn = @(t, y) [y(2);

(P*A)/(rho*a*y(1)) - g];

tspan = [0 7.15e-4]; % simulation time

y0 = [0.001 0]; % assume that 1 mm thick layer of residue in the tube

[t, y] = ode45(odefcn, tspan, y0);

tubeLg = y(:,1);

tubeLg(end) % length of tube filled

plot(t, y(:,1), 'linewidth', 1.5), grid on

xlabel('Time, (seconds)')

ylabel('Length of Tube filled, (meters)')

##### 2 Comments

Sam Chak
on 16 Dec 2022

If the velocity should be decreasing, then the dynamics of system (ODE) must be intrinsically stable based on logical thinking.

However, this is not case with the given ODE:

.

Due to the physical quantity (length of tube filled with liquid aluminum) of y is non-negative, , and is positive, then first term is adding energy to the system in a reciprocal manner until the term is dominated by the energy dissipation from the effect of gravity.

Based on the above prediction, if you simulate the system long enough, you should observe an arc like this:

rho = 2375; % liquid aluminum density

P = 7e5; % pressure

A = 0.015362; % area 153.62 cm^2

Aoa = 581.45; % ratio A/a

a = A/Aoa;

g = 9.8;

odefcn = @(t, y) [y(2);

(P*A)/(rho*a*y(1)) - g];

tspan = [0 5.49696e2]; % simulation time

y0 = [0.001 0]; % assume that 1 mm thick layer of residue in the tube

[t, y] = ode45(odefcn, tspan, y0);

plot(t, y(:,1), 'linewidth', 1.5), grid on

TLDR: The model must be inaccurate. Can you recheck the derivation steps?

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!