The output of ode45 does not make much sense

Dear community members,

I am trying to solve an equation of the following form:

the solution of which must satisfy the integral constraint:

The condition is given at

My approach is to introduce a new function

satisfying the above integral constraint:

Thus, the full system becomes

As an auxiliary problem, I consider the following system:

However, when I attempt to solve it using the following code:

function dndeps = bltzmnsol(eps, n)

m = 9.1094 * 10^-31;
Np = 10^22;
Te = 300;
kB = 1.381 * 10^-23;
EN = 10^-21;
q1 = 1;
M1 = 6.6335 * 10^-26;
kappae = 1.6 * 10^-19;
ee = kappae;
ne = 10^18;

epsil = [ 0 0.03 0.05 0.08 0.10 0.13 0.15 0.18 0.20 0.23 0.25 0.28 0.30 0.38 0.45 0.55 0.68 0.83 1.00 1.98 2.00 2.10 4.23 5.10 5.95 6.00 6.05 7.95 8.00 8.05 10.10 12.20 ...
        15.20 18.20 23.10 26.20 30.90 37.20 45.60 51.00 63.00 79.00 80.00 81.00 94.00 95.00 96.00 118.00 149.00 150.00 151.00 188.00 231.00 284.00 333.00 396.00 483 584 715 895 1000 ];
csec = [ 1.26 * 10^-19 2.82 * 10^-20 1.41 * 10^-20 8.71 * 10^-21 6.21 * 10^-21 4.65 * 10^-21 3.68 * 10^-21 3.02 * 10^-21 2.54 * 10^-21 2.19 * 10^-21 1.91 * 10^-21 1.69 * 10^-21 ...
       1.51 * 10^-21 2.27 * 10^-21 3.17 * 10^-21 4.58 * 10^-21 6.67 * 10^-21 9.63 * 10^-21 1.37 * 10^-20 2.67 * 10^-20 2.74 * 10^-20 2.85 * 10^-20 6.10 * 10^-20 7.65 * 10^-20 ...
       8.93 * 10^-20 9.30 * 10^-20 9.08 * 10^-20 1.19 * 10^-19 1.22 * 10^-19 1.21 * 10^-19 1.40 * 10^-19 1.38 * 10^-19 1.23 * 10^-19 1.02 * 10^-19 7.66 * 10^-20 6.54 * 10^-20 ...
       5.42 * 10^-20 4.43 * 10^-20 3.51 * 10^-20 3.09 * 10^-20 2.53 * 10^-20 2.02 * 10^-20 2.19 * 10^-20 1.98 * 10^-20 1.72 * 10^-20 1.74 * 10^-20 1.68 * 10^-20 1.38 * 10^-20 ...
       1.11 * 10^-20 1.14 * 10^-20 1.09 * 10^-20 9.15 * 10^-21 7.63 * 10^-21 6.32 * 10^-21 5.46 * 10^-21 4.64 * 10^-21 3.83 * 10^-21 3.18 * 10^-21 2.55 * 10^-21 1.95 * 10^-21 ...
       1.70 * 10^-21 ];

sigma1 = @(eps) interp1(epsil, csec, eps, 'linear');

% figure
% plot( epsil, csec, 'o', epsil, sigma1(epsil), ':.');
% title('Cross section linear interpolation');

a1 = Np * ee^2 * EN^2/(3 * m * kappae * ne)/q1;
a2 = m * Np / ne * q1/M1 * 2 /m * kB * Te;
a3 = 2 * m * Np / ne * q1/M1 * 2 /m * kappae;

dndeps = [ n(2);  ( a1 + eps .* sigma1(eps).^2 .* (a2 - a3 .* eps) )/( 2 .* eps .* ( a1 + a2 .* eps .* sigma1(eps).^2 ) ) .* n(2) ];

end

and

eps0 = 10^-15;
l2num = 10^-15;
ninit = [ 0; 1 ];

[eps, n] = ode45(@(eps, n) bltzmnsol(eps, n), [ eps0 20 ], ninit);

the value of the integral

trapz(n(:, 2)) = 8.6458e+08

Moreover,

