Hi!

My first attempt at solving this system would be to use the Keller box scheme. The details of how this ought to be implemented will depend on what your boundary conditions are and where they are. Things could get difficult if the origin is part of the solution domain, so I cannot say more until I know if you are solving in a circle or in an annulus.

I would also use the Keller box scheme in a way which is reminiscent of Crank-Nicolson, which is to discretise the first and second r-derivatives as they stand, rather than by converting this into 8 first order equations in r. The code will run more quickly and still have 2nd order accuracy in r.

I have noticed that there are some extra time-derivatives in a few of the given equations. That can be coped with if using the Keller box scheme.

Given that the system is nonlinear, it will be necessary to find the Jacobian matrix in order to create the Newton-Raphson iterations scheme. Personally, I use numerical derivatives to do this....oter wise you'll have to find this linearisation by hand.

As a young researcher, it took me a long time to create such a general code in Fortran, and many of the "extra features" (such as numerical derivatives, using the equations in 2nd order form) were the result of new thinking over the years.

The change in the amount of encoding in going from explicit to implicit is huge for your system. So I am wondering why you haven't considered using an explicit time-stepping method? Does that exponential cause severe numerical instability problems? One could use Euler timestepping to get from t=0 to t=tmax using dt as your timestep, and then do it again with dt/2 as your timestep, and then use Richardson's extrapolation technique to create second order accurate data in time. Encoding an explicit scheme will be very fast, and with a tiny bit of extra processing, it will be possible to get second order accuracy in time.

Thank you our Prof D. Andrew S. Rees.

In addition you can also use the spectral method.

Dear Prof Wakif,

Do you mean with zeroth order Bessel functions, J_0 and possibly Y_0 if the domain is an annulus? If so, then my thought is that the exponential term, e^{-1/h), could prove to be extremely awkward to write in terms of a sum of Bessel functions. Well, I don't know how to do it!

Regards.

I mean spectral method using chebyshev functions with the gauss lobatto collocation points

Well, I don't think the problem is in the formulation for the spatial discretization but is focused on the time integration.

I agree that the equations can be arranged to adopt a Crank-Nicolson formulation and a suitable linearization.

This problem can be solved with spectral method for space variable and Runge kutta algorithm for temporal variable

Dear all,

thanks a million for your helps. I will try Keller box method.

....although I said that my first attempt at a numerical solution would be to use the Keller box method, I have to say that the only reason I said that was because I have my own near-black-box form. Writing a KB code from scratch is quite involved. If you don't know how the solutions will behave in time, it would be worth using a simple forward Euler in time so that you will then have something to compare with when trying to use a more complicated method. Forward Euler also has the advantage that the dg/dt in the equation for f is known from the equation for g. And the df/dt terms in the third and fourth equations will be known from the equation for f.

What are your boundary and initial conditions?

بالتوفيق سيد سعيد

My best wishes too!

I think that the program is useful in solving non- linear PDE's .

Dear Prof Fuliful,

Your code appears to solve a Laplacian in one variable, which is a linear elliptic steady 2D system, whereas the question was about solving four coupled 1D-in-space nonlinear parabolic equations. Or am I mistaken?

Hello Saed Moradi

It would be preferable and much simpler to solve this problem with the finite element method with the option of using an explicit scheme or implicit schem or Crank-Nickolson scheme.

For a person unused to FEM the only problem is how to form the numerical matrix time marching scheme. But that is quite easy.

Congratulation.

Nice work, sir! Thank you for the feedback.

Sir can you send me your code so that i can solved my similar equation. since i have similar equations but i could not solved.

Md. Roknujjaman See the attached file (Mathematica notebook).