PatrickT

Dr. Patrick T

2153 Reputation

18 Badges

16 years, 96 days

MaplePrimes Activity


These are replies submitted by PatrickT

acer, actually I do intend to use events and eventstop to ensure the simulations don't stray too far off their intended path, so I will do just as you suggest:

"Your final code could do both: check the eventstop and trap errors using try..catch."

acer and preben,

it would be great if try/catch could also catch warnings, with maybe a specific syntax like catch-warnings and catch-errors or something.  In my earlier posts I copy-pasted this message:

Warning, cannot evaluate the solution further left of -.35923730e-2, probably a singularity

This is a warning message about a singularity.

Further down, I copy-pasted the following error message:

Error, (in unknown) cannot evaluate the solution further right of .99999999, probably a singularity

This is an error message about a singularity.

Presumably try/catch will catch the latter but not the former. I'm not sure now why I got different messages, a warning here but an error there, I would need to carefully retrace my steps (most likely there was a call to dsolve here but a call to odeplot there, but I'm not certain), but there isn't a great urgency to find out the reasons in the light of your comment Preben: it is clear that various Maple routines are not in sync when it comes to throwing up warnings or errors. The way I have it set up now it throws up errors and gets caught by try/catch, which is what I want.

When I read the help file I took the word "exception" to mean something like "warning" but upon rereading I see that it might as well be a synonym for "error". Not sure if there are any differences between an exception and an error...

Preben, acer, thanks a lot for your follow-up explanations. The lesson, then, is that "what comes after the colon in catch: does not control what gets caught."  I hadn't paid attention to the colon, to be honest, and didn't realize its importance. The difference between the two cases below is the position of the colon, but I mistakenly took it to mean that the rest of the statement mattered. Thanks for clarifying this, it no doubt would have caught up with me soon.  ;-)

 

Case 1. Works. But the line: StringTools[FormatMessage]("probably a singularity"): is irrelevant.

restart;
ReInit := proc(dsol,tf::realcons)
    local imax, eps, i, pts, ini ;
    imax:=100;
    eps:=.1;
    for i to imax do
    try
        dsol(tf);
        return dsol;
    catch:StringTools[FormatMessage]("probably a singularity"): # Case 1
        pts := Statistics:-Sample(Uniform(-eps,eps),1) ;
        ini := [ rhs(dsol(0)[2])+pts[1]] ;
        dsol('initial'=ini);
        print(i,ini)      
    end try
    end do;
    if i=imax+1 then FAIL end if;     
    end proc :
sol := dsolve({diff(u(t),t)=u(t)^2, u(0)=1},numeric,'output' = listprocedure);
   sol := [t = proc(t)  ...  end;, u(t) = proc(t)  ...  end;]
sol(1):
ReInit(sol,1.0001);
Error, (in unknown) cannot evaluate the solution further right of .99999999, probably a singularity
                     1, [1.06294473727864]
                     2, [1.14410312469376]
                     3, [1.06950048795246]
                              FAIL

Case 2. Fails. There is no colon immediately after catch.

restart;
ReInit := proc(dsol,tf::realcons)
    local imax, eps, i, pts, ini ;
    imax:=100;
    eps:=.1;
    for i to imax do
    try
        dsol(tf);
        return dsol;
    catch "probably a singularity": # Case 2
        pts := Statistics:-Sample(Uniform(-eps,eps),1) ;
        ini := [ rhs(dsol(0)[2])+pts[1]] ;
        dsol('initial'=ini);
        print(i,ini)      
    end try
    end do;
    if i=imax+1 then FAIL end if;     
    end proc :
sol := dsolve({diff(u(t),t)=u(t)^2, u(0)=1},numeric,'output' = listprocedure);
   sol := [t = proc(t)  ...  end;, u(t) = proc(t)  ...  end;]
sol(1):
ReInit(sol,1.0001);
Error, (in unknown) cannot evaluate the solution further right of .99999999, probably a singularity
Error, (in unknown) cannot evaluate the solution further right of .99999999, probably a singularity

Preben, acer, thanks a lot for your follow-up explanations. The lesson, then, is that "what comes after the colon in catch: does not control what gets caught."  I hadn't paid attention to the colon, to be honest, and didn't realize its importance. The difference between the two cases below is the position of the colon, but I mistakenly took it to mean that the rest of the statement mattered. Thanks for clarifying this, it no doubt would have caught up with me soon.  ;-)

 

Case 1. Works. But the line: StringTools[FormatMessage]("probably a singularity"): is irrelevant.