n(20, 1) = 5.5140e-14

Is this something to be expected? Note that ode15s, ode23s, etc. lead to a similar output.

And most importantly, is there a way to determine the value of for which

and

simultaneously?

[Walter Roberson]

Use the two-parameter form of trapz() as your eps information is not spaced at unit intervals.

[Asatur Khurshudyan]

Thank you for correcting. However, the amplitude of is too big to result in unit integral.

Answer 1

[Torsten]

function main
  n0 = 1.0;               %assume value for n(0)
  epsilonmax0 = 10;       %assume value for epsilonmax
  x0 = [n0,epsilonmax0];
  sol = fsolve(@fun,x0)
end

function res = fun(x)
  n0 = x(1);
  epsilonmax = x(2);
  fun_ode = @bltzmnsol;
  [EPS,N] = ode45(fun_ode,[0 epsilonmax],[0 ; n0]);
  res(1) = N(end,1) - 1.0;
  res(2) = N(end,2) - 1e-16;
end

function dndeps = bltzmnsol(eps, n)
  m = 9.1094 * 10^-31;
  Np = 10^22;
  Te = 300;
  kB = 1.381 * 10^-23;
  EN = 10^-21;
  q1 = 1;
  M1 = 6.6335 * 10^-26;
  kappae = 1.6 * 10^-19;
  ee = kappae;
  ne = 10^18;
  epsil = [ 0 0.03 0.05 0.08 0.10 0.13 0.15 0.18 0.20 0.23 0.25 0.28 0.30 0.38 0.45 0.55 0.68 0.83 1.00 1.98 2.00 2.10 4.23 5.10 5.95 6.00 6.05 7.95 8.00 8.05 10.10 12.20 ...
            15.20 18.20 23.10 26.20 30.90 37.20 45.60 51.00 63.00 79.00 80.00 81.00 94.00 95.00 96.00 118.00 149.00 150.00 151.00 188.00 231.00 284.00 333.00 396.00 483 584 715 895 1000 ];
  csec = [ 1.26 * 10^-19 2.82 * 10^-20 1.41 * 10^-20 8.71 * 10^-21 6.21 * 10^-21 4.65 * 10^-21 3.68 * 10^-21 3.02 * 10^-21 2.54 * 10^-21 2.19 * 10^-21 1.91 * 10^-21 1.69 * 10^-21 ...
           1.51 * 10^-21 2.27 * 10^-21 3.17 * 10^-21 4.58 * 10^-21 6.67 * 10^-21 9.63 * 10^-21 1.37 * 10^-20 2.67 * 10^-20 2.74 * 10^-20 2.85 * 10^-20 6.10 * 10^-20 7.65 * 10^-20 ...
           8.93 * 10^-20 9.30 * 10^-20 9.08 * 10^-20 1.19 * 10^-19 1.22 * 10^-19 1.21 * 10^-19 1.40 * 10^-19 1.38 * 10^-19 1.23 * 10^-19 1.02 * 10^-19 7.66 * 10^-20 6.54 * 10^-20 ...
           5.42 * 10^-20 4.43 * 10^-20 3.51 * 10^-20 3.09 * 10^-20 2.53 * 10^-20 2.02 * 10^-20 2.19 * 10^-20 1.98 * 10^-20 1.72 * 10^-20 1.74 * 10^-20 1.68 * 10^-20 1.38 * 10^-20 ...
           1.11 * 10^-20 1.14 * 10^-20 1.09 * 10^-20 9.15 * 10^-21 7.63 * 10^-21 6.32 * 10^-21 5.46 * 10^-21 4.64 * 10^-21 3.83 * 10^-21 3.18 * 10^-21 2.55 * 10^-21 1.95 * 10^-21 ...
           1.70 * 10^-21 ];
  sigma1 = @(eps) interp1(epsil, csec, eps, 'linear');
  % figure
  % plot( epsil, csec, 'o', epsil, sigma1(epsil), ':.');
  % title('Cross section linear interpolation');
  a1 = Np * ee^2 * EN^2/(3 * m * kappae * ne)/q1;
  a2 = m * Np / ne * q1/M1 * 2 /m * kB * Te;
  a3 = 2 * m * Np / ne * q1/M1 * 2 /m * kappae;
  dndeps = [ n(2);  ( a1 + eps^2 .* sigma1(eps).^2 .* (a2 - a3 .* eps) )/( 2 .* eps .* ( a1 + a2 .* eps .* sigma1(eps).^2 ) ) .* n(2) ];
