Matlab gives [ empty sym ] when using dsolve

I have two differential equations which are complicated. Later I will substitute the solution into another differential equation, so I need to keep the analytical form.

syms Ap(z) Ac(z)

ode1 = diff(Ap) == (Ap*(Ac^2*conj(Ac)^2*7.12e+27i - Ac*conj(Ac)*(2.45e+29 + 4.89e+28i) + ...
    Ac*Ap*conj(Ac)*conj(Ap)*7.12e+27i)*1.53e+34i)/(Ac^3*conj(Ac)^3*5.07e+55i + ...
    Ac^2*Ap*conj(Ac)^2*conj(Ap)*1.52e+56i - Ac^2*conj(Ac)^2*6.97e+56i + Ac*Ap^2*conj(Ac)*conj(Ap)^2*1.52e+56i + ...
    Ac*Ap*conj(Ac)*conj(Ap)*1.39e+57i + Ac*conj(Ac)*6.22e+58i + Ap^3*conj(Ap)^3*5.07e+55i + Ap*conj(Ap)*5.98e+58i);
ode2 = diff(Ac) == -(1.53e+34*Ac*Ap*conj(Ap)*(2.45e+29 + Ac*conj(Ac)*7.12e+27i + ...
    Ap*conj(Ap)*7.12e+27i))/(5.07e+55*Ac^3*conj(Ac)^3 - 6.97e+56*Ac^2*conj(Ac)^2 + 5.07e+55*Ap^3*conj(Ap)^3 + ...
    6.22e+58*Ac*conj(Ac) + 5.98e+58*Ap*conj(Ap) + 1.39e+57*Ac*Ap*conj(Ac)*conj(Ap) + ...
    1.52e+56*Ac*Ap^2*conj(Ac)*conj(Ap)^2 + 1.52e+56*Ac^2*Ap*conj(Ac)^2*conj(Ap));

odes = [ode1; ode2];
cond = [Ap(0) == 1.3105, Ac(0) == 13.1046];
dsolve(odes, cond)

I want to solve "Ap(z)" and "Ac(z)", but I get a empty solution:

ans =
[ empty sym ]

How can I do for a such problem? Thanks very much for any advice.

Answer 1

[Samay Sagar]

In cases where the dsolve function in MATLAB returns an empty solution, it means that it was unable to find an analytical solution for the given set of differential equations. This can happen when the equations are too complex or do not have a closed-form solution.

In such situations, you might consider using numerical methods to approximate the solution instead. MATLAB provides a variety of numerical ODE solvers that can handle complex systems of differential equations. You can use the ode45, ode23, or ode15s functions, among others, to numerically solve the equations.

[John D'Errico]

@Samay Sagar

You are completely correct, in that an empty solution does mean dsolve was unable to solve the problem. (And you got an upvote for that from me!)

However, you might also have pointed out that a numerical solve will also almost surely fail. We see constants in those equations that cover a HUGE range of magnitudes. For example, 6.22e+58i, 7.12e+27i, 5.07e+55, etc.

The problem is, these numbers, over a large dynamic range means that any work in double precision arithmetic will almost surely be a problem. As such, I would be very surprised to see a valid numerical solve, and even if tools like ode45 do return some solution, the solution would not be something I would trust without a great deal of time spent in analyzing it. At the very least, this problem would certainly be a stiff one. (Stiff systems are often characterized by problems with widely varying rate constants.) So tools with names that end in s would be the first choice, and even then, I would be leery.

[I-Chun]

@Samay Sagar, @John D'Errico

Ah, that makes sense. I should not be attached to finding analytical solution. Thank you very much.

[Steven Lord]

To restate John's second point, where do those very large magnitude coefficients come from?

Perhaps if this is a physical problem a change of units might be in order (using meters instead of millimeters or vice versa.)

If those came about from fitting a polynomial curve to data, what degree polynomial did you fit? If you had say 100 data points you could fit a degree 99 curve to that data, but that doesn't mean you should do that.

x = 0:100;
y = sind(x);
p = polyfit(x, y, numel(x)-1);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Note that the plot of p takes on values over 1 when x is larger than say 80. But the sine function's value for real inputs is between -1 and 1 inclusive.

plot(x, polyval(p, x))

Perhaps use a lower degree curve. Note the strong agreement between the actual data (the + symbols) and the fitted curve (the o symbols.)

p2 = polyfit(x, y, 5);
figure
plot(x, y, '+', x, polyval(p2, x), 'o')