CODEWITHSUNDEEP
wolfram-mathematicasolver
Duncan Mark Horne
Duncan Mark Horne

Reputation: 35

When adding functions, Mathematica ignores one when the step size is comparably large

Firstly, this problem is much easier to understand visually, I would post images but I'm new on here. I'm solving two coupled differential equations (for temperature and density) using NDSolve. The temperature equation has a heating function added to it in the form of a Gaussian and I want to thin the Gaussian (lower the variance) until it's nearly a delta function, but as I lower the variance to a certain point NDSolve starts ignoring the heating function, presumably something to do with the step size being too large? Here's the code I'm using:

Start with some jargon:

a = 1.99*10^-9;
b = 0.24*10^-3;
d = 1.21*10^-3;
T0 = 1*10^6;
n0 = 0.9*10^9;
ti = -400;
tf = 500;
kB = 1.38*10^-16;

and define the heating function "Qgt" as a Gaussian (with standard deviation "sig" and amplitude "Ag" plus some background "Qb":

Qb = 0.33*10^-3;
sig = 3;
var = sig^2;
Ag = 16.5;
Qg = Ag*Exp[-(t - 10)^2/(2*var)];
Qgt = Qg + Qb;

Then run the solver:

sss = NDSolve[{T'[t] == -(n[t]^-1) T[t]^(7/2) (a) - 
  n[t] T[t]^(-1/2) (b) + Qgt/(3*kB*n[t]), 
  n'[t] == T[t]^(5/2) (a) - (n[t]^2) (T[t]^(-3/2)) (d), T[ti] == T0,
  n[ti] == n0}, {T, n}, {t, ti, tf}];

and finally plot T[t]:

TP = Plot[T[t] /. sss, {t, ti, 400}, PlotRange -> All]

If you run all this you get a working plot where T[t] remains constant initially then the heating spike greatly increases the temperature and the system cools back to constant. BUT when I reduce "sig" to 2, the entire heating function gets ignored and the temperature never spikes up! Please run this quickly to see my problem, I don't understand what's happening and ultimately I want "sig" to be as low as 1 or 0.5. Thank you!

Upvotes: 1

Views: 122

Answers (1)

Chris Degnen
Chris Degnen

Reputation: 8655

Increasing the working precision on the solver seems to do the trick.

a = 1.99*10^-9;
b = 0.24*10^-3;
d = 1.21*10^-3;
T0 = 1*10^6;
n0 = 0.9*10^9;
ti = -400;
tf = 500;
kB = 1.38*10^-16;

Qb = 0.33*10^-3;
sig = 0.5;
var = sig^2;
Ag = 16.5;
Qg = Ag*Exp[-(t - 10)^2/(2*var)];
Qgt = Qg + Qb;

sss = Quiet@
   NDSolve[{T'[t] == -(n[t]^-1) T[t]^(7/2) (a) - n[t] T[t]^(-1/2) (b) + Qgt/(3*kB*n[t]),
     n'[t] == T[t]^(5/2) (a) - (n[t]^2) (T[t]^(-3/2)) (d), T[ti] == T0,
     n[ti] == n0}, {T, n}, {t, ti, tf}, WorkingPrecision -> 24];

Plot[T[t] /. sss, {t, ti, 400}, PlotRange -> All]

enter image description here

Upvotes: 2

Related Questions