1
$\begingroup$

I'm trying to solve the following nonlinear PDE

enter image description here

Note that when n is equal to unit, the equation becomes linear. Thus, I applied the NDSolve function to solve it numerically. Here is my code:

  kd = 1; n = 2; T = 1000; R = 1 + Sqrt[T*Pi]; 
sys = {(-D[u[r, t], r])^(1/n - 1)*((kd^(1/n)/n)*D[u[r, t], r, r] + (kd^(1/n)/r)*D[u[r, t], r]) == D[u[r, t], t]   
      ,
    Derivative[1, 0][u][1, t] == -(2^n/kd), u[R, t] == 0, u[r, 0] == 0 
   }; 
AbsoluteTiming[sol = NDSolve[sys, {u}, {r, 1, R}, {t, 0, T}, 
     Method -> {"MethodOfLines", 
       "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> T}}] ]

The R value represent the influence radius to overcome the infinite boundary condition. It has been valid compared to reference solution (n=1). However, when I'm trying to evaluate the numerical result as n > 1 (in this case the n value is chosen as 2), the massages appear:

Power::infy: Infinite expression 1/0.^1.5 encountered. >>

Power::infy: Infinite expression 1/0.^1.5 encountered. >>

Power::infy: Infinite expression 1/0.^1.5 encountered. >>

General::stop: Further output of Power::infy will be suppressed during this calculation. >>

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

General::stop: Further output of Infinity::indet will be suppressed during this calculation. >>

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >>

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >>

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >>

General::stop: Further output of NDSolve::ndnum will be suppressed during this calculation. >>

Which step I did wrong?

The reference solution is provided as the function

s[r, t]

The code is

  Vi[n_, i_] := 
     Vi[n, i] = (-1)^(i + n/2) Sum[ 
         k^(n/2) (2 k)! /( (n/2 - k)! k! (k - 1)! (i - k)! (2 k - 
                i)! ), { k, Floor[ (i + 1)/2 ], Min[ i, n/2] } ] // N; 
    Stehfest[F_, s_, t_, n_: 16] :=
     If[n > 16, Message[Stehfest::optimalterms, n];
            If[ OddQ[n], Message[Stehfest::odd, n];
                    "Enter an even number of terms",
                    If[n > 32, Message[Stehfest::terms, n];
                        " Try a smaller value for n. Maximum allowable n is 32 ",
                        Log[2]/t Sum[ Vi[n, i]*F /. s -> i Log[2]/t , {i, 1, n} ] ]],
            If[ OddQ[n], Message[Stehfest::odd, n];
                    "Enter an even number of terms", 
        If[n > 32, Message[Stehfest::terms, n];
                        " Try a smaller value for n. Maximum allowable n is 32.",
                        Log[2]/t Sum[ 
           Vi[n, i]*F /. s -> i Log[2]/t , {i, 1, n} ] ]]]  // N; 
  \[Beta]=0; rw=1;s[r_, t_] := 
     Stehfest[(2^(2 + n) r^(1/2 - n/2)
          rw BesselK[(-1 + n)/(-3 + n), -((
           r^(3/2 - n/2) Sqrt[(
            2^(1 + n) n p (\[Beta] + p \[Eta] + \[Beta] \[Eta]))/(
            kd (\[Beta] + p \[Eta]))])/(-3 + n))])/(kd p (rw^2 Sqrt[(
            2^(1 + n) n p (\[Beta] + p \[Eta] + \[Beta] \[Eta]))/(
            kd (\[Beta] + p \[Eta]))]
             BesselK[
             2/(-3 + n), -((
              rw^(3/2 - n/2) Sqrt[(
               2^(1 + n) n p (\[Beta] + p \[Eta] + \[Beta] \[Eta]))/(
               kd (\[Beta] + p \[Eta]))])/(-3 + n))] + 
           rw^2 Sqrt[(
            2^(1 + n) n p (\[Beta] + p \[Eta] + \[Beta] \[Eta]))/(
            kd (\[Beta] + p \[Eta]))]
             BesselK[(
             2 (-2 + n))/(-3 + n), -((
              rw^(3/2 - n/2) Sqrt[(
               2^(1 + n) n p (\[Beta] + p \[Eta] + \[Beta] \[Eta]))/(
               kd (\[Beta] + p \[Eta]))])/(-3 + n))] + 
           2 (-1 + n) rw^((1 + n)/2)
             BesselK[(-1 + n)/(-3 + n), -((
              rw^(3/2 - n/2) Sqrt[(
               2^(1 + n) n p (\[Beta] + p \[Eta] + \[Beta] \[Eta]))/(
               kd (\[Beta] + p \[Eta]))])/(-3 + n))])), p, t, 16]
$\endgroup$
  • $\begingroup$ You can read this post to learn how to format the code: mathematica.meta.stackexchange.com/q/1584/1871 Also, your code involves more than one simple mistake, please double check it. $\endgroup$ – xzczd Sep 16 '17 at 14:27
  • $\begingroup$ The containing code has been edited. It looks better right now. I found the boundary condition is wrong. I change the boundary condition from Dirichlet to Neumann type. $\endgroup$ – LingLong Sep 17 '17 at 0:46
  • 1
    $\begingroup$ What is rw here? $\endgroup$ – zhk Sep 17 '17 at 1:48
  • $\begingroup$ sorry for the typo. I revised it as R $\endgroup$ – LingLong Sep 17 '17 at 6:06
  • $\begingroup$ I found that it can evaluate the numerical result when n less than 1. But n larger than 1 the MMA would not give the solution. $\endgroup$ – LingLong Sep 18 '17 at 12:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.