How to solve system of 3rd order differential equations in Matlab

I want to solve the following system of differential equations in Matlab for g_a and g_b. I'm using cylindrical coordinates (r, theta) and h, ? and ? are constants. I've read the documentation but I cannot see how I can proceed. If possible, I would like to get an analytical solution - not numerical. Thank you

[Jan]

What have you tried so far? Why do you assume that this has an analytical solution?

[David Goodmanson]

Hi Emilio,

The way that this equation is written, r is not indicated as being a function of theta. In that case r would effectively be a constant. Is that correct?

[John D'Errico]

Was it necessary to remove your question? That makes the answer meaningless, and useless to anyone who might be interested in learning something.

[Rena Berman]

(Answers Dev) Restored edit

[LuisGarcia]

@John D'Errico None of the replies were satisfactory and I want to re-formulate the question.

[Torsten]

Solve the first ODE for g_a'''(theta) and the second for g_b'''(theta) and then use the usual technique to reduce higher-order ODEs to systems of first-order ODEs:

http://www.math.poly.edu/courses/ma2132/Notes/MA2132EquationsToSystems.pdf

Then use ODE15S to solve.

Best wishes

Torsten.

Answer 1

[David Goodmanson]

Hi Luis,

I am going to assume that r is constant and not a function of theta. One way that could happen is if these equations are modeling small excursions about a fixed r, is that the case?

Anyway, if r is constant, then then an analytic solution is possible. Here is an outline, without doing out the details. Lots of small steps to go through, which is typical. First, in order for the units to work out it looks like both b and mu have units of momentum^2. After multiplying the first equation by (1/h)*r^(-1/2) and the second by (mu/h)*r^(-3/2) and rearranging, you can arrive at

c1*ga' +c2*ga''' = c3*gb +c4*gb''
c5*gb' +c6*gb''' = c7*ga +c8*ga''

where the c's are uncomplicated factors that you can work out. Some of the c's include two dimensionless constants that set the problem,

    C = l/(h*r^2)   and    D = mu/h.

Since all terms include a linear factor of ga or gb, their size is undetermined within an overall multiplicative constant, which could have dimensions. Without including that constant, an analytic solution is

ga = cos(m*theta+phi)   gb = A*sin(m*theta+phi)        (to solve for m and A)

Phi can be determined later by initial conditions (maybe the overall size can be too). After plugging this in and doing the derivatives of the trig functions, you will get a couple of equations

-c1*m +c2*m^3 = A*(c3 -c4*m^2)
 A*(c5*m -c6*m^3) = c7 -c8*m^2

Multiplying the two equations together to eliminate A gives a cubic polynomial in the variable m^2. At this point if not before, you can bring in Matlab to solve this. There will be six possible values for m. For each one you can back substitute into one of the equations to find the appropriate A.

Depending on the constants it could be that m^2 has one positive root which would allow two related physical solutions m = +-sqrt(root) for ga and gb.

[LuisGarcia]

I'm planning to edit the question but your response is appropriate and rather general.

Answer 2

[John D'Errico]

But nobody knows what you have tried. Just saying that you have tried lots of things helps nobody to help you, because that forces someone to write the solution completely for you. In turn, that MIGHT take a great deal of time, because they will need to spend a fair amount of time understanding what you have, what you are trying to solve, what the variables mean, etc.

It appears that you are looking for a solution of this, wanting to find functions g_a(theta) & g_b(theta), for any value of r, as well as the other constants. Clearly there is no solution at r==0, a singularity. And we are given no information about how this process evolves if r were to vary, so effectively we can treat r as a constant, a non-zero one at that since these equations will have no solution at all at the singularity. (You might eventually look at a limiting behavior as r approaches zero though. This is also beyond the scope of Answers, as not being a MATLAB question.)

So one thing you can do is multiply (53) by r^(-1/2). That will simplify it by a fair amount. Or you might try multiplying by just r^(1/2). Again, that will simplify things a bit.

Similarly, do the same for (54). I can't know which is more appropriate, but either might work as well as the other.

I would also be looking for a non-dimensionalization of the problem.

Thus, I have NO idea what h, l, u all mean in context here. What are the units of h, l and u? I might guess that all three have similar units. In fact, I might even guess that these variables have units of length, since they appear in ways where they might also cancel the units of length on r. And since we are told that r is a coordinate variable, it has units of length.

For example, I see the term

(1 + h/u)

Since the constant 1 (one, not l) has no units at all, it makes no physical sense to add h/u to 1, UNLESS h and u have the same units, and therefore h/u is also non-dimensional, thus unit-less itself.

So I would suggest it MIGHT make sense to transform the problem in terms of some fixed ratio for h/u. I might also look at non-dimensionalizing r, so work in terms if r/l or r/h, or even r/u. The correct choice depends up on what all of these variables mean to you, in context.

All of this simplifies the problem. It MIGHT help you to understand some characteristics of the solution, much like understanding how Reynolds number impacts the solution of problems where it applies.

Does any of this apply to how you would solve the problem in MATLAB, or if a solution exists at all? Of course not, beyond possibly making the problem a bit less complex. In fact, nothing about the above comments is at all related to MATLAB. But it will simplify the problem a bit. As well, it will improve your understanding of the problem.

In the end, I might be surprised if an analytical solution exists at all. Or not. That won't happen for me without far more understanding of what you are trying to solve than I am willing to invest the time to gain.

Remember that equations, variables, etc., mean little by themselves. So not until you provide sufficient context and understanding can someone really help you here. And that makes this into something that can only become a long, drawn out consulting issue, but not a question about MATLAB.

Other things you can do are to consider other ways of attacking the problem. For example, can you represent g_a(theta) & g_b(theta) in the form of a series in theta? Think Fourier series, or Chebyshev polynomials for example, or other families of functions that would be appropriate. Again, we are given no context at all, so I am just throwing out some random ideas that MIGHT make sense in the proper context.

The point is, you have much that you can do in the way of mathematical analysis worth investing, BEFORE you just throw this at a differential equation solver and hope a solution magically pops out the end. Again though, that makes this not yet a problem about MATLAB, but a problem of mathematics.