9 years, 99 days

## Supplied worksheet has two problems...

1. You have used x both as the independent variable in all of the ODEs, but also as an index in a for-loop. The latter use will result on 5.2 being assigned to 'x' after the loop has executed. This will effectively substitute 5.2 everywhere you use 'x' in your ODEs. You really do not want this, so I have changed the offending loop index at the point highlighted in the attached
2. Second problem is more subtle. Location is also highlighted in the attached. The loop which solves the ODEs for coefficients of p0, p1,p2,...etc. For p0 one gets an explicit expression fro c(x), whihc is good. However for p1, the relevant ODE contains c(x) (which will be correctly substituted from the previous loop iteration), c(x) and f(x). It is therefore not possible to get an explicit expression for c(x). Maple will provide a 'formal' solution for c(x) in terms of double integrals of the unknown f(x) - but I'm guessing that this is not what you want! It gets worse! This formal soultion for c(x) is used in the ODE for p2 which also contains the two unknowns c(x) and f(x). Maple will provide another 'formal' solution, but this contains the previous 'formal' solution for c(x) - so ones gets an even more complicated formal solution containing a formal solution. This situation gets worse at each loop iteration, and the expression complexity "explodes". I'm guessing thsat this is not what you want, but it is an inevitable consequence of the way the problem is formulated. Can only suggest that you check the logic of your algorithm very carefully > # # added restart, just because it's always a good idea #   restart;
 > # # Added this just so that worksheet will display # correctly ion the MApleprimes website #   interface(rtablesize=10):
 >  (1)
 >  (2)
 >     (3)
 >  (4)
 >  (5)
 >   (6)
 > > >  (7)
 >   > > > >       (8)
 > > Download hpm2.mw

## Since you have an IVP,...

it is possible to use the 'parameters' option with the dsolve/numeric command.

This will produce a solution module with explicit values of 'p', 'lambda' and 'rho' as yet unspecified

These parameters can be subsequently given values, and relevant outputs obtained, as in the following.

Note that depending on the parameter values which are used, the ode system may (or may not) exhibit a singularity in the solution process

 > restart;
 > with(plots):   interface(rtablesize=10):
 > N:= 10:   beta:= 1:   alpha:= (N + 2)/(N - 2):   ics:= {U(0) = beta, V(0) = 1, X(0) = 0, Y(0) = 1}:   sys_ode:= { diff(U(r), r) = r^(N - 1)*rho^(N - 2)*p*X(r),               diff(V(r), r) = r^(N - 1)*rho^(N - 2)*p*Y(r),               diff(X(r), r) = -r^(N - 1)*rho^N*(U(r)^(alpha - 1) + lambda*U(r)),               diff(Y(r), r) = -r^(N - 1)*rho^N*(lambda + (2^alpha - 1)*U(r)^(alpha - 2))*V(r)             }:   sol:= dsolve( [ sys_ode[], ics[] ],                   numeric,                   parameters = [lambda, p, rho]               ):
 > # # List the unknown parameters (just as a check!) #   sol(parameters); # # Supply one set of parameter values. Note # this is dependent on order. # # Then plot the solution #   sol( parameters = [1, 1, 1] ):   odeplot( sol,            [ [r, U(r)], [r, V(r)], [r,X(r)], [r, Y(r)] ],              r=0..1,              color=[red, blue, green ,black]          ); # # Supply another set of parameter values and # plot the solution #   sol( parameters = [2, 2, 2] ):   odeplot( sol,            [ [r, U(r)], [r, V(r)], [r,X(r)], [r, Y(r)]],              r=0..1,              color=[red, blue, green ,black]          ); # # Supply another set of parameter values and # plot the solution #   sol( parameters = [1, 2, 3] ):   odeplot( sol,            [ [r, U(r)], [r, V(r)], [r,X(r)], [r, Y(r)]],            r=0..1,            color=[red, blue, green ,black]          );    > Download odePars.mw

.

## Not convinced...

that this PDE is separable, in either the original or simplified (with supplied time dependency) versions.

Output of the PDETools:-separability() command suggests otherwise.

This command returns the conditions which have to be satisfied in order for the PDE to be separable. So if it returns exactly 0 conditions, then PDE is separable. For your case the original PDE returns 3 (lengthy) conditions: the "simplified" version returns 1 condition