end

I did not test the code, but you'll most probably need stricter tolerances for fsolve and ode45 than those that are predefined.

[Asatur Khurshudyan]

Thank you for your answer. I was getting this error:

Error using fsolve (line 300)
Objective function is returning undefined values at initial point. FSOLVE cannot continue.

I placed initial condition at

eps0 = 10^-16;

and added options to fsolve and ode45 as follows (respecively):

options = optimset('Display', 'iter', 'TolFun', 1e-30);
opts = odeset('RelTol', 1e-10, 'AbsTol', 1e-30);

This resulted in

l1num = 8.2145;

And when I used this piece

l1num = 8.2145;
l2num = 10^-15;
ninit = [ 1; l2num ];
opts = odeset('RelTol', 1e-10, 'AbsTol', 1e-30);
[epsil, n] = ode45(@(eps, n) bltzmnsol(eps, n), [ l1num eps0 ], ninit, opts);

figure(1)
plot(epsil, n(:, 1), 'b');
figure(2)
plot(epsil, n(:, 2), 'r');

to check the correctness of the root, I get

[Walter Roberson]

format long g

bltzmnsol(0, [0;1.0])

     0
     0
     1

ans = 2×1

1 Inf

Look at your code. eps is time in your code, and your final line has

dndeps = [ n(2); ( a1 + eps^2 .* sigma1(eps).^2 .* (a2 - a3 .* eps) )/( 2 .* eps .* ( a1 + a2 .* eps .* sigma1(eps).^2 ) ) .* n(2) ];

with initial time eps = 0, that code gives a division by 0.

function sol = main
  n0 = 1.0;               %assume value for n(0)
  epsilonmax0 = 10;       %assume value for epsilonmax
  x0 = [n0,epsilonmax0];
  sol = fsolve(@fun,x0)
end

function res = fun(x)
  n0 = x(1);
  epsilonmax = x(2);
  fun_ode = @bltzmnsol;
  [EPS,N] = ode45(fun_ode,[0 epsilonmax],[0 ; n0])
  res(1) = N(end,1) - 1.0;
  res(2) = N(end,2) - 1e-16;
end

