Allan Wittkopf

240 Reputation

7 Badges

18 years, 255 days

MaplePrimes Activity


These are answers submitted by Allan Wittkopf

Could you possibly attach a worksheet or text form of the example that caused this error? Thanks, Allan
Unfortunately Maple does not currently support Elliptic PDE numerical solutions. It is currently restricted to finite interval hyperbolc and parabolic PDE problems. If this is not what you are looking for, then please ignore this message. - Allan
After a little work with the problem it appears to me as though the solution P(x)=0 is unique for the stated conditions. The shooting problem - the ODE set up as an IVP with D(P)(0)=0 and P(0) varying, produces a smooth monotonic distribution of values of D(P)(8e-8) that is only zero when P(0)=0 and hence P(z)=0. I just uploaded a worksheet that demonstrates this Download 233_BVP_shoot.mws
View file details. Could you provide some information on the MatLab solution? Does the solution appear in the same way regarless of the number of points used in the domain, or is it specific to a fixed number of points (say 100)? Do you have a gif of a plot of the solution? Do you have a value of P(0) from the MatLab result? Thanks, Allan
Can I assume that you meant to write: (D@@2)(x)(t)=(-g/v)*sqrt(D(x)(t)^2+D(y)(t)^2)*D(x)(t) (D@@2)(y)(t)=g-g/v*sqrt(D(x)(t)^2+D(y)(t)^2)*D(y)(t) where I changed (D@@2)(x)(t) to (D@@2)(y)(t) in the second equation? Also, what values are you using for 'g' and 'v'? Thanks, Allan
The problem is actually a bug in the initial solution profile generation. This bug only occurs when the same variable occurs in multiple mixed boundary conditions, as it does in your example. A simple workaround for this, at least for this problem, is to solve the mixed boundary conditions for the left endpoint variables via the command: cc := op(solve({cc},{ww[1,2,0](0),ww[1,2,1](0),ww[1,2,2](0),ww[1,2,3](0)})); prior to the call to dsolve numeric. Having done this, dsolve numeric returns a computed solution. This bug will be fixed.
Can you try renaming the libc.so.6 to libc.so.6.old and try this again? Sincerely, Allan Wittkopf
If you didn't create one, then you likely have nothing to worry about. By default on an install, the ini file is not installed. If you want to check, then look in your Users directory (I believe it is a subdirectory of the Maple installation directory) for a file called maple.ini. If you see a file called maple?.ini (where ? is a version number), then just ignore that, as is contains version-specific initialization information. I have not used windows in so long (Linux!), I am not certain about the name of the directory, but as I recall it was something like "User" or "Users". So you were able to run both worksheets successfully in Maple 10? You mentioned that you had Maple 9.5, but was that where you were seeing the problems? I was able to run both problems successfully in Maple 9.5, but I was running on Linux. Windows behaves a bit differently, so I'll have to get access to a windows machine, and try this out if this is the case.
As for the second example, I was able to repeat the problem you show in Maple 9, and I even recall the cause of the issue (there was a mismatch between the number of mesh points that could be set internally by the routine, and the maximum number of mesh points allowed for use), and this was fixed for Maple 9.5. In Maple 9.5 and Maple 10, the problem completes with the solution as seen in Download 233_2847_problem-2.mws
View file details. Am I correct in guessing that you are using a version of Maple that is older than Maple 9.5? If so, I strongly recommend upgrading, as there were several improvements made since 9.0 with respect to various issues with numerical bvp solution, validation that the desired accuracy goal was reached across the whole solution, better handling of continuation problems (continuation parameter), efficiency improvements, etc. Sincerely, Allan Wittkopf
Regarding the first issue, I expect that this should not happen. Without the DE system available, I cannot tell you why it occurs, or how to avoid the problem. Could you please post the system, either in worksheet or text format? As for the second, the BVP solver effectively uses a mesh doubling technique, so you may get better results by setting initmesh to 32*128. If this is the same system, then I reiterate, please post the system. Note that using this high a number of initial mesh points will make the process quite slow, and you only want to use this high a number if your problem is quite badly behaved (you can set maxmesh high though). Sincerely, Allan Wittkopf
Numerical integration, along with several other numerical features in Maple, attempt to perform the integration in what is called 'hardware floats' when the Digits environment variable is set quite low. Consider:
> xpr := -20.01782577*I*exp(-490.0000000-50.*x)*(-0.1625801243*10^212
> *I+erf(0.2258769757*I*(9801.+980.*x)^(1/2)-0.2258769757*I));
                                                                 212
xpr := -20.01782577 I exp(-490.0000000 - 50. x) (-0.1625801243 10    I

                                          1/2
     + erf(0.2258769757 I (9801. + 980. x)    - 0.2258769757 I))

> f := unapply(xpr,x);
f := x -> -20.01782577 I exp(-490.0000000 - 50. x) (

                    212
    -0.1625801243 10    I

                                          1/2
     + erf(0.2258769757 I (9801. + 980. x)    - 0.2258769757 I))

> evalhf(f(10));
                          Float(undefined) I

> plot(Re(xpr),x=0..30);

0.35                                                                  
  +A                                                                  
  +AA                                                                 
0.3AA                                                                 
  +AA                                                                 
  +AA                                                                 
