## 245 Reputation

19 years, 32 days

## The Matlab Result...

I am unable to repeat the Matlab result, and if I understand the problem correctly, for the given input system, the Matlab result must be wrong. Once again, back to the ODE argument. For the Matlab solution, it appears as though the value P(0) is in the region of 0.135. Now when I apply a shooting method with P(0) ~= 0.135, and D(P)(0)=0, the resulting plot is actually oscillitory, it does not vary smoothly as the Matlab plot does. In fact, if I obtain plots for P(0)=0.134 and P(0)=0.136, there is very little difference between them. Both functions are oscillitory, have 4 peaks and 5 valleys, and neither has D(P)(8e-8) even close to 0. Now since the ODE does not become undefined in the region z=0..8e-8, there is a unique solution for a given value of D(P)(0) and P(0) for this problem, and since there is little variation around P(0)=0.135, then the Matlab solution is not a solution of the IVP taken from z=0, therefore it is not the solution to the problem. The fact that there is such a dramatic difference makes me believe that either: 1) The solution in Matlab is an artificial solution introduced by the discretization. or 2) There is a very small error (perhaps a sign error) in the setup of the problem. Sorry I was unable to help, Allan

## Double Check...

Just to confirm then, the *exact* problem you are solving is: A:=abs(3.0164e14*z-2.8071e7)*(-0.82500e-2*z+3e-9)^2: alpha_z:=3.0164e14*z-2.8071e7: beta_z:=9.9000e15*z-5.400e8: gamma_z:=6.6e9: # (taken as constant for now) bc := D(P)(0)=0, D(P)(8e-8)=0; deq := A*diff(P(z),z,z)-alpha_z*P(z)-beta_z*P(z)^3-gamma_z*P(z)^5=0; Is every coefficient and term in the above formulation correct to the specified number of decimal points? It takes some time to thoroughly analyze these problems, so getting it right would be helpful. Thanks, Allan

## No Worksheet...

I don't seem to be able to upload the worksheet to the site, so here are the commands I used to demonstrate that the problem is single valued: # Define problem A:=abs(3.0164e14*z-2.8071e7)*(-0.82500e-2*z+3e-9)^2: alpha_z:=3.0164e14*z-2.8071e7: beta_z:=9.9000e15*z-5.400e8: gamma_z:=6.6e9: # (taken as constant for now) bc := D(P)(0)=0, D(P)(8e-8)=0; deq := A*diff(P(z),z,z)-alpha_z*P(z)-beta_z*P(z)^3-gamma_z*P(z)^6=0; # Generate a solution for D(P)(0)=0 and some arbitrary value for P(0) Digits := 30: dsn := dsolve({deq,D(P)(0)=0,P(0)=0},numeric,abserr=1e-20,relerr=1e-20): # Create a function that uses dsolve/numeric to compute the value D(P)(8e-8) fr> fcn := proc(P0) local val; if type(P0,'numeric') then dsn(initial=[0,P0/1e16,0]); try val := op([3,2],dsn(8e-8)); catch: val := Float(undefined); end try; val; else 'procname'(P0); end if; end proc: # Plot this function over several ranges: plot(1e13*fcn(z/1e6),z=-1..1); plot(1e10*fcn(z/1e3),z=-1..1); fcn(1e-3); fcn(1e-2); fcn(1e-1); fcn(1e0); fcn(1e1); # So from the plots, and from the above trend, it appears as though # the only value of P(0) that gives D(P)(8e-8)=0 is 0

## Solving 'A' w.r.t. 'z'...

The reason I did this is as follows: The ODE in question is leading linear, and the coefficients are purely polynomial in the interior of the region 0..8e-8. Since this is true, if the coefficient of the leading derivative of the problem is non-vanishing, then the initial value problem over the specified range is well defined. This suggests that use of a shooting approach will be a possibility for obtaining a solution. When I attempted to apply a shooting approach, varying the value of P(0), (with the incorrrect value for 'A'), I obtained a problem where the value of D(P)(8e-8) varied monotonically w.r.t. P(0), and the only value for which I could get D(P)(8e-8)=0 was for P(0)=0. Now that I have the corrected coefficient, the problem is behaving differently. I'll keep looking. - Allan

## The problem...

As it is defined: A:=abs(3.0164*10^14*z-2.8071*10^7)*(-0.82500e-2*z+3.*10^(9))^2; alpha_z:=3.0164*10^14*z-2.8071*10^7; beta_z:=9.9000*10^15*z-5.400*10^8; gamma_z:=6.6e9; # (taken as constant for now) bcs := D(P)(0)=0,D(P)(8e-8)=0; deq := A*diff(P(z),z,z)-alpha_z*P(z)-beta_z*P(z)^3-gamma_z*P(z)^6=0; the solution (unless I am mistaken) is single-valued over the speccified interval. If you look at the coefficient of diff(P(z),z,z) with the above definitions, you get: > lc := abs(.3016400000e15*z-28071000.00)*(-.82500e-2*z+3000000000.)^2; 15 8 lc := | 0.3016400000 10 z - 0.2807100000 10 | 10 2 (-0.0082500 z + 0.3000000000 10 ) > fsolve(lc=0,z=0..8e-8); 15 8 fsolve(| 0.3016400000 10 z - 0.2807100000 10 | 10 2 -7 (-0.0082500 z + 0.3000000000 10 ) = 0, z, 0 .. 0.8 10 ) So as far as I can tell, there should only be one solution, and that solution should be zero. Did I make the typo correction incorrectly in the provided value for A? - Allan

## Older maple...