See the attached

 > > interface(rtablesize=10):
 > > >  (1)
 > >
 > > > > # # Check separability of original PDE. Returns # three conditions which would have to be met # for PDE to be separable #   sep:=[PDETools:-separability(PDE, f(t,r,phi), `*`)]:   numelems(sep); (2)
 > # # Simplify the PDE using an explicit time function #   PDE2:=simplify(expand(subs( f(t,r,phi)=exp(-lambda^2*alpha*t)*g(r,phi), PDE))): # # Check the separability of this simplified PDE # # Still doesn't separate - one condition left #   sep:=[PDETools:-separability(PDE2, g(r,phi), `*`)]:   numelems(sep); (3)
 >

Download PDEIssue.mw

## You cannot use...

the boundary condition

theta(infinity) = 0

when looking for the nummeric solution to an ODE. Just change it ot something like

theta(100) = 0

as in the attached

 > restart;
 > interface(rtablesize=10):
 > BVP1 := { diff( theta(eta), eta, eta) +2 * eta * diff(theta(eta), eta), theta(infinity) = 0, theta(0) = 1 };   #Analytically:     sol1 := dsolve(BVP1);   plot( rhs(sol1), eta=0..10);   > # # Define "infinity" as 100 to make # the boundary condition # # theta(infinity) = 0 # # usable in a "numerical" context #   inf:=100:   BVP2 := { diff( theta(eta), eta, eta) +2 * eta * diff(theta(eta), eta), theta(inf) = 0, theta(0) = 1 }:   sol:=dsolve(BVP2, numeric, maxmesh=1024);   plots:-odeplot(sol, [eta, theta(eta)], eta=0..10);  >

Download odeProb.mw

## Code executes...

with no errors - see the attached

Problem *may* be a Maple version issue. What version are you running?

 > restart:
 > interface(rtablesize=10): kernelopts(version); Physics:-Version(); with(geometry): with(plots): S:=segment: L:=line: Per:=PerpendicularLine: R:=5: xA:=0: yA:=0: point(A,xA,yA): xI:=R/3: yI:=0: point(I1,xI,yI): circle(C,[A,R]): quadri:=proc(t)               local xM,yM,xN,yN,xE,yE,dr1:               xM:=evalf(R*cos(t)):               yM:=evalf(R*sin(t)):               point(M,xM,yM):               L(lMI,[M,I1]):               intersection('h',C,lMI,[M,N]):               L(lAM,[A,M]):               L(lAN,[A,N]):               Per(lME,M,lAM):               Per(lNE,N,lAN):               intersection(E,lME,lNE):               S(sAM,[A,M]):               S(sAN,[A,N]):               S(sME,[M,E]):               S(sNE,[N,E]):               dr1:=draw({lMI(color=blue),sAM(color=black),sAN(color=black),                          sME,sNE}):               display({dr1}):          end:     dr:=draw({C},view=[-6..17,-10..6]):  display([dr,quadri(0.7),quadri(1),quadri(1.2)],view=[-6..17,-10..6]);   >

Download geo.mw

## A Caveat...

The help page at ?fsolve/details states (my emphasis), for the option 'fulldigits'

fulldigits
Prevent fsolve from decreasing Digits for intermediate calculations at high settings of Digits.  With this option, fsolve may escape ill-conditioning problems, but is slower.

without specifying how "high" the setting of "Digits" has to be before reduction of accuracy for intermediate calculations is applied. I doubt that this is likely to affect your calculations, but it is probably safer to use the 'fulldigits' for fsolve() in your calculation.

Based on the desired number of digits (ND) , you can define an acceptable "tolerance" as 1.0*10-ND, but this will no take account of "rounding" in intermediate calculations, so a better definition of the tolerance might be 1.0*10-ND-2.

The just keep incrementing the setting of Digits until the difference between the value obtained with the current setting and that obtained with the previous setting is less than the specified tolerance. Once this criterion is met, round the value obtained to the desired number of digits. See the attached.

 > restart:
 > interface(rtablesize=10): # # Specify required number of digits numDig # so the "best case error" will be # # 1.0*10^(-numDig) # # However the worst case error may be # substantially greater than this, due to # rounding at every stage of the calculation # so define the "desired" error to be a couple # of digits better, as in #   numDig:= 10:   eps:= 1.0*10^(-numDig-2):   f:= 2*sin(x^2/2)-sin(1*x/2)^2:   df:= diff(f,x):   dff:= diff(f,x,x): # # Define a function whih computes the required # quantity.Use the 'fulldigits' option to ensure # that fsolve() is always using the current # setting of Digits internally #   do_dff:= ()-> eval                 ( dff,                   x = fsolve                       ( df,                         x = 1..2,                         fulldigits                       )                 ): # # Compute the desired answer using default # Digits setting. #   old:= do_dff(): # # Increment 'Digits' and repeat the calculation # until the difference between the value on the # current setting of Digits and that on the previous # setting is less than the desired error #   for j from 1 do       Digits:= Digits+1:       new:= do_dff();       if   abs(new-old) < eps       then break       else old:= new:       fi:   od: # # Out of idle curiosity, check the value of Digits # where the error criterion was met, and the value # obtained. #   Digits;   new; # # Round the obtained value to the required number # of Digits #   evalf[numDig](new);   (1)
 >

