tomleslie

13876 Reputation

20 Badges

15 years, 169 days

MaplePrimes Activity


These are replies submitted by tomleslie

When I look at a problem like this, I tend to start by deleting "options" that I wouldn't generally use. The two which stood out for me were Digits=18, and initmesh=30000.

The help page for dsolve numeric/BVP states that Maple's default value for initmesh is 8 and any user-supplied value should be between 8 and 8192. I have no idea why you were using 30000 - and so I deleted it. If I were really conceerned about accuracy in this problem, then after solving it on the defaults, I might try initmesh=16, 32 or 64, just to see if it made much difference - but 30000, come on! Why did you set this value for this parameter - honestly, what were you trying to achieve??

Digits=18 also seemed a little pointless - the accuracy of any dsolve/numeric calculation depends primarily on the settings of abserr and relerr within the dsolve() command - unless you change these I doubt if setting Digits=15 does very much for you - it might provide 15-digit output, but siince the default value of abserr is 1e-6, what are you getting from all of these digits except numerical noise?

(I also suspect that decreasing abserr and relerr would improve accuracy more than playing with initmesh (up to a point) so again, why play with initmes??)

I'm a big fan of the "Keep it simple" approach - you can translate this as "leave as much as possible on defaults". Only when thta fails should one start playing with available options, and as a general rule, move these "slowly" and "one at a time"

I'm slightly surprised that intsolve() is doing anything! Some observations

  1. Are the names 's' and 'r' in the integrand meant to be related to the previously defined functions s() and r()?
  2. if they are supposed to refer to the previously defined functions then they would have to be written as functions of either the integration variable, or the limit variable - so any of s(x), s(y), r(x), r(y) would make some kind of sense in the integrand.
  3. However if they are functions of the limit variable then they can be moved outside the integral, which doesn't make a *lot* of sense
  4. On the other hand if they are functions of the integration variable then you should be aware that they are undefined for y=0 (ie within the range of integration), because they contain terms which will evaluate to 0/0 - not good!
  5. Is there a missing multiplication sign immediately to the left of u() in the integrand. Since you are using 2D input, and post pictures (rather then code), it is difficult to be certain of this
  6. Do you really want to multiply the integral term by 1 - if so why?
  7. In the definition of the function s() there may be a missing multiplication immediately to the right of 0.0647
  8. Should one assume that you are interested in solutions of myeq=0 - surely solved by u(x)=0 for any s(), r(), so makes no sense! Are we missinf the rhs of myeq??

Maple tells me that this system cannot be solved when I force the use of the RKF method. On the other hand, if I jlet Maple pick its own solution methods, I get the result table

   S    M   Sq   Nt   Nb   Le   Pr    D(theta)(1)   D(phi)(1)
 0.0  0.5  1.0  0.2  0.5  1.0  1.0    1.00563,    1.34606
 0.5  0.5  1.0  0.2  0.5  1.0  1.0    0.96157,    1.29726
 1.0  0.5  1.0  0.2  0.5  1.0  1.0    0.91916,    1.24918
 0.0  0.0  1.0  0.2  0.5  1.0  1.0    1.00553,    1.34595

