Carl Love

Carl Love

28085 Reputation

25 Badges

13 years, 97 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

It'll be much easier to use calculus than to use Maple's CDF command for this. Indeed, I can do the calculus in my head.

Hint: The integral of any PDF over its support equals 1. (The support of a PDF is the values of its independent variable (in this case) where the PDF is not 0.) So, it can be done with two lines of Maple code: the first uses solve and int to get the value of c. The second uses int and piecewise to get the final answer.

This is usually a sign of a licensing problem. You'll probably need to directly contact Maplesoft's customer support to get it fixed.

One of the basic techniques of computer algebra is to solve the same problem several times over different prime fields and then "lift" those solutions to the integers. There are several algorithms for the lifting. Here, I use the oldest and most well known: Chinese remaindering (see help page ?chrem).

The code below finds an 20-digit (67-bit) value of d in 2.3 seconds. That's larger than wordsize. This only requires 13 primes, all less than 80. It attempts only 619 factorizations over finite fields and 13 irreducs over the integers.
 

restart
:

First, generate a suitably complicated example:

p1,p2:= 'randpoly([x,y], coeffs= rand(-2^33..2^33))'$2;

7292118554*x^2*y^3-6262890291*x^4+551415575*x*y^3-5966975316*x^2*y-6044637241*x*y^2+4637448791, 6492083754*x^3*y^2+1792855492*x^2*y^3+5169280949*x^4-4787552440*x*y^3-3142438976*y^3+7094595294*x*y

Check that both factors are irreducible (or else the problem is trivial):

irreduc(p1), irreduc(p2);

true, true

pp:= expand(p1*p2);

47341044396665371716*x^5*y^5+13073714797853998568*x^4*y^6-40659208311285432414*x^7*y^2+26466552265028799574*x^6*y^3+3579836096160068550*x^4*y^5-33922791533958883860*x^3*y^6-32374639466943366159*x^8-5903765788563875549*x^5*y^3-49940215697038518186*x^4*y^4-10837161074674577572*x^3*y^5-25554968543247613704*x^2*y^6-30844971824152054884*x^6*y-31246428133517221709*x^5*y^2+19680750552850382016*x^4*y^3+80301837210034055916*x^3*y^4+28939017772064418040*x^2*y^5-1732789794853451200*x*y^6-44432671985366890554*x^5*y+22662926145261620466*x^2*y^4+18994903661899505216*x*y^5-12226569040249721490*x^3*y^2-34569979390122633682*x^2*y^3+23972275687279382659*x^4-22202029274727100040*x*y^3-14572899830042478016*y^3+32900822368794589554*x*y

Check for integer factors:

content(pp);

1

Select a coefficient to become d:

m:= max(coeffs(pp));

80301837210034055916

p:= subs(m=d, pp);

47341044396665371716*x^5*y^5+13073714797853998568*x^4*y^6-40659208311285432414*x^7*y^2+26466552265028799574*x^6*y^3+3579836096160068550*x^4*y^5-33922791533958883860*x^3*y^6+d*x^3*y^4-32374639466943366159*x^8-5903765788563875549*x^5*y^3-49940215697038518186*x^4*y^4-10837161074674577572*x^3*y^5-25554968543247613704*x^2*y^6-30844971824152054884*x^6*y-31246428133517221709*x^5*y^2+19680750552850382016*x^4*y^3+28939017772064418040*x^2*y^5-1732789794853451200*x*y^6-44432671985366890554*x^5*y+22662926145261620466*x^2*y^4+18994903661899505216*x*y^5-12226569040249721490*x^3*y^2-34569979390122633682*x^2*y^3+23972275687279382659*x^4-22202029274727100040*x*y^3-14572899830042478016*y^3+32900822368794589554*x*y

Need this utility procedure for Cartesian products:

CartProdSeq := proc(L::seq(list))
local Seq,i,j;
option `Copyright (C) 2007, Joseph Riel. All rights reserved.`;
    eval([subs(Seq=seq, foldl(Seq
                              , [cat(i,1..nargs)]
                              , seq(cat(i,j)=L[j],j=nargs..1,-1)
                             ))]);
end proc
:

This procedure finds suitable values of d over small prime fields by brute force and uses Chinese remaindering to "lift" the ds to the integers.