> Dear Allan! > > thanks for your post. I can imagine what amount of work that is, > and I am not sure if I want to endure this, and if I really could > do it. I know by now, that is essentially the java and the > dynamical loader which cause the problems and I was hoping for a > sort of a direct fix. Well Maple 6-7 had nothing to do with java at all, while Maple 8 only used java for Maplets, so really the problem is dynamical loading and the changes made on Linux to gcc. > I wonder what an 'unlimited' (in time) license does mean. Sure > enough, if one would have the sources then another make, make > install would do the trick. I see the legitimate interest and > need of maple soft to sell new products, but at some times I > would just be happy to be able to get an old version running > _painlessly_ It does run painlessly on an older O/S :-) Seriously, there's just no way to tell what will happen next on Linux. When the change to the behavior of glibc happened, I expect the developers worked to change usage of the relevant bits of glibc so that the binary would become independent of that change, but this could not be done until it was detected (which of course was after it happened) and corresponds to after Maple 8 was released. > Even maple 10 does not run out of the box on a decent linux > machine, like SuSE 10.1, and there is no Maple 11 around yet > which the should run (??) This is more glibc changes (sigh...) combined with binary incompatibilities wrt. those changes for java. Did you see the exceptionally good post about modifying the installer to comment out the LD_ASSUME_KERNEL in the installation script? > I think that I will possibly just run an old retired computer > without internet connection as a sort of old maple resort home > thing for testing programs across platforms. That's definitely the easier way (easier than the way I do it anyway...) - Allan

## Another Method...

Another method of determining the point at which the solution attains a certain value is to use stop conditions: > dEq := diff(theta(t), t, t) = sin(theta(t)): > soln := dsolve({dEq, theta(0) = 0, (D(theta))(0) = 1}, theta(t),numeric, > stop_cond=[theta(t)-2]); soln := proc(x_rkf45) ... end proc > soln(10); Warning, cannot evaluate the solution further right of 1.4849919, stop condition #1 violated [t = 1.48499193139580199, theta(t) = 2.00000000000000044, d -- theta(t) = 1.95763036755598319] dt - Allan

## Another Method...

Another method of determining the point at which the solution attains a certain value is to use stop conditions: > dEq := diff(theta(t), t, t) = sin(theta(t)): > soln := dsolve({dEq, theta(0) = 0, (D(theta))(0) = 1}, theta(t),numeric, > stop_cond=[theta(t)-2]); soln := proc(x_rkf45) ... end proc > soln(10); Warning, cannot evaluate the solution further right of 1.4849919, stop condition #1 violated [t = 1.48499193139580199, theta(t) = 2.00000000000000044, d -- theta(t) = 1.95763036755598319] dt - Allan

## Update...

This bug is now fixed in Maple 11 (now in beta) As an aside, I have noticed that in your worksheets you often use a *very* fine initial mesh, even when it is not necessary. The resulting solution is obtained quite slowly. Use of a coarser initial mesh will provide a solution with the same accuracy much more rapidly in the cases where a solution can be obtained. In the example you attached to this post you use 'initmesh=8*128', and the computation runs 5 seconds on my machine. Instead choosing 'initmesh=128' produces a solution with the same accuracy in 0.7 sec, and a solution 1000 times more accurate in 0.73 sec. For a smaller problem (as in this example), the difference is not huge, but for larger problems with many equations, the difference can be much more substantial. Sincerely, Allan Wittkopf

## Update...

This bug is now fixed in Maple 11 (now in beta) As an aside, I have noticed that in your worksheets you often use a *very* fine initial mesh, even when it is not necessary. The resulting solution is obtained quite slowly. Use of a coarser initial mesh will provide a solution with the same accuracy much more rapidly in the cases where a solution can be obtained. In the example you attached to this post you use 'initmesh=8*128', and the computation runs 5 seconds on my machine. Instead choosing 'initmesh=128' produces a solution with the same accuracy in 0.7 sec, and a solution 1000 times more accurate in 0.73 sec. For a smaller problem (as in this example), the difference is not huge, but for larger problems with many equations, the difference can be much more substantial. Sincerely, Allan Wittkopf

## Investigated...

OK, this is a bug. Given where it occurs in the code, it will only arise when: 1) Continuation is in use 2) The problem requires mesh adaptation to reduce the error 3) No other steps occur between the completion of the continuation and the start of the mesh adaptation. The reason that this works when you set tolerance to 1e-2 is because mesh adaptation is not required to reach that error tolerance. This also shows up on Linux with a tolerance of 1e-5 for the same problem (Linux gets one more Digit of accuracy than Windows). Strangely enough, even dodging the bug (in the code), you still do not get a successful computation, as the requested tolerance is too small for hardware Digits (the solution is O(1e6)). This will be fixed for Maple 11. - Allan

## Investigated...

OK, this is a bug. Given where it occurs in the code, it will only arise when: 1) Continuation is in use 2) The problem requires mesh adaptation to reduce the error 3) No other steps occur between the completion of the continuation and the start of the mesh adaptation. The reason that this works when you set tolerance to 1e-2 is because mesh adaptation is not required to reach that error tolerance. This also shows up on Linux with a tolerance of 1e-5 for the same problem (Linux gets one more Digit of accuracy than Windows). Strangely enough, even dodging the bug (in the code), you still do not get a successful computation, as the requested tolerance is too small for hardware Digits (the solution is O(1e6)). This will be fixed for Maple 11. - Allan

## The two worksheets...

The problem with the first worksheet is a bug, and I notice that it is still in Maple 10 if you keep the tolerance at 1e-4. I will investigate.

## The two worksheets...

The problem with the first worksheet is a bug, and I notice that it is still in Maple 10 if you keep the tolerance at 1e-4. I will investigate.
 1 2 3 4 5 6 Page 4 of 6
﻿