restart;
ReInit := proc(dsol,tf::realcons)
    local imax, eps, i, pts, ini ;
    imax:=100;
    eps:=.1;
    for i to imax do
    try
        dsol(tf);
        return dsol;
    catch:StringTools[FormatMessage]("probably a singularity"): # Case 1
        pts := Statistics:-Sample(Uniform(-eps,eps),1) ;
        ini := [ rhs(dsol(0)[2])+pts[1]] ;
        dsol('initial'=ini);
        print(i,ini)      
    end try
    end do;
    if i=imax+1 then FAIL end if;     
    end proc :
sol := dsolve({diff(u(t),t)=u(t)^2, u(0)=1},numeric,'output' = listprocedure);
   sol := [t = proc(t)  ...  end;, u(t) = proc(t)  ...  end;]
sol(1):
ReInit(sol,1.0001);
Error, (in unknown) cannot evaluate the solution further right of .99999999, probably a singularity
                     1, [1.06294473727864]
                     2, [1.14410312469376]
                     3, [1.06950048795246]
                              FAIL

Case 2. Fails. There is no colon immediately after catch.

restart;
ReInit := proc(dsol,tf::realcons)
    local imax, eps, i, pts, ini ;
    imax:=100;
    eps:=.1;
    for i to imax do
    try
        dsol(tf);
        return dsol;
    catch "probably a singularity": # Case 2
        pts := Statistics:-Sample(Uniform(-eps,eps),1) ;
        ini := [ rhs(dsol(0)[2])+pts[1]] ;
        dsol('initial'=ini);
        print(i,ini)      
    end try
    end do;
    if i=imax+1 then FAIL end if;     
    end proc :
sol := dsolve({diff(u(t),t)=u(t)^2, u(0)=1},numeric,'output' = listprocedure);
   sol := [t = proc(t)  ...  end;, u(t) = proc(t)  ...  end;]
sol(1):
ReInit(sol,1.0001);
Error, (in unknown) cannot evaluate the solution further right of .99999999, probably a singularity
Error, (in unknown) cannot evaluate the solution further right of .99999999, probably a singularity

Thanks pagan, things are clearer now, to put it in simple terms, if I understand correctly, each time I call dsol with a new set of initial conditions it "overwrites" the previous initial conditions. My confusion comes from the fact that I write dsol on the left-hand side of an assignment and dsolve(system,initial_conditions,etc) on the right-hand side, and so I expected that dsol would forever be associated with the initial_conditions used to define it, but it doesn't work that way. Knowing that, I have rewritten my procedure to take as separate arguments dsol and the initial conditions and each time being explicit about which initial conditions are used with dsol. I hope this is the way forward. Thanks for your help in clarifying these issues.

Thanks Preben. As far as I can tell,   

catch:StringTools:-FormatMessage("probably a singularity"):

or

catch: StringTools:-FormatMessage(lastexception[2..-1]);

have the same effect in my specific application. But the latter will be more general and will catch other types of errors (other than singularity errors) if such happen. So it is a more useful syntax for my purpose. Thanks.

 

P.S. For some reason the loop is not working today, it's only catching once and redefining the initial conditions once. Strange. I could swear it was working last night. I must investigate this ...

P.S.2 : I had misplaced a line in the code pasted above (now fixed) and this had caused the loop to break. It's now fixed. It works exactly as desired. It's pretty fast. No complaints. Thanks again.

Thanks pagan, things are clearer now, to put it in simple terms, if I understand correctly, each time I call dsol with a new set of initial conditions it "overwrites" the previous initial conditions. My confusion comes from the fact that I write dsol on the left-hand side of an assignment and dsolve(system,initial_conditions,etc) on the right-hand side, and so I expected that dsol would forever be associated with the initial_conditions used to define it, but it doesn't work that way. Knowing that, I have rewritten my procedure to take as separate arguments dsol and the initial conditions and each time being explicit about which initial conditions are used with dsol. I hope this is the way forward. Thanks for your help in clarifying these issues.

Thanks Preben. As far as I can tell,   

catch:StringTools:-FormatMessage("probably a singularity"):

or

catch: StringTools:-FormatMessage(lastexception[2..-1]);

have the same effect in my specific application. But the latter will be more general and will catch other types of errors (other than singularity errors) if such happen. So it is a more useful syntax for my purpose. Thanks.

 

P.S. For some reason the loop is not working today, it's only catching once and redefining the initial conditions once. Strange. I could swear it was working last night. I must investigate this ...

