Change of boundary condition

9 visualizaciones (últimos 30 días)
Hai Nguyen
Hai Nguyen el 19 de Sept. de 2020
Editada: Hai Nguyen el 5 de Oct. de 2020
Attached is a long code of a mathematical problem in modeling.
If I want to change r(1)=ya(2)-t^alpha with alpha given, can you please help define t so that it is the same as t in the main loop while?

Respuesta aceptada

Walter Roberson
Walter Roberson el 25 de Sept. de 2020
Change
Phig=gamma*(Vol+2*C*Leak)-t;Phib=GS(y(:,1),y(:,N),tipVal);
to
Phig = gamma*(Vol+2*C*Leak) - t;
Phib = GS(y(:,1), y(:,N), tipVal, t);
Change
function r=GS(ya,yb,tipVal)
r(1)=ya(2)-1;r(2)=yb(1)-tipVal;end
to
function r = GS(ya, yb, tipVal, t)
alpha = 1/3; %or 3/8. Or whatever is appropriate
r(1) = ya(2) - t.^alpha;
r(2) = yb(1) - tipVal;
end
  8 comentarios
Hai Nguyen
Hai Nguyen el 27 de Sept. de 2020
I got it now. Thank you for your clear explanation.
Walter Roberson
Walter Roberson el 28 de Sept. de 2020
I do not have a complete program for you yet, but I decided I would post the current version of the code for you.
I implemented looping over alpha values.
I also did investigations of some of the failures, and put in smooth controls when some of the problems are detected, as I was able to track down where some of the problems happened, even if I do not yet understand why in all cases. So far the controls seem to have caught all of the cases of (for example) trying to spline over complex values.
My investigations appeared at first to show that there was some sharp cut-off of alpha values, that any alpha up to some particular value about 0.4 worked, and that every alpha larger than that led to failure. However, when I narrowed down the failure point down to the 0.0001 level, I did find that there are some alphas that work after some alphas have failed. The range over which this happens is pretty narrow, though, and in my testing so far, it at least looks like there is some alpha value beyond which nothing works.
The current program outputs summary information that can look similar to this:
Jacobian was dubious sometimes, rcond down to 3.73266e-13, but proceededed using \ operator anyhow; at alpha = 0.401732
Alpha = 0.401732 succeeded!
b went complex at t = 0.419431 alpha = 0.401733, iter = 5
No convergence at t = 0.419431
Jacobian was dubious sometimes, rcond down to 3.73266e-13, but proceededed using \ operator anyhow; at alpha = 0.401733
Alpha = 0.401733 stopped early due to failure, reason = "b went complex"
and
lambda_u exceeded at t = 0.0262153 alpha = 1.96, iter = 10
No convergence at t = 0.0262153
Jacobian was near singular, rcond down to 4.19951e-15, substituted pinv for \ operator at times; at alpha = 1.96
Alpha = 1.96 stopped early due to failure, reason = "lambda_u exceeded"
There are several key points about this:
  • there are alpha values that succeed with no warning messages needed
  • there are alpha values that succeed with warning -- the warning tells you that something very questionable is going on with the mathematics and yet something at least looks like success came out of it
  • there are alpha values for which convergence fails because the b matrix for J\b becomes complex. When b becomes complex, then that "pollutes" the remaining calculations and would inevitably lead to an error about trying to take spline over complex values. There were no cases in which the calculation was able to succeed after b became complex
  • there are alpha values for which convergence fails because lambda_u exceeds the built-in safety point that apparently has something to do with the equation being considered to rise too rapidly (I do not know what that is about)
  • there are alpha values for which the jacobian ends up with a reciprocal condition number that is small enough to be of distinct concern, that a single bit of round-off difference could result in a change about 3E12 larger than the value associated with the bit. Those cases happened regularly enough that the entire calculation chain should be reviewed as being questionable. But still, I was able to find a boundary where when the reciprocal condition number was above that boundary, the calculations sometimes succeeded. For this range of reciprocal condition number, I just set up a warning, but still have the code path use the J\b calculation.
  • there are alpha values for which the jacobian ends up small enough that MATLAB was warning about likely failure of the J\b operation. I was able to find a boundary where when the reciprocal condition number was smaller than the boundary, that the overall calculation has never succeeded. Failure was not typically immediate in such cases, but left to enough iterations, there has always been failure if the reciprocal condition number was smaller. However, when I detect that case, I do not terminate iteration, but I do avoid the singularity warning by switching that one operation from J\b to pinv(J)*b
  • 20 iterations was too small. I changed the limit. There were some calculations going at least as high as 55 iterations, and successes over 30.
  • So far, every alpha 0.401767 or higher has failed. 0.401766 succeeded even though some entries slighty lower failed.
  • In the case where jacobian (near-) singularities were detected, the failures distribute fairly evenly between exceeding lambda_u or b going complex (I thought I had detected a pattern that lambda_u was "nearly always" the failure cause in this situation, but further testing showed that they are pretty balanced.)
As to why b sometimes goes complex: I traced it down to some calculations of the form sqrt(b - spline interpolation involving gamma and tau and different x locations), which produces a complex value if the interpolated value of the spline exceeds b. The calculation involved three different terms similar to that, added together; by splitting up the code my experience was that the first and third term both produced complex values and (so far in my testing) the second term was real-valued.
I have not tracked down the circumstances under which the interpolation gets too large. In at least some circumstances it gets triggered when the new gamma value created in the main loop and inserted into the last gammah vector entry, is lower than the second-last entry. As the gammah values are effectively x values in the interpolation, that corresponds to having a series of x values that doubles back on itself, leading to a distinct hook in the gammah vs tau plot. My testing earlier seemed to indicate that some doubling back was tolerated by the calculation, that the spline did not end up projecting a value too large, but that the doubling back "too much" was a problem in other cases.The next step would be to examine that more closely to try to figure out what was safe and what was not.
I am not comfortable with the robustness of the algorithm. Reciprocol condition numbers do not have exact failure boundaries, but it is generally considered that by the time you get down to 1e-10 that your algorithm has probably failed (design error ==> needs more robust code at the very least, but quite possibly that the algorithm is Wrong) -- or else that your equations have a singularity in that area. rcond like 3.73266e-13 tells you that errors are being magnified by a factor of a trillion. You should not trust any answers involving such calculations. And yet the algorithm produces those jacobians often, including for alpha = 0 (equivalent to the original version of the code with no time-dependent modification.) The algorithm needs a good going-over.
I am also not comfortable with the way it handles alp value, switching between 1/3 and 3/8 as it goes, but passing the last one it happened to be using into the plotting routine as-if it had always been that value. And I am not comfortable that there is a definite difference between that alp value and the alpha value that you asked to have put in: I worry that they are really the same factor and that the code should not be switching between 1/3 and 3/8 on the fly.
Current code version is attached.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Performance and Memory en Help Center y File Exchange.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by