Find_d:= proc(p::polynom(integer), d::name)
local
    pf:= 2, Es:= Array(1..0),Em:= Array(1..0), Ms:= Array(1..0),
    Mm:= Array(1..0), C:= {coeffs(p)}, Ep, i, ML, e, dt, EL
;
    `mod`:= mods;
    do
        pf:= nextprime(pf);
        if 0 in C mod pf then next fi;
        Ep:= Array(1..0);
        for i from (1-pf)/2 to (pf-1)/2 do
            if (Factor(eval(p, d= i)) mod pf)::{`*`, `^`} then
                Ep,= i
            fi
        od;
        if numelems(Ep)=0 then
            return "Impossible to factor"
        elif numelems(Ep)=1 then
            Es,= seq(Ep); Ms,= pf
        else
            Em,= [seq(Ep)]; Mm,= pf
        fi;
        ML:= [seq(Ms), seq(Mm)]; EL:= seq(Es);
        for e in CartProdSeq(seq(Em)) do
            dt:= chrem([EL,e[]], ML);
            if not irreduc(eval(p, d= dt)) then
                userinfo(1, Find_d, ML, [EL,seq(Em)]);
                return dt
            fi
        od
    od
end proc
:   

infolevel[Find_d]:= 1:

CodeTools:-Usage(Find_d(p, d));

Find_d: [19 23 29 31 37 41 43 53 59 61 71 73 79] [-2 10 -14 -10 -9 -15 21 21 -19 -9 -10 0 19]

memory used=210.06MiB, alloc change=0 bytes, cpu time=2.30s, real time=2.11s, gc time=421.88ms

80301837210034055916

m;

80301837210034055916

ilog10(m);

19

ilog2(m);

66

 


 

Download ChineseLifting.mw

If the result containing xs were named res, you could do

eval(res, x= x(t))

There are some shorter and cruder possibilities for your particular example, but they won't work for more complicated examples.

When you get a non-smooth implicitplot, the first thing to try is the gridrefine option. It's a positive integer defaulting to 1. Here, I've set gridrefine= 3 for your plot:

The minimum is 0 at a=0; the value of is irrelevant. It's not a local minimum because a=0 is one of the boundaries. This can be confirmed by simple plotting:

plot3d(f, 2..5, -1..0)

An iron-clad proof that this is the minimum will be harder to obtain.

So, the Optimization package was more accurate in this case.

Here is VV's algorithm as a predicate (i.e., returns true or false) procedure. This doesn't rely on the vertices being numeric and tries to do a minimal amount of copying:

IsOuterplanar:= proc(G::Graph)
uses GT= GraphTheory;
    GT:-IsPlanar(GT:-GraphJoin(G, GT:-PathGraph(1)))
end proc
:

 

The cases that don't work via a direct built-in distributive property can be done with map:

e:= a <= b;
map(exp, e);
map(`/`, e, -2);

Caution: This works regardless of whether the results are mathematically valid.

Many of these distributive properties are documented at ?evalapply.
 

Using solve(..., parametric), a definitive programmatic solution can be obtained in Maple 2019:

f:= 2*x^5-x^3*y+2*x^2*y^2-x*y^3+2*y^5;
solve({f, x < 0, y < 0}, parametric);

      [ ]

Note that those brackets won't be empty if either or both inequalities are reversed.

The parametric option only works for polynomial systems and systems that can be easily reduced to polynomials.

I find ColumnGraph difficult to use when you have non-categorical horizontal-axis data (such as real-number measurements), and it becomes clear that this wasn't its intended usage. Instead, it can be done with a regular plot command like this:

plot(
   ((x,y)-> [[x,0],[x,y]])~(yr, sle), 
   thickness= 70, xtickmarks= yr, view= [0..0.35, DEFAULT]
);
 

The colors can be easily changed if you want.

Just include fclose("input.txt") in the procedure.

I think that in most languages, most opened files stay open until they are explicitly closed.

 

Here is a much-faster and slightly more accurate way to solve the BVP by finite-difference methods. This is a method that I just came up with: Start by approximating the function value at the midpoint of the interval. So now we have two subintervals of length (c-a)/2. On each iteration:

  1. Divide all subintervals in half.
  2. Approximate the function value at each odd-numbered node (the s[k]) using the even-numbered nodes (s[k-1] and s[k+1]) adjacent to it, which have already been approximated. The odd-numbered nodes are the midpoints of the new subintervals.
  3. Use fsolve on the entire set of finite-difference equations to refine the approximations (both odd and even nodes), using the existing approximations as initial guesses to fsolve.