P.S.2 : I had misplaced a line in the code pasted above (now fixed) and this had caused the loop to break. It's now fixed. It works exactly as desired. It's pretty fast. No complaints. Thanks again.

In trying to move forward, I took a slightly different approach. The code is a little longer and it's a little less convenient, but it's working. I proceeded in two steps. First, I set up a procedure to compute new initial conditions. Secondly, I use the new initial conditions with dsolve. Thus, I'm computing the solution twice (not counting the intermediate computations), once to tweak the initial conditions and once to compute the trajectories. I wouldn't be surprised if there were a simpler, more efficient way.

Some  further tweaking will be needed to get convergence at all points, I now have about 32 out of 40 points that converge, an improvement from the initial 26...

### Procedure MakeInits searches for initial conditions that converge

MakeInits := proc( dsol, {ini::list:=[0,0,0]}, {endpoint::realcons:=0}, {verbose::boolean:=false} )
    local i, j, n, e, sol, pts, ini0, ini1 ;
    n := 200 ;
    e := 1e-5 ;
    ini0 := map(rhs,ini) ;
    for i to n do
    try
        dsol('initial'=ini0) ;
        dsol(endpoint) ;
        if verbose=true then print("tried") ; end if ;
        return ini0 ;
    catch:StringTools[FormatMessage]("probably a singularity"):
        pts := Statistics:-Sample(Normal(0,e),3) ;
        # pts := Statistics:-Sample(Uniform(-e,e),3) ;
        ini1 := [ ini0[1]+pts[1]
                   , ini0[2]+pts[2]
                   , ini0[2]+pts[3]
                   ] ;
        dsol('initial'=ini1) ;
        if verbose=true then print(i,ini1) ; end if ;
        if verbose=true then print("caught") ; end if ;
    end try
    end do ;
    if i=n+1 then FAIL ;  end if ;
    return ini1 ;
    end proc :

### Create New Initial Conditions

Inits := [ seq( MakeInits(
    'dsol' = Sol(candidates[i]), 'ini' = ini[i], 'endpoint' = -10, 'verbose' = true
    ) , i = 1 .. 40 ) ] ;

### Solve System With New Initial Conditions

etc.

EDIT : Fixed the misplaced line "return ini1;" in the code above. Removed the line "sol:=dsol" since pagan explained (below) that it doesn't help.

In trying to move forward, I took a slightly different approach. The code is a little longer and it's a little less convenient, but it's working. I proceeded in two steps. First, I set up a procedure to compute new initial conditions. Secondly, I use the new initial conditions with dsolve. Thus, I'm computing the solution twice (not counting the intermediate computations), once to tweak the initial conditions and once to compute the trajectories. I wouldn't be surprised if there were a simpler, more efficient way.

Some  further tweaking will be needed to get convergence at all points, I now have about 32 out of 40 points that converge, an improvement from the initial 26...

### Procedure MakeInits searches for initial conditions that converge

MakeInits := proc( dsol, {ini::list:=[0,0,0]}, {endpoint::realcons:=0}, {verbose::boolean:=false} )
    local i, j, n, e, sol, pts, ini0, ini1 ;
    n := 200 ;
    e := 1e-5 ;
    ini0 := map(rhs,ini) ;
    for i to n do
    try
        dsol('initial'=ini0) ;
        dsol(endpoint) ;
        if verbose=true then print("tried") ; end if ;
        return ini0 ;
    catch:StringTools[FormatMessage]("probably a singularity"):
        pts := Statistics:-Sample(Normal(0,e),3) ;
        # pts := Statistics:-Sample(Uniform(-e,e),3) ;
        ini1 := [ ini0[1]+pts[1]
                   , ini0[2]+pts[2]
                   , ini0[2]+pts[3]
                   ] ;
        dsol('initial'=ini1) ;
        if verbose=true then print(i,ini1) ; end if ;
        if verbose=true then print("caught") ; end if ;
    end try
    end do ;
    if i=n+1 then FAIL ;  end if ;
    return ini1 ;
    end proc :

### Create New Initial Conditions

Inits := [ seq( MakeInits(
    'dsol' = Sol(candidates[i]), 'ini' = ini[i], 'endpoint' = -10, 'verbose' = true
    ) , i = 1 .. 40 ) ] ;

### Solve System With New Initial Conditions

etc.

EDIT : Fixed the misplaced line "return ini1;" in the code above. Removed the line "sol:=dsol" since pagan explained (below) that it doesn't help.

If I understand the original question, I'm having basically the same problem.