function dndeps = bltzmnsol(eps, n)
  disp([eps; n])
  m = 9.1094 * 10^-31;
  Np = 10^22;
  Te = 300;
  kB = 1.381 * 10^-23;
  EN = 10^-21;
  q1 = 1;
  M1 = 6.6335 * 10^-26;
  kappae = 1.6 * 10^-19;
  ee = kappae;
  ne = 10^18;
  epsil = [ 0 0.03 0.05 0.08 0.10 0.13 0.15 0.18 0.20 0.23 0.25 0.28 0.30 0.38 0.45 0.55 0.68 0.83 1.00 1.98 2.00 2.10 4.23 5.10 5.95 6.00 6.05 7.95 8.00 8.05 10.10 12.20 ...
            15.20 18.20 23.10 26.20 30.90 37.20 45.60 51.00 63.00 79.00 80.00 81.00 94.00 95.00 96.00 118.00 149.00 150.00 151.00 188.00 231.00 284.00 333.00 396.00 483 584 715 895 1000 ];
  csec = [ 1.26 * 10^-19 2.82 * 10^-20 1.41 * 10^-20 8.71 * 10^-21 6.21 * 10^-21 4.65 * 10^-21 3.68 * 10^-21 3.02 * 10^-21 2.54 * 10^-21 2.19 * 10^-21 1.91 * 10^-21 1.69 * 10^-21 ...
           1.51 * 10^-21 2.27 * 10^-21 3.17 * 10^-21 4.58 * 10^-21 6.67 * 10^-21 9.63 * 10^-21 1.37 * 10^-20 2.67 * 10^-20 2.74 * 10^-20 2.85 * 10^-20 6.10 * 10^-20 7.65 * 10^-20 ...
           8.93 * 10^-20 9.30 * 10^-20 9.08 * 10^-20 1.19 * 10^-19 1.22 * 10^-19 1.21 * 10^-19 1.40 * 10^-19 1.38 * 10^-19 1.23 * 10^-19 1.02 * 10^-19 7.66 * 10^-20 6.54 * 10^-20 ...
           5.42 * 10^-20 4.43 * 10^-20 3.51 * 10^-20 3.09 * 10^-20 2.53 * 10^-20 2.02 * 10^-20 2.19 * 10^-20 1.98 * 10^-20 1.72 * 10^-20 1.74 * 10^-20 1.68 * 10^-20 1.38 * 10^-20 ...
           1.11 * 10^-20 1.14 * 10^-20 1.09 * 10^-20 9.15 * 10^-21 7.63 * 10^-21 6.32 * 10^-21 5.46 * 10^-21 4.64 * 10^-21 3.83 * 10^-21 3.18 * 10^-21 2.55 * 10^-21 1.95 * 10^-21 ...
           1.70 * 10^-21 ];
  sigma1 = @(eps) interp1(epsil, csec, eps, 'linear');
  % figure
  % plot( epsil, csec, 'o', epsil, sigma1(epsil), ':.');
  % title('Cross section linear interpolation');
  a1 = Np * ee^2 * EN^2/(3 * m * kappae * ne)/q1;
  a2 = m * Np / ne * q1/M1 * 2 /m * kB * Te;
  a3 = 2 * m * Np / ne * q1/M1 * 2 /m * kappae;
  dndeps = [ n(2);  ( a1 + eps^2 .* sigma1(eps).^2 .* (a2 - a3 .* eps) )/( 2 .* eps .* ( a1 + a2 .* eps .* sigma1(eps).^2 ) ) .* n(2) ];
end

[Asatur Khurshudyan]

Thank you. Yes, I have already fix that by starting from

eps0 = 10^-16;

However, the root does not seem to provide a correct answer. The shape of the solution seems reasonable, but the amplitude is too big.

[Torsten]

Start simple:

Replace your dndeps at the end of your function "bltzmnsol" by

dndeps = [n(2); 2*n(2)];

The result looks unexpected to me.

And I'm not sure whether your problem from above is really well-defined.

For

dndeps = [n(2); -2*n(2)];

e.g., you can easily see that no solution exists.

[Asatur Khurshudyan]

Yes, that is unexpected to me as well. Most probably, the first equation cannot be handled because there is no in the second equation (?). I am not sure.

On the other hand, is there a way to solve the actual problem?

[Torsten]

What I gave you is exactly the mathematical solution for this problem.

Unknowns are

n(0) and epsilonmax, 

equations are

n(epsilonmax) = 1e-16,  integral_{0}^{epsilonmax} n deps = 1

Two algebraic equations in two unknowns.

Solver is "fsolve". ODE45 is only used to evaluate n and integral n at eps=epsilonmax.

[Asatur Khurshudyan]

I run your code. It returns

fsolve stopped because the relative size of the current step is less than the
value of the step size tolerance squared, but the vector of function values
is not near zero as measured by the value of the function tolerance.

and

sol =

0.0000    8.0000

where 8 is the value I specify for epsilonmax0.

[Torsten]

I never claimed that the problem is easy and solvable in one attempt, but the implementation is correct.

I doubt that n(0)=0 is a good choice as a starting guess for n0 since this will lead to n=0 for all epsilon.

You'll have to play with the starting values for n0 and epsilonmax and with the tolerances for "fsolve" and "ode45".

[Torsten]

n(0) = 3.033620709655687e-10 yields an integral of 1 for n being integrated between 1e-16 and Infinity.

The value for epsilonmax cannot be determined for accuracy reasons. It is some value >=5 (where n can no longer be distinguished from zero).

[Asatur Khurshudyan]

I came to almost the same conclusion, which is good, because that means we can consider, let's say, and both conditions will be satisfied.