This method is limited by the fickleness of fsolve. In actual practice, I'd code a Newton's method solver to take the place of fsolve for this. Nonethless, the worksheet below gets a 127-point solution in a few seconds with a maximum absolute error of 6e-7, (or 9 significant digits). That's pretty good!
 

restart
:

Digits:= 15:

sy := 2400.:
Hp := 0.1e9:
(P,Q):= (-900., -200): #boundary values
(a,c):= (20,21): #solution interval
N:= 7: #Meaning 2^7-1 interior points
Et := (r*(diff(sr(r), r))-sy)/Hp:
Er := (sy-r*(diff(sr(r), r)))/Hp:
ode := simplify(Er-Et-r*(diff(Et, r))*(1+(diff(Et, r))*r/(2+4*Et))):

#Solve for highest-order derivative:
ODEs:= solve(ode, diff(sr(r),r$2));

(-3.*r*(diff(sr(r), r))-99995200.+2.*(-1.*r^2*(diff(sr(r), r))^2+4800.*r*(diff(sr(r), r))+0.249999999424000e16)^(1/2))/r^2, (-3.*r*(diff(sr(r), r))-99995200.-2.*(-1.*r^2*(diff(sr(r), r))^2+4800.*r*(diff(sr(r), r))+0.249999999424000e16)^(1/2))/r^2

ODE:= ODEs[1]; #Select one solution of interest.

(-3.*r*(diff(sr(r), r))-99995200.+2.*(-1.*r^2*(diff(sr(r), r))^2+4800.*r*(diff(sr(r), r))+0.249999999424000e16)^(1/2))/r^2

FD_ODE:= subs( #Substitute finite-difference approximations.
    diff(sr(r),r$2)= (s[k+1]-2*s[k]+s[k-1])/h^2,
    diff(sr(r),r)= (s[k+1]-s[k-1])/2/h,
    sr(r)= s[k],
    diff(sr(r),r$2) = ODE
);

(s[k+1]-2*s[k]+s[k-1])/h^2 = (-1.50000000000000*r*(s[k+1]-s[k-1])/h-99995200.+2.*(-.250000000000000*r^2*(s[k+1]-s[k-1])^2/h^2+2400.00000000000*r*(s[k+1]-s[k-1])/h+0.249999999424000e16)^(1/2))/r^2

sk:= solve(FD_ODE, s[k]); #Solve for s[k] ~ sr(r).

-.250000000000000*(4.*(-.250000000000000*r^2*(s[k+1]-1.*s[k-1])^2/h^2+2400.*r*(s[k+1]-1.*s[k-1])/h+0.249999999424000e16)^(1/2)*h^2+3.*h*r*s[k-1]-3.*h*r*s[k+1]-2.*r^2*s[k-1]-2.*r^2*s[k+1]-199990400.*h^2)/r^2

#Use s[k] and b.v.s to approximate midpoint of interval:
F(a+(c-a)/2):= eval(sk, [h= (c-a)/2, r= a+(c-a)/2, s[k-1]=P, s[k+1]=Q]);

-538.621994036408

F(a):= P:  F(c):= Q:

#Iteratively divide subintervals in half, approximating midpoints, and refining the
#approximations with fsolve:
for j from 2 to N do
    h:= (c-a)/2^j;
    for i to 2^(j-1) do
        r:= a+(2*i-1)*h;
        F(r):= eval(sk, [s[k-1]= F(r-h), s[k+1]= F(r+h)])
    od;
    r:= 'r';
    fsol:= fsolve(
        eval(
            {seq(eval(s[k]=sk, [k= _k, r= a+h*_k]), _k= 1..2^j-1)},
            [s[0]= P, s[2^j]= -200]
        ),
        {seq(s[_k]= F(a+_k*h), _k= 1..2^j-1)}
    );
    assign(seq('F'(a+_k*h)= eval(s[_k], fsol), _k= 1..2^j-1))
od;

1/4

r

{s[1] = -716.331363276185, s[2] = -538.623740093766, s[3] = -366.600201340521}

1/8

r

{s[1] = -807.402783510995, s[2] = -716.331702654157, s[3] = -626.750542820446, s[4] = -538.624176550255, s[5] = -451.918524364314, s[6] = -366.600517235586, s[7] = -282.638060624027}