http://www.mapleprimes.com/questions/134162-Procedure-To-Modify-Initial-Conditions-In-Dsolve#comment134231&rid=437061

If I understand the original question, I'm having basically the same problem.

http://www.mapleprimes.com/questions/134162-Procedure-To-Modify-Initial-Conditions-In-Dsolve#comment134231&rid=437061

@Christopher2222 

Earlier in the thread, Will said "Votes on comments will be different than votes on Answers, in that they won't cause the comments to be re-ordered."

Personally I don't find the reordering of answers useful at all. The idea probably comes from sites like stackoverflow, where it makes sense, but I don't think it serves a purpose here. In my humble opinion, the number of votes is a good enough signal to identify the useful content. Being able to follow a conversation in the correct order is essential.

I do have a follow-up question. It's probably a simple question, but it comes up in a fairly messy setup, so let me try to explain the problem in simple terms. If the problem statement is unclear, upon request I will attach a worksheet in which the problem arises.

I "dsolve" several initial conditions at once, along the lines of:

  seq( plots:-odeplot( Sol(i)
       , [ u(t), v(t), w(t) ]
       , t = -100 ..0
       ) , i = 1..40
     ) :

and obtain warning messages like these:


Warning, cannot evaluate the solution further left of -38.078909, event #4 triggered a halt
...
Warning, cannot evaluate the solution further left of -51.084118, event #4 triggered a halt
Warning, cannot evaluate the solution further left of -.35923730e-2, probably a singularity
...
Warning, cannot evaluate the solution further left of -.19887345e-1, probably a singularity


The first few trajectories do converge but the last few don't.

Let's focus on what happens to trajectories 26 and 27 (of a total of 40), where 26 converges but 27 does not.

Take trajectory 26, which converges, and apply the procedure ReInit, this then simply returns the procedure unchanged. The output looks like this:

sol26 := Sol(26):
sol26r := ReInit(sol26,-10):
sol26r(-100);
Warning, cannot evaluate the solution further left of -51.084118, event #2 triggered a halt
Warning, cannot evaluate the solution further left of -51.084118, event #2 triggered a halt
Warning, cannot evaluate the solution further left of -51.084118, event #2 triggered a halt
Warning, cannot evaluate the solution further left of -51.084118, event #2 triggered a halt
         [t(-100) = HFloat(-51.084118951402644),
           u(t)(-100) = HFloat(-7.333333333999999),
           r(t)(-100) = HFloat(-0.033719542481319885),
           q(t)(-100) = HFloat(0.07082105340966839)]

Take trajectory 27, which does not converge with the original initial conditions, and apply the procedure ReInit to search for initial conditions until convergence obtains (in this example, it worked). The output looks like this:

sol27 := Sol(27):
sol27r := ReInit(sol27,-10):
sol27r(-100);
Warning, cannot evaluate the solution further left of -53.186755, event #2 triggered a halt
Warning, cannot evaluate the solution further left of -53.186755, event #2 triggered a halt
Warning, cannot evaluate the solution further left of -53.186755, event #2 triggered a halt
Warning, cannot evaluate the solution further left of -53.186755, event #2 triggered a halt
         [t(-100) = HFloat(-53.186755254308736),
           u(t)(-100) = HFloat(-7.333333333999995),
           r(t)(-100) = HFloat(-0.033719585093739166),
           q(t)(-100) = HFloat(0.07082134219745041)]

In the above, the procedure succeeded in either returning the solution untouched or finding a new list of initial conditions that do converge. This is great news.

The problem I have, however, is that the new set of initial conditions "contaminates" subsequent computations, in the sense that the statement dsol('initial'=ini) inside the procedure clears the original initial conditions of Sol(i) even when the i index is incremented. In other words, after ReInit has determined a new list of initial conditions for simulations number 27, it then applies it to subsequent simulation, 28, 29, etc..

To illustrate the problem, if I odeplot the solution, only the last trajectory is plotted, with plots:-odeplot(sol26r) completely oblivious of the previous computation supposedly stored in sol26r (it is as if sol27r overwrites sol26r).

Thus, (I think) I'm looking for a way to not alter my original call. Or maybe I need to store the original initial conditions and re-initialize again?

I tried to have sol:=dsol: at the beginning of the procedure and dsol:=sol: at the end of it, so that the new set of initial conditions would be applied to the copy and not to the original, but that didn't work (why?).

As I said, if my question is obscure, I can upload a worksheet. I can also branch this off if it drifts off topic.

