Trying to explain my problem better, I need to start with the diffusion Equation, which is the one that I am trying to solve.

The Diffusion Equation

\begin{equation} \dfrac{\partial u}{\partial t}=\dfrac 3 x\dfrac \partial{\partial x}\left(\sqrt x \frac \partial {\partial x}(\nu(T(u,x))u\sqrt x)\right) \end{equation}

Here $u$ is the surface density of a gas disk. $T$, is the temperature of that disk. My aim is calculate the evolution of that disk using an initial condition. After treat this equation with MOL (Method of lines) in Matlab I am able to solve this equation. But for that I need to use another equations that will relate some of the variables in my problem, so finally I have:

\begin{equation} \nu\propto \frac{T}{\Omega } \end{equation}

$\Omega$ is just the keplarian velocity of the disk, therefore is a function that depends on $x$, and also I have

\begin{equation} \Xi(T,\Omega)T^3=u^2\Omega \end{equation}

and $\Xi$ is another function which depend also in $T$ and $\Omega$

Then, since I have an initial condition for $u$ and of course I define my dominion in $x$ ($x=1..10^3$ for instance). I just need to solve the last equation and find the value for T given some position and a value of $u$. I do this with a method to find zeros.


I need to add advection to my problem. In terms of equations I have this now:

\begin{equation} \Xi(T,\Omega)T^3=u^2\Omega + \frac{T^4\;u\;H(T,\Omega)}{\Omega x^2}\left|\frac{d\;\log\;H(T,\Omega)}{d\;\log\;x} - \frac{d\;\log\;T}{d\;\log\;x} - \frac{d\;\log\;u}{d\;\log\;x} \right| \end{equation}

\begin{equation} H^2\Omega^2-H\frac{T^4}{u}-T=0 \end{equation}

But I only need one root (the real and positive one), and of course the method to solve this (the one that I'm trying to find) will be the same using any root of H.

My problem here is the derivatives, without them, I can use the same method that before to find T. But now I can't do that. Because I need to calculate those derivates.

Please ask me anything about the problem.


Well, after doing a "little" bit of algebra, I realize that the equation is of this type:


So it can be solve numerically as a regular ODE. With two possible solutions

$$T'+A(T)+B(T)=f(T)$$ and, $$T'+A(T)+B(T)=-f(T)$$

Thank you very much for all the help!

  • $\begingroup$ What's the application domain of this problem? $\endgroup$ Commented May 10, 2012 at 20:37
  • $\begingroup$ Well T goes from 0.. finite value. (Is that what you ask?) If not please tell me. $\endgroup$
    – Nikko
    Commented May 10, 2012 at 20:39
  • $\begingroup$ @Nikko: I think he's asking what physical phenomenon you're using this equation to model. $\endgroup$
    – Dan
    Commented May 10, 2012 at 21:13
  • $\begingroup$ Oh ok then. I am trying to calculate the temperature of a gas disk, given the surface density of the disk. This equation will give me that temperature with the initial surface density profile of the disk. $\endgroup$
    – Nikko
    Commented May 10, 2012 at 21:21

2 Answers 2


Your equation is an implicit ordinary differential equation.

Depending on the boundary conditions:

  • If you have an initial value problem, you could use a code that solves differential-algebraic equations like DASSL, DASPK, or IDA. These problems all use BDF methods; in particular, IDA has a MATLAB interface, since it is part of SUNDIALS. In the community of people I work with, these methods tend to be the most "standard," because BDF methods handle stiffness well, and I tend to work with stiff problems.
  • If you have an initial value problem, you could use also use other numerical methods, like Radau collocation, or spectral collocation.
  • If you have a boundary value problem, in principle, you could adapt methods for explicit ordinary differential equations to implicit ordinary differential equations, according to Numerical Solution of Boundary Value Problems for Ordinary Differential Equations (See section 11.2 for details). It's not entirely clear how to follow the authors' advice; presumably, a literature search would turn up appropriate methods.
  • You could treat it as a boundary value problem for differential-algebraic equations; again, you'd have to do a literature search to find appropriate methods. Some type of collocation should work here also. One pertinent code is COLDAE, by Ascher, et al.

You'd have to supply more information for me to give you more specific advice.

  • $\begingroup$ I edit the information about my problem, is it better now? $\endgroup$
    – Nikko
    Commented May 12, 2012 at 17:48
  • $\begingroup$ It's somewhat better. I'd still like to see the set of differential and algebraic equations that you pass to an appropriate solve, along with initial conditions, so that I can give you better advice. I apologize for the late response; my hard drive went down this past weekend, and I spent a lot of time repairing my computer and restoring the data from a backup. $\endgroup$ Commented May 14, 2012 at 23:00
  • $\begingroup$ It's ok. Thank you. Actually I'm really close to solve it now. After some algebra the equations turns out to be this form: $T'=f(T,x)$ The only problem is I only have boundary conditions and no initial conditions. So IDA will do this right?. I am reading also the BVP solver of matlab. What do you think is better? $\endgroup$
    – Nikko
    Commented May 15, 2012 at 11:05
  • $\begingroup$ IDA will only solve initial value problems (IVPs). In principle, you could implement a method like multiple shooting to solve a boundary value problem (BVP) using an IVP solver, you're probably better off calling a BVP solver directly. Depending on your problem, the MATLAB BVP solvers may or may not fit your needs. Try them out; if they don't fit your needs, you should hunt around for a better BVP solver that will give you sufficiently good results. $\endgroup$ Commented May 15, 2012 at 17:56
  • $\begingroup$ The thing that I'm confused. My equation is highly non-linear (because of f(T,x)). But since is of first order (only with T') maybe is an IVP for x or not? $\endgroup$
    – Nikko
    Commented May 15, 2012 at 20:56

Instinctively, I'd say use Chebfun, which uses spectral collocation and automatic differentiation to solve non-linear BVPs. In any case, I'd replace $R$ with $\exp(r)$ to make the derivatives more reasonable.

I'd try it myself, but there is quite a bit of information missing, e.g. what boundary conditions are you using and what are your constants?

  • $\begingroup$ Thanks Pedro!. Ok so the thing is that is not obvious to me the boundary conditions. Because what defines it is on the G(R) function. $\endgroup$
    – Nikko
    Commented May 10, 2012 at 20:52
  • $\begingroup$ I try to improve the information of my problem. $\endgroup$
    – Nikko
    Commented May 12, 2012 at 17:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.