1/16

r

{s[1] = -853.508345065787, s[2] = -807.402833913394, s[3] = -761.678835056789, s[4] = -716.331787495739, s[5] = -671.357199424526, s[6] = -626.750646968260, s[7] = -582.507772946101, s[8] = -538.624285660747, s[9] = -495.095957713546, s[10] = -451.918624844611, s[11] = -409.088184797350, s[12] = -366.600596206804, s[13] = -324.451877511242, s[14] = -282.638105886437, s[15] = -241.155416202110}

1/32

r

{s[1] = -876.705615839058, s[2] = -853.508351878262, s[3] = -830.407622484913, s[4] = -807.402846513878, s[5] = -784.493447266133, s[6] = -761.678852447755, s[7] = -738.958494129353, s[8] = -716.331808705933, s[9] = -693.798236857205, s[10] = -671.357223508305, s[11] = -649.008217790943, s[12] = -626.750673004960, s[13] = -604.584046580309, s[14] = -582.507800039417, s[15] = -560.521398959977, s[16] = -538.624312938115, s[17] = -516.816015551953, s[18] = -495.095984325560, s[19] = -473.463700693282, s[20] = -451.918649964453, s[21] = -430.460321288467, s[22] = -409.088207620233, s[23] = -387.801805685979, s[24] = -366.600615949426, s[25] = -345.484142578309, s[26] = -324.451893411252, s[27] = -303.503379924991, s[28] = -282.638117201931, s[29] = -261.855623898047, s[30] = -241.155422211120, s[31] = -220.537037849301}

1/64

r

{s[1] = -888.340631566191, s[2] = -876.705616723052, s[3] = -865.094881842957, s[4] = -853.508353581365, s[5] = -841.945958875511, s[6] = -830.407624943098, s[7] = -818.893279281000, s[8] = -807.402849663968, s[9] = -795.936264143343, s[10] = -784.493451045786, s[11] = -773.074338971996, s[12] = -761.678856795455, s[13] = -750.306933661165, s[14] = -738.958498984399, s[15] = -727.633482449457, s[16] = -716.331814008430, s[17] = -705.053423879971, s[18] = -693.798242548061, s[19] = -682.566200760802, s[20] = -671.357229529205, s[21] = -660.171260125981, s[22] = -649.008224084346, s[23] = -637.868053196831, s[24] = -626.750679514094, s[25] = -615.656035343739, s[26] = -604.584053249152, s[27] = -593.534666048323, s[28] = -582.507806812698, s[29] = -571.503408866015, s[30] = -560.521405783163, s[31] = -549.561731389036, s[32] = -538.624319757401, s[33] = -527.709105209765, s[34] = -516.816022314253, s[35] = -505.945005884493, s[36] = -495.095990978500, s[37] = -484.268912897571, s[38] = -473.463707185189, s[39] = -462.680309625923, s[40] = -451.918656244346, s[41] = -441.178683303944, s[42] = -430.460327306049, s[43] = -419.763524988761, s[44] = -409.088213325886, s[45] = -398.434329525870, s[46] = -387.801811030752, s[47] = -377.190595515111, s[48] = -366.600620885024, s[49] = -356.031825277028, s[50] = -345.484147057086, s[51] = -334.957524819564, s[52] = -324.451897386208, s[53] = -313.967203805121, s[54] = -303.503383349760, s[55] = -293.060375517928, s[56] = -282.638120030772, s[57] = -272.236556831790, s[58] = -261.855626085840, s[59] = -251.495268178153, s[60] = -241.155423713356, s[61] = -230.836033514497, s[62] = -220.537038622073, s[63] = -210.258380293071}

1/128

r