Thanks a lot for suggestions!

I do have a follow-up question. It's probably a simple question, but it comes up in a fairly messy setup, so let me try to explain the problem in simple terms. If the problem statement is unclear, upon request I will attach a worksheet in which the problem arises.

I "dsolve" several initial conditions at once, along the lines of:

  seq( plots:-odeplot( Sol(i)
       , [ u(t), v(t), w(t) ]
       , t = -100 ..0
       ) , i = 1..40
     ) :

and obtain warning messages like these:


Warning, cannot evaluate the solution further left of -38.078909, event #4 triggered a halt
...
Warning, cannot evaluate the solution further left of -51.084118, event #4 triggered a halt
Warning, cannot evaluate the solution further left of -.35923730e-2, probably a singularity
...
Warning, cannot evaluate the solution further left of -.19887345e-1, probably a singularity


The first few trajectories do converge but the last few don't.

Let's focus on what happens to trajectories 26 and 27 (of a total of 40), where 26 converges but 27 does not.

Take trajectory 26, which converges, and apply the procedure ReInit, this then simply returns the procedure unchanged. The output looks like this:

sol26 := Sol(26):
sol26r := ReInit(sol26,-10):
sol26r(-100);
Warning, cannot evaluate the solution further left of -51.084118, event #2 triggered a halt
Warning, cannot evaluate the solution further left of -51.084118, event #2 triggered a halt
Warning, cannot evaluate the solution further left of -51.084118, event #2 triggered a halt
Warning, cannot evaluate the solution further left of -51.084118, event #2 triggered a halt
         [t(-100) = HFloat(-51.084118951402644),
           u(t)(-100) = HFloat(-7.333333333999999),
           r(t)(-100) = HFloat(-0.033719542481319885),
           q(t)(-100) = HFloat(0.07082105340966839)]

Take trajectory 27, which does not converge with the original initial conditions, and apply the procedure ReInit to search for initial conditions until convergence obtains (in this example, it worked). The output looks like this:

sol27 := Sol(27):
sol27r := ReInit(sol27,-10):
sol27r(-100);
Warning, cannot evaluate the solution further left of -53.186755, event #2 triggered a halt
Warning, cannot evaluate the solution further left of -53.186755, event #2 triggered a halt
Warning, cannot evaluate the solution further left of -53.186755, event #2 triggered a halt
Warning, cannot evaluate the solution further left of -53.186755, event #2 triggered a halt
         [t(-100) = HFloat(-53.186755254308736),
           u(t)(-100) = HFloat(-7.333333333999995),
           r(t)(-100) = HFloat(-0.033719585093739166),
           q(t)(-100) = HFloat(0.07082134219745041)]

In the above, the procedure succeeded in either returning the solution untouched or finding a new list of initial conditions that do converge. This is great news.

The problem I have, however, is that the new set of initial conditions "contaminates" subsequent computations, in the sense that the statement dsol('initial'=ini) inside the procedure clears the original initial conditions of Sol(i) even when the i index is incremented. In other words, after ReInit has determined a new list of initial conditions for simulations number 27, it then applies it to subsequent simulation, 28, 29, etc..

To illustrate the problem, if I odeplot the solution, only the last trajectory is plotted, with plots:-odeplot(sol26r) completely oblivious of the previous computation supposedly stored in sol26r (it is as if sol27r overwrites sol26r).

Thus, (I think) I'm looking for a way to not alter my original call. Or maybe I need to store the original initial conditions and re-initialize again?

I tried to have sol:=dsol: at the beginning of the procedure and dsol:=sol: at the end of it, so that the new set of initial conditions would be applied to the copy and not to the original, but that didn't work (why?).

As I said, if my question is obscure, I can upload a worksheet. I can also branch this off if it drifts off topic.

Thanks a lot for suggestions!


follow-up on the above. I got the solution from the help pages, which state:

"You can use StringTools[FormatMessage] to format the string produced from the lastexception value."

http://www.maplesoft.com/support/help/Maple/view.aspx?path=error


follow-up on the above. I got the solution from the help pages, which state:

"You can use StringTools[FormatMessage] to format the string produced from the lastexception value."

http://www.maplesoft.com/support/help/Maple/view.aspx?path=error

Dear Preben and others interested in catching, the following works:

catch:StringTools[FormatMessage]("probably a singularity"):

P.S. Only did preliminary tests, wanted to write back asap, if further problems arise I will write again.

Thanks again Preben!

First 16 17 18 19 20 21 22 Last Page 18 of 93