0.25 A                                                                
  +A A                                                                
  +A A                                                                
0.2A  A                                                               
  *A  A                                                               
  *    A                                                              
0.15   A                                                              
  *     A                                                             
0.1     AA                                                            
  *      AA                                                           
  *       AA                                                          
0.05       A                                                          
  *                                                                   
  *                                                                   
  *-+-+-+--+-+-+-+-+--+-+-+-+--+-+-+-+-+--+-+-+-+--+-+-+-+-+--+-+-+-+ 
0            5         10         15         20         25         30
So notice that the evaluation at 10 is undefined, and the plot simply stops just before 5. This is because the individual quantities inside the equation become so large that they cannot even be represented in a hardware floats (the limitation is around 1e310). There is a solution for this, and that is to force the use of additional Digits to avoid use of hardware floats, and to specify the desired tolerance independently as follows:
> Digits := 20:
> evalf(Int(xpr,x=0..30,epsilon=1e-9));
                        0.99995862345984971643
So in this case we are working with 20 Digit floats, and requesting an error bound of 1e-9 on the solution. - Allan
You can obtain a plot of the transform numerically, for some fixed value of 'k', as follows: Since JacobiSN is O(1), you can choose to terminate the integration early and still get a reasonable precision, so if you want a precision of, say 1e-4 in the plot, you could do:
> trf := s->evalf(Int(exp(-s*t)*JacobiSN(t,1/3),t=0..-ln(1e-4)/s)):
> plot(trf,1/2..5);
  *                                                                           
0.8A                                                                          
  + A                                                                         
  +  A                                                                        
  +  AA                                                                       
  +    A                                                                      
0.6    AA                                                                     
  +      A                                                                    
  +       A                                                                   
  +        AA                                                                 
  +         AA                                                                
0.4           AA                                                              
  +            AAA                                                            
  +              AAA                                                          
  +                 AAA                                                       
  +                   AAAA                                                    
0.2                       AAAA                                                
  +                          AAAAAA                                           
  +                                AAAAAAAAA                                  
  +                                         AAAAAAAAAAAAAAAA                  
  +-+--+--+--+---+--+--+---+--+--+---+--+--+--+---+--+--+--****************** 
          1                2               3                4               5 
> trf := s->evalf(Int(exp(-s*t)*JacobiSN(t,9/10),t=0..-ln(1e-4)/s)):
> plot(trf,1/2..5);                                     
  *                                                                           
  *A                                                                          
1 +A                                                                          
  + A                                                                         
  + AA                                                                        
0.8  A                                                                        
  +   A                                                                       
  +    A                                                                      
  +     A                                                                     
0.6      A                                                                    
  +       A                                                                   
  +        AA                                                                 
  +         AA                                                                
0.4           AA                                                              
  +             AA                                                            
  +               AAA                                                         
  +                 AAAAA                                                     
0.2                     AAAAAA                                                
  +                           AAAAAAAA                                        
  +                                   AAAAAAAAAAAAAAAAA                       
  +-+--+--+--+---+--+--+---+--+--+---+--+--+--+---+--+-********************** 
          1                2               3                4               5 

There are several problems with the input form, and I will describe these in turn: 1) The error you receive should suggest reading the help page, from which you can discover that multiple variables need to be entered in a list or set, so change the call line to: dsolve({RF,ic},{x[1](t),x[2](t)},type=numeric); 2) Even with this, it will not work, as the initial conditions are over-specified. The equations are first order, and as such you cannot dictate the values of the derivatives at the initial point, so the initial conditions should be changed to: ic:=x[1](0)=5, x[2](0)=1; 3) Even with this, it will not work, but the remaining problem is obvious, and you will get the error message: Error, (in dsolve/numeric) x[1](t) and x[1] cannot both appear in the given ODE Note: x[1](t) is a function of 't', while x[1] is a constant. I assume when you wrote x[1] you meant x[1](t), so replace the standalone x[1] and x[2] with x[1](t) and x[2](t), and then: ... > dsolve({RF,ic},{x[1](t),x[2](t)},type=numeric); proc(x_rkf45) ... end proc a solution.
The problem you have presented is an implicit complex DAE boundary value problem, but Maple is currently restricted to explicit real valued boundary value problems, so you need to do some manual work to get the input into a form that can be solved in Maple. I have uploaded a workup in Maple using some manipulation and techniques described in ?numeric_bvp,advanced. Note that the posed problem is *hard*. Note that it looks as though the boundary conditions for beta were incorrect. View 233_BVPWorkup.mw on MapleNET or Download 233_BVPWorkup.mw
If you look at ?dsolve,numeric,IVP you will see (procopts) that there is a capability for specifying a procedure that can be used to evaluate the derivative (instead of specifying the differential system as a set of equations). You can use this capability to construct a procedure that does root finding, computes an integral, etc. the possibilities are endless. Then when dsolve,numeric performs the integration, it calls the procedure to compute the derivatives, the root finding can then occur, and the resulting derivative values are returned.
Given that you can transform the PDE problem into an ODE problem, that seems to indicate that the problem is either parabolic or hyperbolic, which means that the error message may be false. Can you provide the system being used for the PDE solve, and the system being used for the ODE solve? We may be better able to help.
1 2 3 Page 2 of 3