{s[1] = -894.167267028199, s[2] = -888.340631678682, s[3] = -882.520084721372, s[4] = -876.705616943965, s[5] = -870.897219151880, s[6] = -865.094882168244, s[7] = -859.298596833823, s[8] = -853.508354007005, s[9] = -847.724144563739, s[10] = -841.945959397506, s[11] = -836.173789419279, s[12] = -830.407625557475, s[13] = -824.647458757920, s[14] = -818.893279983809, s[15] = -813.145080215662, s[16] = -807.402850451286, s[17] = -801.666581705737, s[18] = -795.936265011275, s[19] = -790.211891417331, s[20] = -784.493451990461, s[21] = -778.780937814309, s[22] = -773.074339989569, s[23] = -767.373649633940, s[24] = -761.678857882098, s[25] = -755.989955885647, s[26] = -750.306934813087, s[27] = -744.629785849762, s[28] = -738.958500197837, s[29] = -733.293069076249, s[30] = -727.633483720675, s[31] = -721.979735383484, s[32] = -716.331815333711, s[33] = -710.689714857010, s[34] = -705.053425255618, s[35] = -699.422937848320, s[36] = -693.798243970405, s[37] = -688.179334973634, s[38] = -682.566202226199, s[39] = -676.958837112685, s[40] = -671.357231034033, s[41] = -665.761375407506, s[42] = -660.171261666648, s[43] = -654.586881261242, s[44] = -649.008225657282, s[45] = -643.435286336932, s[46] = -637.868054798488, s[47] = -632.306522556342, s[48] = -626.750681140946, s[49] = -621.200522098774, s[50] = -615.656036992290, s[51] = -610.117217399900, s[52] = -604.584054915927, s[53] = -599.056541150570, s[54] = -593.534667729869, s[55] = -588.018426295669, s[56] = -582.507808505583, s[57] = -577.002806032956, s[58] = -571.503410566831, s[59] = -566.009613811911, s[60] = -560.521407488525, s[61] = -555.038783332590, s[62] = -549.561733095584, s[63] = -544.090248544494, s[64] = -538.624321461797, s[65] = -533.163943645418, s[66] = -527.709106908693, s[67] = -522.259803080339, s[68] = -516.816024004417, s[69] = -511.377761540296, s[70] = -505.945007562620, s[71] = -500.517753961271, s[72] = -495.095992641337, s[73] = -489.679715523079, s[74] = -484.268914541892, s[75] = -478.863581648273, s[76] = -473.463708807788, s[77] = -468.069288001034, s[78] = -462.680311223613, s[79] = -457.296770486090, s[80] = -451.918657813961, s[81] = -446.545965247625, s[82] = -441.178684842342, s[83] = -435.816808668206, s[84] = -430.460328810109, s[85] = -425.109237367706, s[86] = -419.763526455385, s[87] = -414.423188202231, s[88] = -409.088214751994, s[89] = -403.758598263055, s[90] = -398.434330908399, s[91] = -393.115404875574, s[92] = -387.801812366662, s[93] = -382.493545598249, s[94] = -377.190596801385, s[95] = -371.892958221557, s[96] = -366.600622118662, s[97] = -361.313580766961, s[98] = -356.031826455056, s[99] = -350.755351485856, s[100] = -345.484148176547, s[101] = -340.218208858556, s[102] = -334.957525877521, s[103] = -329.702091593261, s[104] = -324.451898379740, s[105] = -319.206938625042, s[106] = -313.967204731331, s[107] = -308.732689114827, s[108] = -303.503384205772, s[109] = -298.279282448394, s[110] = -293.060376300885, s[111] = -287.846658235362, s[112] = -282.638120737838, s[113] = -277.434756308193, s[114] = -272.236557460144, s[115] = -267.043516721208, s[116] = -261.855626632680, s[117] = -256.672879749594, s[118] = -251.495268640697, s[119] = -246.322785888420, s[120] = -241.155424088842, s[121] = -235.993175851664, s[122] = -230.836033800178, s[123] = -225.683990571238, s[124] = -220.537038815227, s[125] = -215.395171196027, s[126] = -210.258380390992, s[127] = -205.126659090915}

interface(rtablesize= 2^N+1):

#Put final approximation in two-column form:
FDM_Sol:= Matrix([lhs,rhs]~(sort(evalf([indices(op(4, eval(F)), pairs)]), key= lhs)));