None of which agree with the solutions you provided. So either

  1. I can't type in the equations properly - although since you are too lazy to do this, I don't feel inclined to spend time checking
  2. I can't type in the boundary conditions correctly - although since you are too lazy to do this, I don't feel inclined to spend time checking
  3. The equations you supply don't correspond with the results you supply, no way I can check this
  4. Problems with one or more of the above may explain why Maple thinks that the Runge-Kutta method cannot be applied to this problem??
  5. Looks like you will have to do some work after all :-(

I would probably avoid creating then saving an animation: rather, create individual plots and export them to individual files using the exportplot() command, probably in jpg format.

Having created all the plots, assuming I know which ones I need for the animation, I'd use the

plots[surfdata](image=filename)   - see the last example on the surfdata help page

to retrieve individual plots and store them as elements in a list.

Then use display() with the 'insequence' option to display this list of plots as an animation

 

Maple's garbage collection is pretty efficient - you can force a garbage collection with gc(), but as the help page says - "Note, as of Maple 16 the use of gc is discouraged."

Need a numeric value for eta (unless it is the independent variable?)

I have been using MATLAB to examine the file which you uploaded

It contains 4 low-level `structs` of key/value pairs, combined into a cell array, which is then combined with a string into a higher level 'struc't, so if I want to examine any of the low-level entries say r_hyp, then after

A=load(('zdata.mat'), I have to use access commands such as

A.database_1{1}.r_hyp

A.database_1{2}.r_hyp

A.database_1{3}.r_hyp

To make matters worse the low-level structs have a real mixture of key:value pairs - some of the value entries are scalar floats, some are strings, and some are quite large matrices.

I really can't see any way that Maple will be able to import a data structure this complicated using simple Import() or ImportMatrix() commands with a 'mat' option

I also can't see any way of using low-level (fscanf?) commands because the file format is 'mat' - which is binary

What you esssentiall have is a database: and Maple's 'DataBase' package can be used to extract data from external databases using SQL (more-or-less). So one option would be if Matlab could export data iin a format whcih support JDBC. I believe that the matlab DATABASE toolbox has the functionality to create such a database: but I can't check this because I only have a bare matlab with no toolboxes:-(

From this point I think you have four options

  1. Create the original data in Maple (if possible) - ie bypass Matlab altogether
  2. From Matlab write out the individual structure fields into a set of files, and then read these individually in Maple - this would definitely work, but sounds like a fair amount of work
  3. Given a database toolbox for matlab, then write out the data in a format which supports JDBC connectivity, then use Maple's database() package to access this external database
  4. Give up ;-)

If you upload your matlabfile (*.mat) and the Maple code your are using to import it, I might be able to help: without either, no chance!

Don't even think about executing the attached worksheet until you have read th following very carefully

  1. I should have spotted the missing minus sign in eta definition - I apologise, now fixed
  2. When you said earlier that "J should be jump about 2 time rather than 1/4 times in your code." - I finally tracked this down. I had a premultiplying factor of 1/2*Pi rather than 1/(2*Pi) - I have now fixed this and the scale factor on J now seems to be OK. My screw-up, I apologise again
  3. Debugging the latest problem "but it doesn't show more resualts after 3 period." was more interesting. The first thing that took me a while to find was that you have changed the definition of omega from  omega:=x-> 1+eps*x, to omega:=1+sin( eps*x ). When I have to check code line by line, you have no idea how long it takes to find something like this. It would be so much easier if you listed the changes you make to my code!!!! Just remember next time you change something without telling me, I just might not bother to try and find it
  4. The "important" issue - why isn't any data generated beyong a fixed time-value - Basically Maple is trying to protect the user9 (you/me) from inadvertently(?) asking for *huge* numbers of function evaluations in the solutions of ODEs - this limit can be removed by setting the option maxfun=0 in the dsolve() calculation, which I have done.
  5. I have also updated the limits on required plots from 100 to 8*Pi/epsilon. Since the maxfun limit has been disabled, these will now be calculated all the way to the final point in the attached worksheet.
  6. Something interesting happens shortly after t=1000 - I'd like to think that this is "real" but I am aware that the calculation is numerical, and the only defined values were at t=0. The further ones gets from t=0, the more likely one is to see some kind of *weird* numerical artefact due a build-up of numerical errors over time. So I hope you were anticipating something interesting happening around t=1000 - if not, then you have to ask, "Is it real?"
  7. You should be aware that this calculation now takes about 60secs on my Win7, 64-bit system: my rig is not "state-of-the-art", but is pretty good: anthing better than 30sec would impress me and anything up to 5mins would not surprise me - just depends what hardwar you are running on, so wait
  8. I have seen your subsequent post about changing the way the fundamental time variable is defined. |I have no intention of even looking at this until you determine whether or not the attached code fulfils is providing you with the answers you expect.

hw2corr3.mw

Well there are a couple of differences in function definitions between your original code above  hw2_final.mw

(on which my original response was based) and the "newer" version behivor_of_I_(pendulum).mw

I have updated my original code to reflect these changes (see attached), but it doesn't seem to have much impact on the results which are generated

hw2corr2.mw

The original code


E[p]:= sqrt
             ( h*add
                 ( (sol[2](j*h)-sol[1](j*h))^2,
                   j=0..1/h
                 )
             ):

forms the relevant value - what you call the "norm" of the quantity sol[2](j*h)-sol[1](j*h) - so if you want the same "norm" of either of these solutions, then just use either sol[1](j*h) or sol[2](j*h), as in

E[p]:= sqrt
             ( h*add
                 ( sol[1](j*h)^2,
                   j=0..1/h
                 )
             ):

or

E[p]:= sqrt
             ( h*add
                 ( sol[2](j*h)^2,
                   j=0..1/h
                 )
             ):

Still not sure what the "meaning" of this quantity is - but then I wasn't sure of the "meaning" of the equivalent quantity in your original problem

This *looks* like a 2-D math input issue - If I convert to 1-D math input, then everything seems to work (at least in Maple 2015.1). See attached below

labelProb.mw

Unfortunately the remainder of my life is too short to debug Maple's 2-D math input problems

You *might* be in a text entry mode - in other words just entering pure, non-executable, text. From the menus, try

Insert->Execution group->After Cursor

This should produce a "command prompt" in the form of a '>' symbol. Type a command at this symbol, finishing with a semi-colon (';'), then C/R and see what happens

You state

"It's the same question"

In which case I suggest that you use the "same" answer. In other words, the solution I posted before, relinked here for your convenience

hw2corr.mw

This works for epsilon =1/20, 1/40, 1/60, 1/80, 1/100, 1/200, 1/400, 1/600, 1/800/ 1/1000, 1/2000 at which point I got bored checking.

If you provide the circumstances in which this code fails, I *might* get interested.

I have no intention of debugging more of your code in order to solve a problem which has already been solved

so if I gave you the equation

a*x2+b*x+c=0

and you "split the equation by its coefficients", what conclusion do you draw??? Perhaps

a*x2=0

b*x=0

c=0

Because if you believe this, then there is reslly nothing I can do to help you.

If you believe that _y1 represents a differential then I suggest you try

_y1(x)=2*x:
dsolve(%);

and contemplate why dsolve() reports that you do not have a differential equation - of course if your "notes" say otherwise, then your "notes" must be correct and Maple must be wrong - or perhaps not

( If I wanted to set up a definition so that _y1(x) was equivalent to diff(y1(x), x) then I could do so - however you do not, so for your worksheet, it simply is not true )

 

People have (in the past had difficulty intalling the help for the DirectSearch package - in case this happens to you, check out

http://www.mapleprimes.com/questions/206551

First 167 168 169 170 171 172 173 Last Page 169 of 207