How to solve a highly non-linear ODE containing multiple powers?

I'm working on a time-dependent equation from fluid mechanics for the velocity (u) in a loop. The equation I'm trying to solve is:

Where c1, c2, c3, and c4 are all positive, known constants. I've tried pretty much all of the ODE solvers in MATLAB and they all blow up to infinity and return the following error: "Warning: Failure at t=3.237651e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.150245e-14) at time t."

Can anyone give me some insights into other solvers that might work for this type of problem, or modifications to the solvers I've tried that can resolve the integration tolerance issue?

Answer 1

[Star Strider]

It probably depends on what the c values are. When I run it, it usually takes off to occasionally at some point, however in others, it just becomes asymptotic —

c = rand(4,1)

c = 4x1

0.9593 0.0150 0.0763 0.6945

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

DE = @(t,u) [-u^1.825 u^2 u 1] * c

DE = function_handle with value:
    @(t,u)[-u^1.825,u^2,u,1]*c

[t,u] = ode23s(DE, [0 10], 0);

figure
semilogy(t, u)
grid

.

[Star Strider]

Using your constants —

c = [0.1982; 1.8115; 1.3373; 0.0162]

c = 4x1

0.1982 1.8115 1.3373 0.0162

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

DE = @(t,u) [-u^1.825 u^2 u 1] * c

DE = function_handle with value:
    @(t,u)[-u^1.825,u^2,u,1]*c

[t,u] = ode23s(DE, [0 10], 0);

Warning: Failure at t=3.230728e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.147785e-14) at time t.

figure
semilogy(t, u)
grid

An analytic solution is not an option for a nonlinear differential equation. You can onbly integrate it up to the poiint that it saturates, so anything less than that ‘t’ value(about 3.23) would produce valid results.

.

[Star Strider]

As always, my pleasure!

Answer 2

[Sam Chak]

Hi @Sean

You mentioned that c1, c2, c3, and c4 are all positive known constants. Now, let's consider a special case of your differential equation where c1, c2, and c4 are set to zero. In this simplified scenario, the solution of the differential equation becomes an exponential function of time. However, even if c3 is as small as 10, the response will still exponentially grow and reach a value of 10^13 within 3 seconds.

To better understand the situation, could you please provide the actual values of c1, c2, c3, and c4? This information will help in accurately assessing the behavior of the system.

syms u(t) c3
ode     = diff(u,t) == c3*u;
init    =      u(0) == 1;       % non-zero initial value
uSol    = dsolve(ode, init)

uSol = 

c3val   = 10;       % some positive value
u       = subs(uSol, c3, c3val)

u = 

fplot(u, [0 3]), grid on, xlabel('t'), ylabel('u(t)')

[Sean]

@Star Strider, @Sam Chak

Thanks for the insights! I should have assumed the constants would play a role. The values I currently have are:

c1 = 0.1982

c2 = 1.8115

c3 = 1.3373

c4 = 0.0162

[Sam Chak]

Hi @Sean,

onsidering the given constants for c1 to c4, we can observe that as u moves away from the origin point, the rate of change du/dt becomes increasingly larger. Due to the nonlinearity of the curve, the system will never reach a point where du/dt equals zero. This is why the system experiences a blow-up behavior.

Case 1: Original system

c1  = 0.1982;
c2  = 1.8115;
c3  = 1.3373;
c4  = 0.0162;
u   = linspace(0, 1, 1001);
du  = - c1*u.^1.825 + c2*u.^2 + c3*u + c4;
figure(1)
plot(u, du), grid on
xlabel('u'), ylabel('du/dt')
title('Original system behavior')

However, if we were to increase c1 by 34 times while keeping c2, c3, and c4 unchanged and positive, the system would reach an equilibrium state at u = 0.2. This adjustment would help stabilize the system and prevent it from blowing up.

Case 2: If c1 is increased 34 times...

c1  = 0.1982*34;
du  = - c1*u.^1.825 + c2*u.^2 + c3*u + c4;
figure(2)
plot(u, du), grid on
xl  = xline(0.2, '--', 'equilibrium');
xl.LabelVerticalAlignment   = 'middle';
xl.LabelHorizontalAlignment = 'center';
xlabel('u'), ylabel('du/dt')
title('If c1 is increased 34 times')

[Sam Chak]

Hi @Sean

If your goal is to make u reach zero, then ideally you would need the following conditions:

• c1 > 0,
• c2 < 0,
• c3 < 0, and
• c4 = 0.

By setting c1 to a positive value, c2 to a negative value, c3 to a negative value, and keeping c4 equal to zero, you can create the necessary conditions for u to converge to zero.

Case 3: Ideal coefficients for u(t) to converge to zero

c1  =  1;
c2  = -1;
c3  = -1;
c4  =  0;
u   = linspace(0, 1, 1001);
du  = - c1*u.^1.825 + c2*u.^2 + c3*u + c4;
figure(1)
plot(u, du), grid on
xlabel('u'), ylabel('du/dt')
title('Ideal: c1 > 0, c2 < 0, c3 < 0, c4 = 0')

Answer 3

[Sean]

Both of y'all have been incredibly helpful! I appreciate the insights. Because u represents a fluid speed, it should be bounded between 0 and much less than c, but I think I can adjust the variables that contribute to the constants to achieve this. Thanks again!