Simple ODE solution from MATLAB and Mathematic disagree

Hi, I have in the last 2 months been trying to solve the rather simple ODE:

aY''(x) + ibY'(x) -c = 0

initial condition Y(0) = B

and after checking, MATLAB and Wolfram alpha give two different solutions:

Can someone point me to a reproducible solution to this simple ODE ?

here are the respective codes and outputs:

MATLAB Online:

if true
  % code
end
syms a b c Y(x) x B
eqn = (a)*diff(Y,x, 2) + (i*b)*diff(Y,x) == c;
cond = Y(0) == B;
Y(x) = dsolve(eqn, cond)

ANSWER: C + exp(-(b*x*1i)/a)*(B - C11) - (c*x*1i)/b

Wolfram alpha online input:

a*y''(x)+i*b*y'(x) -c = 0 initial condition y(0)=B

output:

(e^(-(ibx)/a)*(e^((ib*x)/a)*(ib*B - i*a*x) - i*a*C*e^((ibx)/a) + i*a*C))/c

As you can see, these are very different, and have completely different limits and behavior.

What is the true solution to this simple ODE?

and why can't such a simple procedure agree between two major programs as these used? Thanks!

Answer 1

[John D'Errico]

First of all, I sincerely doubt that is the true output directly from Alpha. I think you did some editing. In some places I see ib*x, in some I see ibx, in some places I see i*a*C. I'd be careful with edited results, since they might contain bugs.

In fact, if I paste in the code you wrote into MATLAB (R2017b), I get a somewhat different result!

syms a b c Y(x) x B
eqn = (a)*diff(Y,x, 2) + (i*b)*diff(Y,x) == c;
cond = Y(0) == B;
Y(x) = dsolve(eqn, cond)

Y(x) =
C11 + exp(-(b*x*1i)/a)*(B - C11) - (c*x*1i)/b

So it looks like you did some editing here too, changing at least the first C11 into a C. If you will claim the results to be what you got, you need to learn to use copy and paste.

Next, these look different, but they are not that different. That exp(-i*b*x/a) term out front kills off the exp(i*b*x/a) terms inside the parens. So you have no exp(i*b*x/a) terms at all. The Alpha expression now looks way more like the MATLAB one. In fact, it appears that it will have very much the SAME behavior.

You should recognize that computer algebra systems don't always choose the same forms as might you. They don't always see obvious simplifications, at least obvious ones to your eyes.

Next, you have a second order ODE. But you posed only ONE initial condition! So effectively the ODE solver needs to make some arbitrary choices. They will potentially be different choices in the solution. For example, if I provide no initial condition at all, then dsolve tells me:

Y(x) = dsolve(eqn)

Y(x) =
C12 + C13*exp(-(b*x*1i)/a) - (c*x*1i)/b

This is in fact pretty much what Alpha produces, when posed with that ODE, but no initial condition at all.

So two unknowns. When you supply only one condition, the solver will decide arbitrarily which of the unknown parameters to eliminate.

I think the problem is not with the solvers, but how you are using them, and how you are viewing the result, and even whether you were careful in editing the result. (Lack of care in one place may equate to lack of care in how you used the solver.)

[Sergio Manzetti]

Thanks John. I can't say the editing made such a big problem, I was trying to make it easier to read here.

In case for the conditions, having only one condition is fine, because the next condition is treated by normalization.

About the behaviour: When I plot these two results, they look completely different in MATLAB.

[John D'Errico]

The editing is only an issue (and a minor one) because I can't be sure what you are comparing.

As far as the initial condition goes, while you say you are normalizing things, again, I can't see what you did. You can't normalize a curve with an unknown parameter in there. You certainly cannot plot it. So I am not at all positive what you did. And if I don't know exactly what you did, then I can't reproduce it to see what you are seeing, and to know what is the problem. I am fairly confident that the two solutions are effectively the same, even without doing the algebra, because just looking at the Alpha solution, if you do some basic simplification, the two look quite similar.

I would suggest that instead of post normalizing things as you claim to be doing, that you might consider providing that as a condition for the solver.

One way to test that both are producing valid results is to take the MATLAB solution, and substitute it into the ODE. Then do the same with the Alpha solution. If you don't get zero both times, then they are not solutions to the same ODE.

And By the way, be very careful in the numerics here. Depending on the values of a and b, there might be numerical problems with those exponentials. So it is possible that your normalization is not doing what you expect it to have done. Again, I can only hazard wild guesses here, because I have insufficient information to know anything.

[Sergio Manzetti]

Thanks again. In fact I Have been having quite a lot of problems with the normalization of the C. B is given by some initial condition, i.e. \hbar. C works, nevertheless I had some problems that could only be fixed by restarting MATLAB on the machine, then it gave the expected C. I will give this a further thought!

Thanks !

[Sergio Manzetti]

Dear John, is there a way the numerical normalization can be checked and verified by someone here at MATLAB ? Its for a paper so I'd prefer not to post it on the open forum. The reason is that, as you mentioned, I experienced some problems from day to day on the same code, and had to reboot matlab and got a different answer.

[John D'Errico]

Restarting MATLAB and getting a completely different result is not a MATLAB issue, but always a question of that you changed something when you restarted. Some parameter, something in the equations was different. I recall doing an lengthy analysis (long ago, in APL) where I got a really pretty result. But then I came back and tried to replicate what I had done, and I never could do so. I must have done something different the first time. For example, you can use a script, putting all of your work in the script. Then running it a second time will give fully repeatable results. You can split the script into sections by %% lines. Then only need to build one section at a time until your code works. I frequently use such a method when I am doing some analysis, or working on a new idea.

Unfortunately, while you can use the consulting services provided by TMW, I don't think this is a free service. And Answers is a fully public forum, not a bulletin board where you might advertise for someone to provide free (but private) consulting.

Answer 2

[Sergio Manzetti]

Yes, that is what I experienced too. However, this is actually a bug, and I have reported it to MATLAB Support. It happened also with another command during the same session.

Thanks!