Download dff10.mw

## Version issue...

the code

assume(x, posint);
binomial(x, -1);
assume(x, Non(integer));
binomial(x, -1);

returns

0
0

in Maple 2017 and later versions: but returns

0
binomial(x~, -1)

in earlier versions. Before investigating furthr for a workaround, it would be useful to know precisely which Maple version your are running

## As well as checking the addresses...

you can use the about() function which will show (for your example) that the two instances of the variable 'b' have different associated assumptions.

Should additionally() work this way - well I can think of arguments for and against. But it is certainly a trap for the unwary!

A couple of useful general rules

1. If you are going to make assumptions on variables in a worksheet, do all of them right at the start
2. Certain commands (such as solve() ) will ignore pre-existing assumptions unless you invoke the 'useassumptions' option

## Like this ?...

 > restart;
 > # # Ignore this execution group. It is only necessary # to ensure that the worksheet displays "inline" on # the Mapleprimes website #   interface(rtablesize=10):
 > # # Define a vector of data. Becuase it is being # generated one has to specify a "length" for the # vector. However all subsequent commands in # this worksheet work for any "length" of vector # # Note that rather than generating data, one could # import a vector from Excel using ExcelTools:-Import(), # or from a lot of other formats using ImportVector() #   V:=LinearAlgebra:-RandomVector(2880): # # Do it the easy way with commands from the # statistics package #   Statistics:-Mean(V);   Statistics:-Variance(V);   Statistics:-StandardDeviation(V); # # Do it the hard way - from "first principles" #   evalf(add(V)/numelems(V));   evalf(add((i-%)^2, i in V)/(numelems(V)-1));   sqrt(%);      (1) Download stats1.mw

In order to have a matrix with "blank" entries, you need to use 'NULL', as in th following

 > restart:
 > # # Ignore this execution group. It is only necessary # to ensure that the worksheet displays "inline" on # the Mapleprimes website #   interface(rtablesize=10):
 > # # Define a matrix with 'NULL' entries #   M:=Matrix(3,3,(i,j)->`if`(i=j,i+j, NULL)); (1)
 > # # Access entries of M #   M[1,1]; # returns 2   M[1,2]; # returns nothing at all (2)
 >

Download nullMat.mw

## Errrrrr...

and you can display results within a Maple worksheet, or (assuming that they are 2-dimensional) can be easily exported to Excel, so I don't really understand the issue. A few possible pointers

1. If commands are terminated with a colon character (ie ':') then output will be generated but not displayed, which controls the problem of "screenfuls" of output. If you definitely want out to be displayed then terminate commands with a semicolon character (ie ';')
2. By default Maple will not display 2D structures grater than 10 X 10. Again this is to avoid the issue of "screenfuls" of output. This limit can be overridden by using the interface(rtablesize) command
3. The ExcelTools:-Export() command will export most 2-D data structures in Excel format

What else do you need??

## An explanation...

of the thought process behind my original post.

I assumed that you wanted (somehow) to come with a "simple polynomial" which would closely approximate the graph dislayed in your original post

A more detailed explanation of the steps I originally performed is given in the attached worksheet. The following just summarises the basic issues