Matrix(129, 2, {(1, 1) = 20., (1, 2) = -900., (2, 1) = 20.0078125000000, (2, 2) = -894.167267028199, (3, 1) = 20.0156250000000, (3, 2) = -888.340631678682, (4, 1) = 20.0234375000000, (4, 2) = -882.520084721372, (5, 1) = 20.0312500000000, (5, 2) = -876.705616943965, (6, 1) = 20.0390625000000, (6, 2) = -870.897219151880, (7, 1) = 20.0468750000000, (7, 2) = -865.094882168244, (8, 1) = 20.0546875000000, (8, 2) = -859.298596833823, (9, 1) = 20.0625000000000, (9, 2) = -853.508354007005, (10, 1) = 20.0703125000000, (10, 2) = -847.724144563739, (11, 1) = 20.0781250000000, (11, 2) = -841.945959397506, (12, 1) = 20.0859375000000, (12, 2) = -836.173789419279, (13, 1) = 20.0937500000000, (13, 2) = -830.407625557475, (14, 1) = 20.1015625000000, (14, 2) = -824.647458757920, (15, 1) = 20.1093750000000, (15, 2) = -818.893279983809, (16, 1) = 20.1171875000000, (16, 2) = -813.145080215662, (17, 1) = 20.1250000000000, (17, 2) = -807.402850451286, (18, 1) = 20.1328125000000, (18, 2) = -801.666581705737, (19, 1) = 20.1406250000000, (19, 2) = -795.936265011275, (20, 1) = 20.1484375000000, (20, 2) = -790.211891417331, (21, 1) = 20.1562500000000, (21, 2) = -784.493451990461, (22, 1) = 20.1640625000000, (22, 2) = -778.780937814309, (23, 1) = 20.1718750000000, (23, 2) = -773.074339989569, (24, 1) = 20.1796875000000, (24, 2) = -767.373649633940, (25, 1) = 20.1875000000000, (25, 2) = -761.678857882098, (26, 1) = 20.1953125000000, (26, 2) = -755.989955885647, (27, 1) = 20.2031250000000, (27, 2) = -750.306934813087, (28, 1) = 20.2109375000000, (28, 2) = -744.629785849762, (29, 1) = 20.2187500000000, (29, 2) = -738.958500197837, (30, 1) = 20.2265625000000, (30, 2) = -733.293069076249, (31, 1) = 20.2343750000000, (31, 2) = -727.633483720675, (32, 1) = 20.2421875000000, (32, 2) = -721.979735383484, (33, 1) = 20.2500000000000, (33, 2) = -716.331815333711, (34, 1) = 20.2578125000000, (34, 2) = -710.689714857010, (35, 1) = 20.2656250000000, (35, 2) = -705.053425255618, (36, 1) = 20.2734375000000, (36, 2) = -699.422937848320, (37, 1) = 20.2812500000000, (37, 2) = -693.798243970405, (38, 1) = 20.2890625000000, (38, 2) = -688.179334973634, (39, 1) = 20.2968750000000, (39, 2) = -682.566202226199, (40, 1) = 20.3046875000000, (40, 2) = -676.958837112685, (41, 1) = 20.3125000000000, (41, 2) = -671.357231034033, (42, 1) = 20.3203125000000, (42, 2) = -665.761375407506, (43, 1) = 20.3281250000000, (43, 2) = -660.171261666648, (44, 1) = 20.3359375000000, (44, 2) = -654.586881261242, (45, 1) = 20.3437500000000, (45, 2) = -649.008225657282, (46, 1) = 20.3515625000000, (46, 2) = -643.435286336932, (47, 1) = 20.3593750000000, (47, 2) = -637.868054798488, (48, 1) = 20.3671875000000, (48, 2) = -632.306522556342, (49, 1) = 20.3750000000000, (49, 2) = -626.750681140946, (50, 1) = 20.3828125000000, (50, 2) = -621.200522098774, (51, 1) = 20.3906250000000, (51, 2) = -615.656036992290, (52, 1) = 20.3984375000000, (52, 2) = -610.117217399900, (53, 1) = 20.4062500000000, (53, 2) = -604.584054915927, (54, 1) = 20.4140625000000, (54, 2) = -599.056541150570, (55, 1) = 20.4218750000000, (55, 2) = -593.534667729869, (56, 1) = 20.4296875000000, (56, 2) = -588.018426295669, (57, 1) = 20.4375000000000, (57, 2) = -582.507808505583, (58, 1) = 20.4453125000000, (58, 2) = -577.002806032956, (59, 1) = 20.4531250000000, (59, 2) = -571.503410566831, (60, 1) = 20.4609375000000, (60, 2) = -566.009613811911, (61, 1) = 20.4687500000000, (61, 2) = -560.521407488525, (62, 1) = 20.4765625000000, (62, 2) = -555.038783332590, (63, 1) = 20.4843750000000, (63, 2) = -549.561733095584, (64, 1) = 20.4921875000000, (64, 2) = -544.090248544494, (65, 1) = 20.5000000000000, (65, 2) = -538.624321461797, (66, 1) = 20.5078125000000, (66, 2) = -533.163943645418, (67, 1) = 20.5156250000000, (67, 2) = -527.709106908693, (68, 1) = 20.5234375000000, (68, 2) = -522.259803080339, (69, 1) = 20.5312500000000, (69, 2) = -516.816024004417, (70, 1) = 20.5390625000000, (70, 2) = -511.377761540296, (71, 1) = 20.5468750000000, (71, 2) = -505.945007562620, (72, 1) = 20.5546875000000, (72, 2) = -500.517753961271, (73, 1) = 20.5625000000000, (73, 2) = -495.095992641337, (74, 1) = 20.5703125000000, (74, 2) = -489.679715523079, (75, 1) = 20.5781250000000, (75, 2) = -484.268914541892, (76, 1) = 20.5859375000000, (76, 2) = -478.863581648273, (77, 1) = 20.5937500000000, (77, 2) = -473.463708807788, (78, 1) = 20.6015625000000, (78, 2) = -468.069288001034, (79, 1) = 20.6093750000000, (79, 2) = -462.680311223613, (80, 1) = 20.6171875000000, (80, 2) = -457.296770486090, (81, 1) = 20.6250000000000, (81, 2) = -451.918657813961, (82, 1) = 20.6328125000000, (82, 2) = -446.545965247625, (83, 1) = 20.6406250000000, (83, 2) = -441.178684842342, (84, 1) = 20.6484375000000, (84, 2) = -435.816808668206, (85, 1) = 20.6562500000000, (85, 2) = -430.460328810109, (86, 1) = 20.6640625000000, (86, 2) = -425.109237367706, (87, 1) = 20.6718750000000, (87, 2) = -419.763526455385, (88, 1) = 20.6796875000000, (88, 2) = -414.423188202231, (89, 1) = 20.6875000000000, (89, 2) = -409.088214751994, (90, 1) = 20.6953125000000, (90, 2) = -403.758598263055, (91, 1) = 20.7031250000000, (91, 2) = -398.434330908399, (92, 1) = 20.7109375000000, (92, 2) = -393.115404875574, (93, 1) = 20.7187500000000, (93, 2) = -387.801812366662, (94, 1) = 20.7265625000000, (94, 2) = -382.493545598249, (95, 1) = 20.7343750000000, (95, 2) = -377.190596801385, (96, 1) = 20.7421875000000, (96, 2) = -371.892958221557, (97, 1) = 20.7500000000000, (97, 2) = -366.600622118662, (98, 1) = 20.7578125000000, (98, 2) = -361.313580766961, (99, 1) = 20.7656250000000, (99, 2) = -356.031826455056, (100, 1) = 20.7734375000000, (100, 2) = -350.755351485856, (101, 1) = 20.7812500000000, (101, 2) = -345.484148176547, (102, 1) = 20.7890625000000, (102, 2) = -340.218208858556, (103, 1) = 20.7968750000000, (103, 2) = -334.957525877521, (104, 1) = 20.8046875000000, (104, 2) = -329.702091593261, (105, 1) = 20.8125000000000, (105, 2) = -324.451898379740, (106, 1) = 20.8203125000000, (106, 2) = -319.206938625042, (107, 1) = 20.8281250000000, (107, 2) = -313.967204731331, (108, 1) = 20.8359375000000, (108, 2) = -308.732689114827, (109, 1) = 20.8437500000000, (109, 2) = -303.503384205772, (110, 1) = 20.8515625000000, (110, 2) = -298.279282448394, (111, 1) = 20.8593750000000, (111, 2) = -293.060376300885, (112, 1) = 20.8671875000000, (112, 2) = -287.846658235362, (113, 1) = 20.8750000000000, (113, 2) = -282.638120737838, (114, 1) = 20.8828125000000, (114, 2) = -277.434756308193, (115, 1) = 20.8906250000000, (115, 2) = -272.236557460144, (116, 1) = 20.8984375000000, (116, 2) = -267.043516721208, (117, 1) = 20.9062500000000, (117, 2) = -261.855626632680, (118, 1) = 20.9140625000000, (118, 2) = -256.672879749594, (119, 1) = 20.9218750000000, (119, 2) = -251.495268640697, (120, 1) = 20.9296875000000, (120, 2) = -246.322785888420, (121, 1) = 20.9375000000000, (121, 2) = -241.155424088842, (122, 1) = 20.9453125000000, (122, 2) = -235.993175851664, (123, 1) = 20.9531250000000, (123, 2) = -230.836033800178, (124, 1) = 20.9609375000000, (124, 2) = -225.683990571238, (125, 1) = 20.9687500000000, (125, 2) = -220.537038815227, (126, 1) = 20.9765625000000, (126, 2) = -215.395171196027, (127, 1) = 20.9843750000000, (127, 2) = -210.258380390992, (128, 1) = 20.9921875000000, (128, 2) = -205.126659090915, (129, 1) = 21., (129, 2) = -200.})

