How to solve a nonlinear differential equation of this type

Question

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.

NOW MY PROBLEM:

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.

%%%UPDATE%%%%

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

$$\left|T'+A(T)+B(T)\right|=f(T)$$

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!

Answer 1

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.

Answer 2

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?