1. You choose 11 conditions, which requires a 10-th order polynomial in order to produce enough variables for fitting purposes.
2. One of these t(2.8)=0 is incorrect: I assume you meant t(-2.8)=0 which corresponds to the graph
3. Your (corrected) fitting points extend over the range [-2.8..3.5], but the graph extends over the range [-3.5, 4.0]. You should never EVER expect to fir something on one range and then be able to extrapolate. This is pretty much guaranteed to fail in the extrapolated ranges, ie [-3.5..-2.8] and [3.5..4.0] - and the higher the order of the polynomial, the bigger the FAIL
4. It is further true that whiilst you will fit accurately all the desired roots and extremae, a 10-th order polynomial will have further roots and extremae, and you cannot guarantee that these will not occur in the range where you interested
5. Adding in the two end-points (ie x=-3.5 and x=4) and increasing the polynomila order to 12 addresses (3) above, but not (4), so doing this still result in an very bad fit
6. Since the basic problem at this stage was the "spurious" extrema and structure arising  becuase of (4) above, I decided to get "innovative", with no purely mathematical justification. This involved using the lowest-order polynomials possible which were likely to fit the available data - more or less why cubic splines exist. I don't propose to get into a discussion here on using higher (or even lower!) order splines
7. Using splines in this way has two obvious problems
1. The resultant function is piecewise, although it will be zeroth- and first-order continuous
2. There is no way to specify that a particular point is a maximum/minimum: for example it is trivial to ensure that the resulting curve will pass through the point [-2,2], there is no criterion (I know of) whihc will force this point to be a maximum
8. As things turned out, some of the extremae are slightly "off" - but really not by much, and the overall fit of the spline function is much better than any other fit I have seen posted here
9. Having obtained a spline fit, it is then realtively trivial to generate an arbitrary number of points in the desired range, and then experiment with polynomial orders whihc produce a (more-or-less) equally good fit. In these experiments, a 10-th order polynomial isn't too bad, but a 12-th order polynomial is better. Increasing the polynomila order beyond this doesn't (subjectively) improve things very much

Check the attached - a somewhat extended version of my original response - although the answers are more or less the same.

 > restart;   interface(rtablesize=10); # # Op's original code, rewritten slightly for convenience #   f:=z-> local j; add( a[j]*z^j, j=0..10);   sol:=solve( [ f(-2)=2,   f(-1)=0,  f(0)=-2, f(1)=0,                 f(2.5)=-3, f(-2.8)=0, f(3.6)=0,                 D(f)(-2)=0, D(f)(0)=0, D(f)(1)=0, D(f)(2.5)=0             ]           );   assign(sol);   plot(f, -2..3.6);   plot(f, -3.5..4.0)     > # # The above is obviously crap outside the range # [-2..3,6] and isn't very good within this range! # # There is a "spurious' maximum arounf x=3.2 and # some "doubtful" extra "structure" between x=1 # and x=2. # # My first thought for *improving* this fit was to # incorporate the "endpoints" at x=-3.5 and x=4 # # If anything this probably made things worse! There # is a spurious minimum at -3.1, a spurious maximum # at 3.3 and some "doubtful* extra structure between # x=1 and x=2 #   restart;   f:=z-> local j; add( a[j]*z^j, j=0..12);   sol:=solve( [ f(-3.5)=-5, f(4)=3,    f(-2)=2,   f(-1)=0, f(0)=-2,                 f(1)=0,     f(2.5)=-3, f(-2.8)=0,  f(3.6)=0,                 D(f)(-2)=0, D(f)(0)=0, D(f)(1)=0, D(f)(2.5)=0             ]           );   assign(sol);   plot(f, -3.5..4);   plot(f,  -2..3.6);    > # # Having failed completely with "simple" polynomial fits # I decided to be a little "creative". Since the curve # presented in the OP's original post is "simple", ie nice # and smooth between various critical points (ie roots and # extrema), I decide to investigate what would would # happen if I used the lowest order polynomials possible # which would fit all the relevant cirical values, and # including the end points. # # The obvious choice was to use cubic spline - originally # I just wanted to see how close I could get to original # posted curve using splines, so I came up with # # One problem with using splines is that there is no way # (I know of!) to specify extrema - eg one can require # that the curve goes through the point [-2, 2], but # cannot ensure that this point will actually be a maximum. # Therefore one can anticipate that exact positioning of # maximae and minimae will be a little "off", but (I # thought) what the hell - let's see how close one can get # # And the answer was pretty damn close over the whole # fitting range! #   restart;   sol:= unapply         ( CurveFitting:-Spline           ( [ [-3.5,-4], [-2.8,0], [-2,2], [-1,0], [0,-2], [1,0], [2.5,-3], [4,3] ],             x           ),           x         );   p1:= plot( sol,              -3.5..4,              color=red,              title = typeset("Spline Fit"),              titlefont = [times, bold, 20],              size = [500,500]            );  > # # The above spline fit relies on a piecewise definition, # which will be zeroth and first-order continuous at the # piecewise boundaries. But the OP "appeared" to want a # polynomial fit. Therefore I decide to use the above # spline fit to generate a "gazillion" points acros the # desired range [-3.5, 4], and the use these points to fit # a polynomial - rather than relying on eight or so # obvious critical points # # There is little "mathematical" justification for this. I # just wanted to see how close I could get. Hence the # following. Note that since maximae and minimae position # will be a little "off" in the spline fit, they will # almost certainly be slightly "off" in the polynomial fit # # Generate points to fit from the spline equation above #   XYData:= Matrix            ( [ seq                ( [i, sol(i)],                  i=-3.5..4, 0.1                )              ]            ): # # Specify the order of the polynomial to fit #   deg:=12:   poly:= unapply          ( Statistics:-Fit            ( add              ( a[j]*x^j,                j = 0..deg              ),              XYData,              x            ),            x          ); # # And plot the resulting polynomial #   p2:= plot( poly,              -3.5..4,              color = red,              title = typeset("Polynomial Fit"),              titlefont = [times, bold, 20],              size = [500, 500]            );  > # # Just for comparison purposes, plot the spline # and polynomial fits on the same graph #   plots:-display          ( [ p1, p2],            color = [red, blue],            legend = [ typeset("Spline"),                       typeset(deg, "-th order Polynomial")                     ],            legendstyle = [ font=[times, bold, 16] ],            title = typeset("Spline and Polynomial Fit"),            titlefont = [times, bold, 20],            size = [500, 500]       ); >