#Do same problem with regular Maple BVP solver for comparison:
Sol1:= dsolve({diff(sr(r),r$2) = ODE, sr(a) = P, sr(c) = Q}, numeric, abserr=1e-7):

#difference of two solution methods:
Err:= Vector(2^N+1, i-> eval(sr(r), Sol1(FDM_Sol[i,1])) - FDM_Sol[i,2]):

#maximum absolute error:
LinearAlgebra:-Norm(Err, infinity);

HFloat(5.695120535165188e-7)

 


 

Download BVP_FDM.mw

This is how I like to do a finite product (and it works just as well for add and seq as for mul) that uses all indices except for a specified one. The specified one in this example is i.

i:= 3:  n:= 5:
mul(1/(x[i]-x[j]), j= [$1..i-1, $i+1..n]);

This works even when i is 1 or n.

Another version:

i:= 3:  n:= 5:
mul(1/(x[i]-x[j]), j= {$1..n} minus {i});

 

Your ode is a quadratic expression in diff(sr(r), r$2)---it's highest-order derivative---and thus it's really two odes, one for each solution of the quadratic. Once you isolate these two branches, it's trivial and very fast to solve each with Maple's BVP solver available through dsolve(..., numeric).

Your solution is sufficiently accurate (shown below), although your technique is doomed to be slow due to symbolic expression swell.
 

