Can numerical instabilities occur when using matlab ode45 solver? If so, how can I overcome that issue?

Hello.

I would like to solve 2 unknowns with 2 equations which have the following form:

A1*x^3 + A2*x^2 + A3*x + A4 = 0; ---------- eq1

B1*diff(y,2) + B2*diff(y,t) + B3*y + B4 == 0; ---------- eq2

where, A1, A2, A3, A4, B1, B2, B3 and B4 are terms that contains constants and variables x and y. (please kindly find the code below)

The variables I would like to solve are x and y. I first solved x symbolically using eq1 with matlab 'solve' command. Since it is a cubic equation, it will produce three solutions. Say I know which one is the correct function for x, I was able to choose and then substitute it in the eq2. So, now I got second order ODE with y and t as variables. I then used ode45 to solve that equation. But, somehow the results are not making any sense.

I performed the same calculation with different software (but it is numerical finite element solver) to have some reference solutions and so that I may compare my solutions from MATLAB (kindly find the attached image for x and y results). This way I can check if my results are correct (make sense) or not.

Inside my code, I actually tried three different ways to solve the 'x' value (As have been explained before, I will substitute that symbolic equation of x into Equation 2):

1. (Trial 1) By using solve command and choosing real part of the second solution (the supposed correct solution)
2. (Trial 2) By just copy pasting the equation obtained from solve command (second solution)
3. (Trial 3) By defining a linear function for x (supposing the values of x are known in this case)

Generally speaking, all the above three equations give the same results if I introduce correct 'y' values for the corresponding 't' values. The problem is when I tried with solve command (Trial 1) , matlab is showing errors as follows:

# Warning: Failure at t=2.374146e-04.  Unable to meet integration tolerances without
reducing the step size below the smallest value allowed (4.336809e-19) at time t. #

Solving with the second equation (Trial 2) also does not give satisfactory solutions. When I tried with only linear time function (Trial 3), the calculation is running well and I got more or less satisfactory results. But, I can only do trial 3 because I already know the answers from another software. In fact, I have to solve these two equations simultaneously to get both results of x and y. I try to change the error tolerance or scale the values of y by multiplying with pi/180. It does not change anything. I still have the error. Can you explain me what is happening here OR what should I do to improve the accuracy of my results?

I don't know if my questions are clear enough. What I am thinking is that I am having some numerical problems concerning with my solver. I thank you in advance to all the suggestions and help.

Best Regards

Ye Pyae

My codes are as follows:

if true
syms x y(t)
A1 = -1.571e3;
A2 = 785.398 + 5.978e8*t*y;
A3 = 7.2734e9*t^2 - 98.1748 + 1.4945e7*t*y;
A4 = -4.1099e7*t*y + 2.1946e9*t^2;
% (1) Using solve
eq1 = A1*x^3 + A2*x^2 + A3*x + A4 == 0;
sol = solve(eq1,x,'MaxDegree',3);
x = real(sol(2));
% (2) Using sol(2) in equation form (results is truncated for the display)
% x = (597800000*real(t*y))/4713 - real((((597800000*t*y)/1571 + 6908433867456119/13818662137888768)^2/9 + .... + 6908433867456119/41455986413666304;
% (3) Using linear time function of x *assuming the solution profile is known*
% x = -2e2*t + 0.24;
B1 = 9.1892e-6*(-3*x^2+x/2+1/16);
B2 = -1.8378e-5*(3*x^2-x+1/16)/t;
B3 = 1.4906e4*((0.25+x)/(0.25-x));
B4 = 6.573e8*y*(0.25+x)*(0.25-x) - 6.3868e8*t*(2.75+10*x);
eq2 = B1*diff(y,2) + B2*diff(y,t) + B3*y + B4 == 0;
tspan = [1e-8:1e-4:1e-3];
Init = [0 0];
V = odeToVectorField(eq2); 
F = matlabFunction(V,'vars',{'t','Y'}); 
option = odeset('RelTol',1e-3,'AbsTol',1e-6,'Stats','on');
[T,Y] = ode45(F,tspan,Init,option);
X = real(Y(:,1)); 
figure
plot(T,X)
title('alpha0 Vs time'),xlabel('time (sec)'),ylabel('x (rad)')
grid on
figure
plot(T,Y(:,2))
grid on
end

[Walter Roberson]

I suspect you should be using odeFunction instead of matlabFunction(). See the workflow in the first example for odeFunction()

[ye pyae sone oo]

Thank you Roberson for taking the time to read and reply my question. I have tried with odeFunction(). In accordance with the odeFunction(), I tried to change my workflow by introducing some new lines such as reduceDifferentialOrder(eq2,y(t)), massMatrixForm(eqs,vars), etc. Then I solved my equation using ode15s solver (because I doubt if my equations are rather stiff). But, I still have the following error message:

Warning: Failure at t=2.273824e-04.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed
(4.336809e-19) at time t. 
> In ode15s (line 730)
  In Solving_for_stiff (line 39)

As a consequence, I did not get any results after the timestep 2e-4. Also, the result at that timestep is rather small (about 6 times smaller than what I am expecting). I will try to approach the problem with other solvers too.

Regards

Ye Pyae

[Walter Roberson]

Unfortunately, if there is a symbolic solution for the equation, it takes a long time to solve.

I tried asking a different package for a numeric solution, but it said that it was not able to process complex valued expressions in this context.

[ye pyae sone oo]

In fact, I first tried 'dsolve' to see if I could get any closed form expression. But, since the equations are quite complex in my case, I could not obtain any closed form results. Maybe the way I am doing is wrong here because I am trying to mix both symbolic and numeric solvers. Also I did not understand how matlab can solve the cubic equation symbolically. Anyway, thanks a lot for your reply. I am now thinking that maybe I should solve both coupled equations simutaneously without using any symbolic solvers.