Download doFit2.mw

## Maybe somethinng in the attached...

covers your requirement.

The overall strategy is to identify about six "important" points and fit to these using (cubic) splines. This returns a pretty good fit, but by its nature, teh fitted function will be piecewise. However this piecewise function can be used to generate "lots" of data points. The latter can then be fitted to a (continuous) polynomial function

 > restart;   interface(rtablesize=10):   with(CurveFitting):   with(Statistics):
 > # # Fit the curve using splines, just by # "eyeballing" specific values. This will # lead to a piecewise function # # Accuracy of fit is very dependent on how # well I estimated the following points! #   sol:= unapply         ( Spline           ( [ [-3.5,-4], [-2,2], [-1,0], [0,-2], [1,0], [2.5,-3], [4,3] ],             x           ),           x         ); # # Plot the fit to see how "good" it is #   plot( sol,         -3.5..4,         title = typeset("Spline Fit"),         titlefont = [times, bold, 20],         size = [500,500]       );  > # # Use the spline (piecewise) fit to generate a "large" # set of data points, and then fit these to a polynomial # of degree given by 'deg'. Although the fit with deg:=10 # isn't *too* bad, really eed deg:=12 to get a *good* fit, #   deg:=12: # # Generate points to fit from the spline equation above #   XYData:= Matrix            ( [ seq                ( [i, sol(i)],                  i=-3.5..4, 0.1                )              ]            ):   poly:= unapply          ( Fit            ( add              ( a[j]*x^j,                j = 0..deg              ),              XYData,              x            ),            x          ); # # Plot the spline fit and the polynomial fit on the # same graph #   plot( [ sol, poly ],          -3.5..4,           color = [red, blue],           legend = [ typeset("Spline"),                      typeset(deg, "-th order Polynomial")                    ],           legendstyle = [ font=[times, bold, 16] ],           title = typeset("Spline and Polynomial Fit"),           titlefont = [times, bold, 20],           size = [500, 500]       );  >
 >

Download fitStuff.mw

## Use the appropriate...

plot options, as shown in the attached

 > restart;   interface(rtablesize=10):
 > sys:= tau-> [ diff(x(t), t) = x(t)*(2-(1/20)*x(t))-x(t-tau)*y(t)/(x(t)+10)-10,                 diff(y(t), t) = y(t)*(x(t)/(x(t)+10)-2/3),                 x(0) = 40,                 y(0) = 16               ]:   doplot1:= sol->  plots[odeplot]                    ( sol,                      [[t,x(t)], [t,y(t)]],                      0 .. 10,                      labels=[typeset("t"), typeset("x(t), y(t)")],                      labelfont=[times, bold, 20]                    ):   doplot2:= sol-> plots[odeplot]                   ( sol,                     [x(t), y(t)],                     0 .. 10,                     labels=[typeset("x(t)"), typeset("y(t)")],                     labelfont=[times, bold, 20]                   ):
 > # # Solution for tau=0 #   sol := dsolve( sys(0), numeric): # # Plot of x(t) and y(t) versus t #   doplot1(sol); # # Plot of y(t) versus x(t) #   doplot2(sol);  > # # Solution for tau=0.4 #   sol := dsolve( sys(0.4), numeric): # # Plot of x(t) and y(t) versus t #   doplot1(sol); # # Plot of y(t) versus x(t) #   doplot2(sol);  >
 >

Download delODE.mw

## Maybe...

one of the options in the attached meets your requirements

Once again this website is refusing to display a worksheet - you will have to download to view - sorry

## Something like...

that shown in the attached, which (again!) this website is failing to display.

Download aPlot.mw

﻿