restart
:

sy := 2400.:
Hp := 0.1e9:
P := -900.:
a := 20.:
c := 21.:
N := 6:
s[0] := P:
h := (c-a)/N:
Et := (r*(diff(sr(r), r))-sy)/Hp:
Er := (sy-r*(diff(sr(r), r)))/Hp:
ode := simplify(Er-Et-r*(diff(Et, r))*(1+(diff(Et, r))*r/(2+4*Et))):
#All above code is yours, unaltered.

ODEs:= solve(ode, diff(sr(r),r$2));

(-3.*r*(diff(sr(r), r))-99995200.+2.*(-1.*r^2*(diff(sr(r), r))^2+4800.*r*(diff(sr(r), r))+0.2499999994e16)^(1/2))/r^2, (-3.*r*(diff(sr(r), r))-99995200.-2.*(-1.*r^2*(diff(sr(r), r))^2+4800.*r*(diff(sr(r), r))+0.2499999994e16)^(1/2))/r^2

Sol1:= dsolve(
    {diff(sr(r),r$2) = ODEs[1], sr(a) = P, sr(c) = -200},
    numeric
):

plots:-odeplot(Sol1, gridlines= false);

seq(eval([r, sr(r)], Sol1(a+k/12)), k= 0..12);

[20., HFloat(-899.9999999999995)], [20.08333333, HFloat(-838.0971805211188)], [20.16666667, HFloat(-776.878078709958)], [20.25000000, HFloat(-716.3318146852567)], [20.33333333, HFloat(-656.447706500654)], [20.41666667, HFloat(-597.2152871965709)], [20.50000000, HFloat(-538.6243206036763)], [20.58333333, HFloat(-480.66475368141937)], [20.66666667, HFloat(-423.32673324607947)], [20.75000000, HFloat(-366.6006214884328)], [20.83333333, HFloat(-310.4769500562954)], [20.91666667, HFloat(-254.94643659443858)], [21., HFloat(-199.99999999999994)]

Sol2:= dsolve(
    {diff(sr(r),r$2) = ODEs[2], sr(a) = P, sr(c) = -200},
    numeric
):

plots:-odeplot(Sol2, gridlines= false);

 


 

Download 2branchBVP.mw

Just do

plots:-odeplot(F, [x(t),y(t),z(t)], t= 0..2);

You can change the range 0..2 to any finite real interval.

First 110 111 112 113 114 115 116 Last Page 112 of 395