Issue with solving ill-conditioned symbolic DAE system

Kindly I want to know if there any changes in some solvers like ODE113 or symbolic engines?

Recently I have updated matlab to 2023b version and tried to run a script that I have excuted it before with nothing changed on older version of matlab, the script contains solving a symbolic DAE system Ax = B, with ill-conditioned A matrix, using ODE113 of an older version it was working with conditiong number was around e-21. my results using older version was compared to some other CAEs (figure below) and they were identical, so I do not know what is the issue in 2023b, any help please?

[Torsten]

What MATLAB version did you use prior to R2023b ? I doubt that changes have been made to ode113.

Do you mean you want to solve a DAE system of the form A*xdot = B instead of A*x = b as you state ?

[John D'Errico]

I assume the A*x==b is the algebraic part of the DAE.

[Mohamed Elshami]

@Torsten I was using 2019b to solve symbolic DAE system resulting from Multibody formulation, I found after a lot of trials that they should updated the symbolic engine, I am sure about this conclusion as I returned back to the older version and it works well even the model is ill-conditioned and using the same ode133 solver.

Answer 1

[Mohamed Elshami]

Hi there, If any one inerested. The problem was not the solver type, instead the modeled system was over constrianed, after research I solved this problem by reducing the number of constraints equations by following another concept to remove any redundancy in DAEs set.

[Torsten]

But you said it worked well in the old MATLAB version ...

Answer 2

[John D'Errico]

You have a DAE, with a linear constraint system that happens to be ill-conditioned. I would try reducing the constraint system. You lose nothing in the process. That is, first verify the constraints are consistent. Thus you want to have

rank(A) == rank([A,b])

If not, then you have a fundamental problem. But assuming the above test succeeds, then you can reduce the problem, so the algebraic part is no longer ill-conditioned. (Well, hopefully. That sort of depends on where the ill-conditioning arises, and how.)

I'll give a simple example of a problem with an ill-conditioned constraint.

A = randn(3,2)*randn(2,5)  % Just a randomly generated constraint set.

A = 3x5

0.5715 0.0694 0.3485 0.7862 0.2419 -0.2496 -0.2250 0.2155 -0.2810 0.0926 1.4721 0.8320 -0.3361 1.8158 -0.0422

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

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

Note the rank of A is 2, even though there are 3 constraints.

size(A)

ans = 1x2

3 5

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

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

rank(A)

ans = 2

I'll pick an arbitrary right hand side that is consistent.

b = A*[2;3;5;7;11]

b = 3x1

11.2579 -1.0451 16.0055

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

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

The next test I would perform is to verify the problem is indeed consistent. (Even though I claim that to be the case, you may not believe me.)

Arank = rank(A)

Arank = 2

rank(A) == rank([A,b])

ans = logical
   1

So the algebraic part is indeed consistent, and the rank of A is 2, as I would expect. Next, look at the singular values of A.

svd(A)

ans = 3x1

2.6977 0.6082 0.0000

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

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

As you can see, one of them is massively smaller than the others. That tells me there will be no issue in reducing the problem.

[Q,R] = qr(A,0)

Q = 3x3

-0.3575 -0.8837 -0.3022 0.1561 -0.3756 0.9135 -0.9208 0.2794 0.2722

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

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

R = 3x5

-1.5988 -0.8260 0.2186 -1.9968 -0.0331 0 0.2557 -0.4828 -0.0819 -0.2603 0 0 0.0000 -0.0000 0.0000

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

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

The trick now is to reduce the problem as:

Ahat = R(1:Arank,:);
bhat = Q'*b;
bhat = bhat(1:Arank);

Ahat will no longer be singular. It may still be moderately ill-conditioned though, and that will depend on how A was created, and the source of the ill-conditioning. But you can now use the reduced constraint pair Ahat and bhat instead of A and b.

[Mohamed Elshami]

Thank you very much! your comment enriched my knowledge by your valuable provided information, After a lot of trials I found that the issue is in the symbolic toolbox it seems to be handled completly different in versions newer than 2019b, I returned to that vesrion to check if I forgot somthing and unexpextedly the code worked as first time I wrote it, however itself still not working in newer version 2023,2024.

my DAE system includes a lot of symbolics in A, B matrices that would be first inversed and multiplied ...etc.