Advanced Search

mmcdara

6149 Reputation

17 Badges

9 years, 62 days
Contact mmcdara

MaplePrimes Activity


These are replies submitted by mmcdara

@MaPal93 Meaning of AllTotalDe...

mmcdara 6149
November 15 2023
0 0

@MaPal93 

Meaning of AllTotalDegrees:
AllTotalDegrees is a list of triples, if such a triple is [a, b, c]  then at least one of the equations to solve contains the term

​​​​​​(lambda__1)^a *  (lambda__2)^b *  (lambda__3)^c 



In the second script of yours (where you try to solve for the numerators):
It is you who tried to solve for the numerators, 
The lines

Eqs_redox := eval(Eqs_cal, p = .5); 
EqN := numer(Eqs_cal); nops(%); 
length(EqN); 
EqN_redox := eval(numer(Eqs_redox), lambda__1 = lambda__2); nops(%); 
length(EqN_redox)

are yours, not mine, I just added after them

sol_solve := solve(EqN, Vars):



Why did you use  SolveTools:-PolynomialSystem instead of solve?
(If I'm not mistaken it was a suggestion from @acer but I can't remember why he proposed that?)
Look to this new file from the text  Using SolveTools:-PolynomialSystem
141123_Problem_NoCorrelation_mmcdara_2.mw

Last advice: if you really want to understand what SolveTools:-PolynomialSystem is telling you begin to read these help pages

  • RegularChains
  • PolynomialRing
  • ComprehensiveTriangularize
  • PolynomialSystem

and read some papers or books about the subject... or act simpler.

@Carl Love Thanks Carl, but don'...

mmcdara 6149
November 14 2023
0 0

@Carl Love 

Thanks for your feedback and, by the way, but don't feel obliged to spend too much time my code.

BTW, I wonder if I have not already ask a question about seaching all the occurrences of a specific entry in a hierarchical table, but  couldn't put the finger on it (maybe sevenor eight years ago)

@Carl Love Thanks your comments and...

mmcdara 6149
November 14 2023
0 0

@Carl Love 

Thanks your comments and corrections.

I tried writting a procedure that searches within a hierarchical table (a table whose entries can be tables, and so on) for all the indices where a given entry appears.
Here it is.
 

> 

restart
> 

f := proc(t::table, v::numeric)
  local  r, it, j, depth:
  global here:
  r := select((x->is(rhs(x)=v)), [entries(t,'pairs')]);
  if r <> [] then
    if nargs = 2 then
      here  := (`[]`@lhs)~(r):
      depth := []:
    else
      depth := _passed[3]:
      here  := here, map(i -> [depth[], lhs(i)], r);
    end if:
    it := map(j -> if type(t[j], table) then j end if, [indices](t, nolist));
    if it <> [] then
      for j in it do
        depth := [depth[], j]:
        thisproc(eval(t[j]), v, depth);
      end do
    else
      depth := []
    end if;
  else
    here := []
  end if:
end proc:
 
> 

T := table([ 1=23, 5=36, 3=14, 4=51, 10=51 ]);

here := NULL:
f(T, 51):
map(op, [here]);

table( [( 1 ) = 23, ( 3 ) = 14, ( 4 ) = 51, ( 5 ) = 36, ( 10 ) = 51 ] )

  

[[4], [10]]

 (1)
> 

T := table([ 1=23, 5=36, 3=14, 4=51, 6=table([2=0, 4=51]), 10=51, 0=table([1=8, -3=51, a=table([1=23, r=51, s=51])]) ]):

# From CarlLove
(T-> subsindets(subsindets(`<,>`(op(T)), table, thisproc), list, `<,>`))(T);

here := NULL:
f(T, 51):
map(op, [here]);

Vector(1, {(1) = Vector(7, {(1) = 0 = (Vector(3, {(1) = 1 = 8, (2) = -3 = 51, (3) = a = (Vector(3, {(1) = 1 = 23, (2) = r = 51, (3) = s = 51}))})), (2) = 1 = 23, (3) = 3 = 14, (4) = 4 = 51, (5) = 5 = 36, (6) = 6 = (Vector(2, {(1) = 2 = 0, (2) = 4 = 51})), (7) = 10 = 51})})

  

[[4], [10], [0, -3], [0, a, r], [0, a, s], [0, 6, 4]]

 (2)
> 

 

 

Download Recursive_search.mw

 


Explanation (highlighted is a wrong index, should be [6, 4])

here = [[4], [10]], [[0, -3]], [[0, a, r]], [[0, 6, 4]]

means that the target value (51) is an entry of T[4], T[0][-3], T[0][a][r], ...

Note that I used 

[entries(t,'pairs')]

instead of 

op(op(t))

because op(op(t)) is a list when procedure f is entered the first time, but a sequence when it is entered through thisproc.
(this seems weird to me)


If it's not too much trouble, could you help me correct this procedure, and suggest a more effective version?
Or do you think it would be better if I open a specific thread?

@Carl Love You're right, I didn...

mmcdara 6149
November 14 2023
0 0

@Carl Love 

You're right, I didn't completely answered the question.
I edited my comment and completed it to account your remark.

@dharr My using isolate instead of ...

mmcdara 6149
November 10 2023
0 0

@dharr 

My using isolate instead of solve is mainly a matter of laziness: isolate(equation, v) returns an equality that can directly be used in eval,  as solve (eq, v) requires writing v=solve (eq, v) to be used in eval.

@dharr I'm not at ease with eli...

mmcdara 6149
November 10 2023
0 0

@dharr 

I'm not at ease with eliminate and I often have to do several trials to get what I want.

Nice, I vote yup.

What is a(T) ?...

mmcdara 6149
November 10 2023
0 0

What is a(T) ?

@dharr Thanks...

mmcdara 6149
November 09 2023
0 0

@dharr 

Thanks

Final study...

mmcdara 6149
November 09 2023
0 0

@shashi598

I adapted the Metropolis-Hastings' algorithm (developed to sample an unknown random variable) to find the values of alpha such that D(gamma)(xi__1) = 0.
(
More precisely, as the right boundary condition is gamma(xi__1)+xi__1* D(gamma)(xi__1) = 0, the previous condition is equivalent to gamma(xi__1) = 0, and it is equivalent to search for the values of alpha which make this identity true, or for values of alpha which minimize gamma(xi__1)^2.
I made this latter choice.
)

I got only 2 values (global minima of gamma(xi__1)^2) in the range alpha=1..10:

  1. alpha = 3.7026096348364278322
  2. alpha = 4.2985527068692070036

and 4 local minima (that is minima where gamma(xi__1)^2 is small but not null):

  1. alpha = 5.3002873922249993805
  2. alpha = 6.0544400684250000630
  3. alpha = 7.5717744989250001680
  4. alpha = 8.3303294592000002940

The type (local/global) of the minima are verified graphically.

You will note these values differ from those from the paper you cite.

Reply me if you still read this thread ad are interested to my code.

@Ronan  Thank you Ronan...

mmcdara 6149
November 04 2023
0 0

@Ronan 

Thank you Ronan

No link is active...

mmcdara 6149
November 04 2023
0 0

Could you check them?

A matter of basins of attraction...

mmcdara 6149
November 03 2023
0 0

@shashi598 

Explanation:

Is I'm not mistaken you want to find the value alpha such that diff(gamma(xi), xi) = 0 at xi=xi__1.
Which, give, that gamma(xi) + xi * diff(gamma(xi), xi) = 0 at xi=xi__1, is equivalent to find the value alpha such that gamma(xi) = 0 at xi=xi__1.

Here is the plot, with a logarithmic vertical axis, of gamma(xi__1)^2 in the range 1..6:


Then a focus on the range 3.5..4


and another on the range 4..6

Here are the values of alpha that minimize H(xi) = amma(xi__1)^2 when alpha is assumed to move in diffferent ranges:

This last plot completes the previous one by adding, in gray", the basins of attraction of red points, and in green those of the blue points.
Depending on your starting point in the range 4..6, NLPSolve will get anyone of the 5 points.
By default the starting point (if not specified with the option initialpoint) is the mid point of the searching range, thus 5 here.
Starting from 5 will return the first leftmost red point on this figure (that is your alpha2_result point).

Using DirectSearch may enable jumping out of this basin of attraction (the leftmost green rectangle) for it is not very deep.
Unfortunately the algorithm ends hitting the right boundary of the search domain (where the value of gamma(xi__1)^2 is lower than those at the red point, so it indeed minimizes the objective, at least in some way), but it misses the blue points because of their thinness and a bad tunning (I guess so) of the algorithm. 


Duplicate question?...

mmcdara 6149
November 03 2023
0 0



Isn't it this a duplicate of thar previous question?

BTW, why do you persist in doing the same mistake (unknown=f(x) and parameter mu[f])?
Did you read my answer?


More of this your 4 equations are both identical except for the name of the dependent function, so their solutions will be fundamentaly identical it it is stupid to dsolve 4 equations.
The generic form of these equations is given by coeff(HPMeq, p, 1).
Before doing complex things  ask yourself if this equation can be solved formally; 

# example
test := coeff(HPMEq, p, 0);

# A solution is returned here
Without_bc := dsolve(test);   

# But not here
With_bc := dsolve({test, F[0](1) = 1, F[0](-1) = -1, (D(F[0]))(1) = 0, (D(F[0]))(-1)})

So trying to take the rhs of this last ouput generates the error message you got.

The only way is to dsolve numerically

{test, F[0](1) = 1, F[0](-1) = -1, (D(F[0]))(1) = 0, (D(F[0]))(-1)}

which implies to replace the many parameters by their numerical values.

Make an effort to understand all this and to submit us something that is worth us spending time answering, otherwise it's all a waste of  time for the people here.

odeplot F not f !!!!!!...

mmcdara 6149
November 03 2023
0 0

@Saha 

 

> 

restart;
> 

with(plots):

 

HERE IS A CORRECT WAY TO CODE THE THINGS

Note that I treated ONLY the first ode.

With a little effort you should be capable to generalize this to the two other odes, so stop
doing again and again the same errors and try to learn something.

 

> 

eq1 := mu[hnf]*rho[f]*(diff(F(x), x, x, x, x))/(rho[hnf]*mu[f])+3*alpha*(diff(F(x), x, x))+alpha*eta*(diff(F(x), x, x, x))-2*R*F(x)*(diff(F(x), x, x, x))-rho[f]*M*(diff(F(x), x, x))/rho[hnf] = 0;

mu[hnf]*rho[f]*(diff(diff(diff(diff(F(x), x), x), x), x))/(rho[hnf]*mu[f])+3*alpha*(diff(diff(F(x), x), x))+alpha*eta*(diff(diff(diff(F(x), x), x), x))-2*R*F(x)*(diff(diff(diff(F(x), x), x), x))-rho[f]*M*(diff(diff(F(x), x), x))/rho[hnf] = 0

 (1)
> 

a1 := { phi[1] = 0.1e-1, phi[2] = 0.1e-1, R = 1.0, M = 1, alpha = -1, eta = 1 }:
a2 := { phi[1] = 0.1e-1, phi[2] = 0.1e-1, R = 1.0, M = 1, alpha = -.5, eta = 1 }:
a3 := { phi[1] = 0.1e-1, phi[2] = 0.1e-1, R = 1.0, M = 1, alpha = .5, eta = 1 }:
a4 := { phi[1] = 0.1e-1, phi[2] = 0.1e-1, R = 1.0, M = 1, alpha = 1, eta = 1 }:

AllA := [a1, a2, a3, a4]:

print~(AllA):

{M = 1, R = 1.0, alpha = -1, eta = 1, phi[1] = 0.1e-1, phi[2] = 0.1e-1}

  

{M = 1, R = 1.0, alpha = -.5, eta = 1, phi[1] = 0.1e-1, phi[2] = 0.1e-1}

  

{M = 1, R = 1.0, alpha = .5, eta = 1, phi[1] = 0.1e-1, phi[2] = 0.1e-1}

  

{M = 1, R = 1.0, alpha = 1, eta = 1, phi[1] = 0.1e-1, phi[2] = 0.1e-1}

 (2)
> 

data := {
  rho[f]   = 997.1,
  rho[s1]  = 6320,
  rho[s2]  = 3970,
  c[p][f]  = 4180,
  c[p][s1] = 531.5,
  c[p][s2] = 765,
  k[s1]    = 76.5,
  k[s2]    = 40,
  k[f]     = .613,
  mu[f]   = 1,     # You forgot assigning it, this is my arbitrary choice
  k[nf]   = 1,     # You forgot assigning it, this is my arbitrary choice
  k[s1]   = 1,     # You forgot assigning it, this is my arbitrary choice
  k[s2]   = 1,     # You forgot assigning it, this is my arbitrary choice
  phi[1]  = 0.1,   # You forgot assigning it, this is my arbitrary choice
  phi[2]  = 0.1    # You forgot assigning it, this is my arbitrary choice
}:

AllData := {
  data[]
  , k[nf]    = eval( k[f]*(k[s1]+2*k[f]-2*phi[1]*(k[f]-k[s1]))/(k[s1]+2*k[f]+phi[1]*(k[f]-k[s1])), data)
  , k[hnf]   = simplify(eval( k[nf]*(k[s2]+2*k[nf]-2*phi[2]*(k[nf]-k[s2]))/(k[s2]+2*k[nf]+phi[2]*(k[nf]-k[s2])), data))
  , mu[hnf]  = eval( mu[f]/((1-phi[1])^2.5*(1-phi[2])^2.5), data)
  , rho[hnf] = eval( (1-phi[2])((1-phi[1])*rho[f]+phi[1]*rho[s1])+phi[2]*rho[s2], data)
}:

# Check:

print~(AllData):
 

k[f] = .613

  

k[hnf] = 1.000000000

  

k[nf] = 1

  

k[nf] = .6455375120

  

k[s1] = 1

  

k[s1] = 76.5

  

k[s2] = 1

  

k[s2] = 40

  

mu[f] = 1

  

mu[hnf] = 1.693508780

  

phi[1] = .1

  

phi[2] = .1

  

rho[f] = 997.1

  

rho[hnf] = 397.9

  

rho[s1] = 6320

  

rho[s2] = 3970

  

c[p][f] = 4180

  

c[p][s1] = 531.5

  

c[p][s2] = 765

 (3)
> 

AllB := [ seq(eval(eq1, AllA[i] union AllData), i=1..numelems(AllA)) ]:

print~(AllB):

4.243773825*(diff(diff(diff(diff(F(x), x), x), x), x))-5.505906007*(diff(diff(F(x), x), x))-(diff(diff(diff(F(x), x), x), x))-2.0*F(x)*(diff(diff(diff(F(x), x), x), x)) = 0

  

4.243773825*(diff(diff(diff(diff(F(x), x), x), x), x))-4.005906007*(diff(diff(F(x), x), x))-.5*(diff(diff(diff(F(x), x), x), x))-2.0*F(x)*(diff(diff(diff(F(x), x), x), x)) = 0

  

4.243773825*(diff(diff(diff(diff(F(x), x), x), x), x))-1.005906007*(diff(diff(F(x), x), x))+.5*(diff(diff(diff(F(x), x), x), x))-2.0*F(x)*(diff(diff(diff(F(x), x), x), x)) = 0

  

4.243773825*(diff(diff(diff(diff(F(x), x), x), x), x))+.494093993*(diff(diff(F(x), x), x))+diff(diff(diff(F(x), x), x), x)-2.0*F(x)*(diff(diff(diff(F(x), x), x), x)) = 0

 (4)
> 

bcs := F(-1) = -1, F(1) = 1, (D(F))(-1) = 0, (D(F))(1) = 0;

F(-1) = -1, F(1) = 1, (D(F))(-1) = 0, (D(F))(1) = 0

 (5)
> 

AllC := [ seq( dsolve({AllB[i], bcs}, numeric), i=1..numelems(AllA)) ];

AllC[1](0);

colors  := [red, gold, green, blue]:
legends := [ seq(cat("data set a", i), i=1..numelems(AllA)) ];

display(
  seq(
    odeplot(AllC[i], [x, diff(F(x), x)], -1..1, color=colors[i], legend=legends[i])
    , i=1..numelems(AllA)
  )
)

[proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(10, {(1) = -1.0, (2) = -.8345204915051194, (3) = -.613691680718111, (4) = -.3419593760168745, (5) = -0.30799685620354136e-1, (6) = .261867322088955, (7) = .49531217387834137, (8) = .6834721467651494, (9) = .8471616133630866, (10) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(10, 4, {(1, 1) = -1.0, (1, 2) = .0, (1, 3) = 3.1636830295853944, (1, 4) = -4.607032094044634, (2, 1) = -.9600083071821683, (2, 2) = .4641435042824518, (2, 3) = 2.4670081858897706, (2, 4) = -3.8449519456918044, (3, 1) = -.8038896314203525, (3, 2) = .9217743985030411, (3, 3) = 1.7046802961225183, (3, 4) = -3.1099803799278667, (4, 1) = -.5003436588990865, (4, 2) = 1.2777956453763997, (4, 3) = .93926555275606, (4, 4) = -2.5886165747900933, (5, 1) = -0.6991456323748237e-1, (5, 2) = 1.4489743289330363, (5, 3) = .1685876354270685, (5, 4) = -2.442192850610124, (6, 1) = .35099468269583234, (6, 2) = 1.3910256110120378, (6, 3) = -.5792822530856752, (6, 4) = -2.745091306679072, (7, 1) = .653851403072871, (7, 2) = 1.1762711862773043, (7, 3) = -1.2844409058641224, (7, 4) = -3.3599344632577943, (8, 1) = .8485116918581102, (8, 2) = .8707876943971833, (8, 3) = -1.9881043917473578, (8, 4) = -4.17394509227341, (9, 1) = .9611973502450049, (9, 2) = .4852860212631082, (9, 3) = -2.7494939818683, (9, 4) = -5.1808107093802835, (10, 1) = 1.0, (10, 2) = .0, (10, 3) = -3.632609179737446, (10, 4) = -6.429500616264703}, datatype = float[8], order = C_order); YP := Matrix(10, 4, {(1, 1) = .0, (1, 2) = 3.1636830295853944, (1, 3) = -4.607032094044634, (1, 4) = 5.190185528061891, (2, 1) = .4641435042824518, (2, 2) = 2.4670081858897706, (2, 3) = -3.8449519456918044, (2, 4) = 4.034271280109428, (3, 1) = .9217743985030411, (3, 2) = 1.7046802961225183, (3, 3) = -3.1099803799278667, (3, 4) = 2.657066924436565, (4, 1) = 1.2777956453763997, (4, 2) = .93926555275606, (4, 3) = -2.5886165747900933, (4, 4) = 1.2190298693244515, (5, 1) = 1.4489743289330363, (5, 2) = .1685876354270685, (5, 3) = -2.442192850610124, (5, 4) = -.2762813315156067, (6, 1) = 1.3910256110120378, (6, 2) = -.5792822530856752, (6, 3) = -2.745091306679072, (6, 4) = -1.8524997259979574, (7, 1) = 1.1762711862773043, (7, 2) = -1.2844409058641224, (7, 3) = -3.3599344632577943, (7, 4) = -3.493527624215365, (8, 1) = .8707876943971833, (8, 2) = -1.9881043917473578, (8, 3) = -4.17394509227341, (8, 4) = -5.232027988502383, (9, 1) = .4852860212631082, (9, 2) = -2.7494939818683, (9, 3) = -5.1808107093802835, (9, 4) = -7.134882875694695, (10, 1) = .0, (10, 2) = -3.632609179737446, (10, 3) = -6.429500616264703, (10, 4) = -9.25810567969981}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(10, {(1) = -1.0, (2) = -.8345204915051194, (3) = -.613691680718111, (4) = -.3419593760168745, (5) = -0.30799685620354136e-1, (6) = .261867322088955, (7) = .49531217387834137, (8) = .6834721467651494, (9) = .8471616133630866, (10) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(10, 4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = -0.320965229302303e-6, (1, 4) = 0.20625435773980762e-6, (2, 1) = -0.21662828660741597e-8, (2, 2) = -0.5321047471110476e-7, (2, 3) = -0.29748745205442183e-6, (2, 4) = 0.14474251093437116e-6, (3, 1) = -0.12802224185955755e-7, (3, 2) = -0.12880359685429147e-6, (3, 3) = -0.2754609284674051e-6, (3, 4) = 0.422248919810272e-7, (4, 1) = -0.2488039788282393e-7, (4, 2) = -0.24645040885322146e-6, (4, 3) = -0.2232055894767619e-6, (4, 4) = -0.16612887001043793e-6, (5, 1) = -0.30207191488881113e-8, (5, 2) = -0.31616645991050967e-6, (5, 3) = -0.8442964129401751e-7, (5, 4) = -0.33082012445115147e-6, (6, 1) = 0.17692274965987335e-7, (6, 2) = -0.21929135648323933e-6, (6, 3) = 0.7591113124113682e-7, (6, 4) = -0.16518221277806597e-6, (7, 1) = 0.8997213689790162e-8, (7, 2) = -0.14049504318983955e-6, (7, 3) = 0.15011911446779624e-6, (7, 4) = -0.3775879681945588e-8, (8, 1) = 0.4950006422490701e-9, (8, 2) = -0.8694952700639968e-7, (8, 3) = 0.19486283016094694e-6, (8, 4) = 0.10646699460348378e-6, (9, 1) = -0.233653484917027e-8, (9, 2) = -0.4145693924776069e-7, (9, 3) = 0.23275364039844293e-6, (9, 4) = 0.19144199127325873e-6, (10, 1) = .0, (10, 2) = .0, (10, 3) = 0.26202545849335726e-6, (10, 4) = 0.24531565723179697e-6}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[10] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(3.3082012445115147e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [4, 10, [F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[10] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[10] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(4, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(10, 4, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(4, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(10, 4, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)]'[i] = yout[i], i = 1 .. 4)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[10] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(3.3082012445115147e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [4, 10, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[10] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[10] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(10, 4, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(4, {(1) = 0., (2) = 0., (3) = 0., (4) = 0.}); `dsolve/numeric/hermite`(10, 4, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 4)] end proc, (2) = Array(0..0, {}), (3) = [x, F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)]'[i] = res[i+1], i = 1 .. 4)] catch: error end try end proc, proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(10, {(1) = -1.0, (2) = -.8230223002196657, (3) = -.6120604124702849, (4) = -.365734113634546, (5) = -0.8521083706438085e-1, (6) = .19807440722060754, (7) = .4441585008746424, (8) = .6506759302644971, (9) = .8310728926272699, (10) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(10, 4, {(1, 1) = -1.0, (1, 2) = .0, (1, 3) = 3.2024054605007, (1, 4) = -4.561772646438333, (2, 1) = -.9538821065519519, (2, 2) = .499354216305126, (2, 3) = 2.462188482329671, (2, 4) = -3.834015030700587, (3, 1) = -.7994699363427611, (3, 2) = .9386310747679936, (3, 3) = 1.7246690988739413, (3, 4) = -3.1980477350676435, (4, 1) = -.5235639084761247, (4, 2) = 1.2718031406443107, (4, 3) = .9992635168967455, (4, 4) = -2.739308864691832, (5, 1) = -.13728232826492662, (5, 2) = 1.447997749307266, (5, 3) = .26608580079925903, (5, 4) = -2.5429470673842762, (6, 1) = .27392352756220795, (6, 2) = 1.4206602755546929, (6, 3) = -.46549948973869354, (6, 4) = -2.6793461768366402, (7, 1) = .6025721387687291, (7, 2) = 1.2215183660280955, (7, 3) = -1.1701367490289056, (7, 4) = -3.0985615431983247, (8, 1) = .8251391118646613, (8, 2) = .9098931672344264, (8, 3) = -1.8688868827206906, (8, 4) = -3.7129945299859637, (9, 1) = .9550662052310327, (9, 2) = .5084512958428912, (9, 3) = -2.605114303252651, (9, 4) = -4.49032135449593, (10, 1) = 1.0, (10, 2) = .0, (10, 3) = -3.4418755414950137, (10, 4) = -5.4574567197898505}, datatype = float[8], order = C_order); YP := Matrix(10, 4, {(1, 1) = .0, (1, 2) = 3.2024054605007, (1, 3) = -4.561772646438333, (1, 4) = 4.635306934795669, (2, 1) = .499354216305126, (2, 2) = 2.462188482329671, (2, 3) = -3.834015030700587, (2, 4) = 3.596017463172271, (3, 1) = .9386310747679936, (3, 2) = 1.7246690988739413, (3, 3) = -3.1980477350676435, (3, 4) = 2.4561451443713196, (4, 1) = 1.2718031406443107, (4, 2) = .9992635168967455, (4, 3) = -2.739308864691832, (4, 4) = 1.2964187138341203, (5, 1) = 1.447997749307266, (5, 2) = .26608580079925903, (5, 3) = -2.5429470673842762, (5, 4) = .11608643215970291, (6, 1) = 1.4206602755546929, (6, 2) = -.46549948973869354, (6, 3) = -2.6793461768366402, (6, 4) = -1.1009757814453724, (7, 1) = 1.2215183660280955, (7, 2) = -1.1701367490289056, (7, 3) = -3.0985615431983247, (7, 4) = -2.349548474329587, (8, 1) = .9098931672344264, (8, 2) = -1.8688868827206906, (8, 3) = -3.7129945299859637, (8, 4) = -3.645471486112765, (9, 1) = .5084512958428912, (9, 2) = -2.605114303252651, (9, 3) = -4.49032135449593, (9, 4) = -5.009247180175052, (10, 1) = .0, (10, 2) = -3.4418755414950137, (10, 3) = -5.4574567197898505, (10, 4) = -6.463933479418131}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(10, {(1) = -1.0, (2) = -.8230223002196657, (3) = -.6120604124702849, (4) = -.365734113634546, (5) = -0.8521083706438085e-1, (6) = .19807440722060754, (7) = .4441585008746424, (8) = .6506759302644971, (9) = .8310728926272699, (10) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(10, 4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = -0.21795585591266545e-6, (1, 4) = 0.15859585924547324e-6, (2, 1) = -0.7382788075766052e-9, (2, 2) = -0.3484970048445586e-7, (2, 3) = -0.20310495291492214e-6, (2, 4) = 0.13757794743245046e-6, (3, 1) = -0.4824796191312689e-8, (3, 2) = -0.786371881262705e-7, (3, 3) = -0.17951856691697007e-6, (3, 4) = 0.9262522073134566e-7, (4, 1) = -0.9653802963899059e-8, (4, 2) = -0.13506031575634385e-6, (4, 3) = -0.13410901880448737e-6, (4, 4) = 0.10579088076766957e-7, (5, 1) = -0.5168842837487e-8, (5, 2) = -0.1808930499416432e-6, (5, 3) = -0.55580325882813683e-7, (5, 4) = -0.6730545231057284e-7, (6, 1) = 0.57062913913256185e-8, (6, 2) = -0.15787889363918775e-6, (6, 3) = 0.3988819587012066e-7, (6, 4) = -0.3456323983230098e-7, (7, 1) = 0.36524231994704153e-8, (7, 2) = -0.10923881948943386e-6, (7, 3) = 0.10645399214968022e-6, (7, 4) = 0.5146402980834153e-7, (8, 1) = -0.5261386113252735e-9, (8, 2) = -0.6721228628347843e-7, (8, 3) = 0.1509773553809045e-6, (8, 4) = 0.12360533872225815e-6, (9, 1) = -0.1690561041667948e-8, (9, 2) = -0.30794920851512984e-7, (9, 3) = 0.1838203466569299e-6, (9, 4) = 0.17408963220781308e-6, (10, 1) = .0, (10, 2) = .0, (10, 3) = 0.20174679478711903e-6, (10, 4) = 0.1916208597032606e-6}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[10] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(2.1795585591266545e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [4, 10, [F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[10] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[10] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(4, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(10, 4, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(4, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(10, 4, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)]'[i] = yout[i], i = 1 .. 4)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[10] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(2.1795585591266545e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [4, 10, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[10] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[10] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(10, 4, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(4, {(1) = 0., (2) = 0., (3) = 0., (4) = 0.}); `dsolve/numeric/hermite`(10, 4, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 4)] end proc, (2) = Array(0..0, {}), (3) = [x, F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)]'[i] = res[i+1], i = 1 .. 4)] catch: error end try end proc, proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(11, {(1) = -1.0, (2) = -.8349702453362813, (3) = -.6629989401386771, (4) = -.47695397178592847, (5) = -.271772047084955, (6) = -0.4754132124094439e-1, (7) = .1858007458755159, (8) = .41429059327121887, (9) = .6299012307002473, (10) = .827687111900952, (11) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(11, 4, {(1, 1) = -1.0, (1, 2) = .0, (1, 3) = 3.305926822008611, (1, 4) = -4.523672488554754, (2, 1) = -.9582671849373837, (2, 2) = .4864591352968459, (2, 3) = 2.6039107003814994, (2, 4) = -3.9993592087600627, (3, 1) = -.8393949741649898, (3, 2) = .8774432042162891, (3, 3) = 1.9559445568175164, (3, 4) = -3.5537943498965046, (4, 1) = -.6460071110487181, (4, 2) = 1.1821174091381894, (4, 3) = 1.3306943429352183, (4, 4) = -3.187151863442845, (5, 1) = -.3799225547985622, (5, 2) = 1.3902207781570464, (5, 3) = .707183835292867, (5, 4) = -2.911995514638119, (6, 1) = -0.5579044150654242e-1, (6, 2) = 1.4772264934892843, (6, 3) = 0.7486120014579348e-1, (6, 4) = -2.7510082082307656, (7, 1) = .28515756576090134, (7, 2) = 1.4203581008545112, (7, 3) = -.5612480171204297, (7, 4) = -2.724169791644911, (8, 1) = .5895951934544567, (8, 2) = 1.220374178925209, (8, 3) = -1.1932223381511113, (8, 4) = -2.8287988271332107, (9, 1) = .8201843814448493, (9, 2) = .8959047730799011, (9, 3) = -1.8242537511744776, (9, 4) = -3.0431398450851828, (10, 1) = .9576901660006452, (10, 2) = .4738110215978004, (10, 3) = -2.4535632213564074, (10, 4) = -3.3350595641956406, (11, 1) = 1.0, (11, 2) = .0, (11, 3) = -3.0550877141299635, (11, 4) = -3.656004771833288}, datatype = float[8], order = C_order); YP := Matrix(11, 4, {(1, 1) = .0, (1, 2) = 3.305926822008611, (1, 3) = -4.523672488554754, (1, 4) = 3.448495012654861, (2, 1) = .4864591352968459, (2, 2) = 2.6039107003814994, (2, 3) = -3.9993592087600627, (2, 4) = 2.8945648159393023, (3, 1) = .8774432042162891, (3, 2) = 1.9559445568175164, (3, 3) = -3.5537943498965046, (3, 4) = 2.28816807574519, (4, 1) = 1.1821174091381894, (4, 2) = 1.3306943429352183, (4, 3) = -3.187151863442845, (4, 4) = 1.6612513274817053, (5, 1) = 1.3902207781570464, (5, 2) = .707183835292867, (5, 3) = -2.911995514638119, (5, 4) = 1.0321058465598076, (6, 1) = 1.4772264934892843, (6, 2) = 0.7486120014579348e-1, (6, 3) = -2.7510082082307656, (6, 4) = .414199114412958, (7, 1) = 1.4203581008545112, (7, 2) = -.5612480171204297, (7, 3) = -2.724169791644911, (7, 4) = -.17816998271002774, (8, 1) = 1.220374178925209, (8, 2) = -1.1932223381511113, (8, 3) = -2.8287988271332107, (8, 4) = -.7355628778168632, (9, 1) = .8959047730799011, (9, 2) = -1.8242537511744776, (9, 3) = -3.0431398450851828, (9, 4) = -1.250144245620543, (10, 1) = .4738110215978004, (10, 2) = -2.4535632213564074, (10, 3) = -3.3350595641956406, (10, 4) = -1.6938771934039125, (11, 1) = .0, (11, 2) = -3.0550877141299635, (11, 3) = -3.656004771833288, (11, 4) = -2.0163982799684574}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(11, {(1) = -1.0, (2) = -.8349702453362813, (3) = -.6629989401386771, (4) = -.47695397178592847, (5) = -.271772047084955, (6) = -0.4754132124094439e-1, (7) = .1858007458755159, (8) = .41429059327121887, (9) = .6299012307002473, (10) = .827687111900952, (11) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(11, 4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = -0.9790569767920404e-7, (1, 4) = 0.112448305705116e-6, (2, 1) = 0.45844983361163896e-9, (2, 2) = -0.1094032518009447e-7, (2, 3) = -0.8955075685573662e-7, (2, 4) = 0.1191242510363102e-6, (3, 1) = 0.8009839376081181e-9, (3, 2) = -0.2360070345808789e-7, (3, 3) = -0.7338664365543241e-7, (3, 4) = 0.11390192178403137e-6, (4, 1) = 0.6876293051178672e-9, (4, 2) = -0.35999957113990896e-7, (4, 3) = -0.5141334884606263e-7, (4, 4) = 0.10156126373300032e-6, (5, 1) = 0.7037244555939593e-9, (5, 2) = -0.4632169913045516e-7, (5, 3) = -0.2399556636436062e-7, (5, 4) = 0.8670442182328977e-7, (6, 1) = 0.17494118211695317e-8, (6, 2) = -0.51600905679493496e-7, (6, 3) = 0.7897919546512352e-8, (6, 4) = 0.7578789917080854e-7, (7, 1) = 0.3147374590803252e-8, (7, 2) = -0.4857425579192532e-7, (7, 3) = 0.3987221642558589e-7, (7, 4) = 0.7437579754899402e-7, (8, 1) = 0.33146152582413634e-8, (8, 2) = -0.3828130722331125e-7, (8, 3) = 0.6626767042182201e-7, (8, 4) = 0.7959026322089725e-7, (9, 1) = 0.21983578800699434e-8, (9, 2) = -0.25371684974673274e-7, (9, 3) = 0.8406273273179104e-7, (9, 4) = 0.8317133836751885e-7, (10, 1) = 0.5437855531352124e-9, (10, 2) = -0.12802270154238087e-7, (10, 3) = 0.9379983254400154e-7, (10, 4) = 0.8062674681017602e-7, (11, 1) = .0, (11, 2) = .0, (11, 3) = 0.10052401942370861e-6, (11, 4) = 0.7587793139634221e-7}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[11] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.191242510363102e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [4, 11, [F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[11] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[11] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(4, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(11, 4, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(4, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(11, 4, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)]'[i] = yout[i], i = 1 .. 4)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[11] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.191242510363102e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [4, 11, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[11] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[11] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(11, 4, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(4, {(1) = 0., (2) = 0., (3) = 0., (4) = 0.}); `dsolve/numeric/hermite`(11, 4, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 4)] end proc, (2) = Array(0..0, {}), (3) = [x, F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)]'[i] = res[i+1], i = 1 .. 4)] catch: error end try end proc, proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(13, {(1) = -1.0, (2) = -.8648048849604972, (3) = -.7258144766126813, (4) = -.5806548650009224, (5) = -.4251954881870461, (6) = -.25830200663898156, (7) = -0.7906320968469074e-1, (8) = .10896859691459088, (9) = .30218845783835674, (10) = .4950993608940931, (11) = .6842045337170743, (12) = .8563830920544138, (13) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(13, 4, {(1, 1) = -1.0, (1, 2) = .0, (1, 3) = 3.373722494305436, (1, 4) = -4.540543484968664, (2, 1) = -.9709994051712143, (2, 2) = .4157529846211831, (2, 3) = 2.784924226032603, (2, 4) = -4.175384508203814, (3, 1) = -.8881431843504255, (3, 2) = .7636161458832533, (3, 3) = 2.2284237467793426, (3, 4) = -3.8397480759739473, (4, 1) = -.7557361304421567, (4, 2) = 1.047742309074612, (4, 3) = 1.6935372508295938, (4, 4) = -3.5384993911508076, (5, 1) = -.57456181187085, (5, 2) = 1.2693883369188317, (5, 3) = 1.1648387824551296, (5, 4) = -3.2729827638861617, (6, 1) = -.3489758737997563, (6, 2) = 1.4193148302520373, (6, 3) = .6380114908991168, (6, 4) = -3.050802164777067, (7, 1) = -0.8721290290439132e-1, (7, 2) = 1.4856754446728115, (7, 3) = .10760443033540439, (7, 4) = -2.8784616899570166, (8, 1) = .1908914101879065, (8, 2) = 1.4557957427238546, (8, 3) = -.4218076470447148, (8, 4) = -2.7630465319931523, (9, 1) = .4610068956882461, (9, 2) = 1.3231672424400462, (9, 3) = -.9491453862557799, (9, 4) = -2.704760635382748, (10, 1) = .6953692349713777, (10, 2) = 1.0898540289713925, (10, 3) = -1.4695029587825694, (10, 4) = -2.697678065088815, (11, 1) = .8721442840050412, (11, 2) = .7635868911696608, (11, 3) = -1.9821500585914997, (11, 4) = -2.7295310300633284, (12, 1) = .9719047540926605, (12, 2) = .38161022893967533, (12, 3) = -2.456305521844615, (12, 4) = -2.780622799851062, (13, 1) = 1.0, (13, 2) = .0, (13, 3) = -2.8591443331860087, (13, 4) = -2.829325237274913}, datatype = float[8], order = C_order); YP := Matrix(13, 4, {(1, 1) = .0, (1, 2) = 3.373722494305436, (1, 3) = -4.540543484968664, (1, 4) = 2.8169961287747656, (2, 1) = .4157529846211831, (2, 2) = 2.784924226032603, (2, 3) = -4.175384508203814, (2, 4) = 2.5703447861797697, (3, 1) = .7636161458832533, (3, 2) = 2.2284237467793426, (3, 3) = -3.8397480759739473, (3, 4) = 2.252520951780353, (4, 1) = 1.047742309074612, (4, 2) = 1.6935372508295938, (4, 3) = -3.5384993911508076, (4, 4) = 1.8969145895692585, (5, 1) = 1.2693883369188317, (5, 2) = 1.1648387824551296, (5, 3) = -3.2729827638861617, (5, 4) = 1.5218776963785599, (6, 1) = 1.4193148302520373, (6, 2) = .6380114908991168, (6, 3) = -3.050802164777067, (6, 4) = 1.1463563853206438, (7, 1) = 1.4856754446728115, (7, 2) = .10760443033540439, (7, 3) = -2.8784616899570166, (7, 4) = .784059925028772, (8, 1) = 1.4557957427238546, (8, 2) = -.4218076470447148, (8, 3) = -2.7630465319931523, (8, 4) = .4516205476114067, (9, 1) = 1.3231672424400462, (9, 2) = -.9491453862557799, (9, 3) = -2.704760635382748, (9, 4) = .16021142715357983, (10, 1) = 1.0898540289713925, (10, 2) = -1.4695029587825694, (10, 3) = -2.697678065088815, (10, 4) = -0.7729300109926096e-1, (11, 1) = .7635868911696608, (11, 2) = -1.9821500585914997, (11, 3) = -2.7295310300633284, (11, 4) = -.24793741323626228, (12, 1) = .38161022893967533, (12, 2) = -2.456305521844615, (12, 3) = -2.780622799851062, (12, 4) = -.33242403861147385, (13, 1) = .0, (13, 2) = -2.8591443331860087, (13, 3) = -2.829325237274913, (13, 4) = -.3338159043213642}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(13, {(1) = -1.0, (2) = -.8648048849604972, (3) = -.7258144766126813, (4) = -.5806548650009224, (5) = -.4251954881870461, (6) = -.25830200663898156, (7) = -0.7906320968469074e-1, (8) = .10896859691459088, (9) = .30218845783835674, (10) = .4950993608940931, (11) = .6842045337170743, (12) = .8563830920544138, (13) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(13, 4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = -0.3415082014380156e-7, (1, 4) = 0.48382845169449726e-7, (2, 1) = 0.27985988830463634e-9, (2, 2) = -0.29747826355335317e-8, (2, 3) = -0.3095389907332211e-7, (2, 4) = 0.5147221467081313e-7, (3, 1) = 0.5634000429010082e-9, (3, 2) = -0.6211191988774322e-8, (3, 3) = -0.25930801472913458e-7, (3, 4) = 0.5154159579146692e-7, (4, 1) = 0.7626915129413819e-9, (4, 2) = -0.9346779650031032e-8, (4, 3) = -0.19466852860884187e-7, (4, 4) = 0.49264545198076486e-7, (5, 1) = 0.9192783248897028e-9, (5, 2) = -0.12057649486810488e-7, (5, 3) = -0.11869291988131077e-7, (5, 4) = 0.4550569117408885e-7, (6, 1) = 0.11235020840583405e-8, (6, 2) = -0.14047846787318742e-7, (6, 3) = -0.3354285736686311e-8, (6, 4) = 0.4099263112608454e-7, (7, 1) = 0.14812479427117337e-8, (7, 2) = -0.1499312754963664e-7, (7, 3) = 0.5692827186509835e-8, (7, 4) = 0.3658035870208618e-7, (8, 1) = 0.18709366046929603e-8, (8, 2) = -0.1462195147457879e-7, (8, 3) = 0.144873707517166e-7, (8, 4) = 0.32954590048008254e-7, (9, 1) = 0.20523163456542436e-8, (9, 2) = -0.12985898887961813e-7, (9, 3) = 0.22089109706539154e-7, (9, 4) = 0.3014074072119169e-7, (10, 1) = 0.1742172910215283e-8, (10, 2) = -0.10486608719607884e-7, (10, 3) = 0.2771055675320638e-7, (10, 4) = 0.27458856045370038e-7, (11, 1) = 0.8783277962319171e-9, (11, 2) = -0.765167083971309e-8, (11, 3) = 0.3113401515282755e-7, (11, 4) = 0.2425696806134047e-7, (12, 1) = 0.10072978358224701e-9, (12, 2) = -0.41105690627636155e-8, (12, 3) = 0.33450107039721626e-7, (12, 4) = 0.21040741295422602e-7, (13, 1) = .0, (13, 2) = .0, (13, 3) = 0.35515395588326785e-7, (13, 4) = 0.18772939557195428e-7}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[13] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(5.154159579146692e-8) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [4, 13, [F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[13] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[13] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(4, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(13, 4, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(4, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(13, 4, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)]'[i] = yout[i], i = 1 .. 4)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[13] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(5.154159579146692e-8) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [4, 13, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[13] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[13] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(13, 4, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(4, {(1) = 0., (2) = 0., (3) = 0., (4) = 0.}); `dsolve/numeric/hermite`(13, 4, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 4)] end proc, (2) = Array(0..0, {}), (3) = [x, F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)]'[i] = res[i+1], i = 1 .. 4)] catch: error end try end proc]

  

[x = 0., F(x) = HFloat(-0.02521855623397279), diff(F(x), x) = HFloat(1.4530068811004437), diff(diff(F(x), x), x) = HFloat(0.09321390760388122), diff(diff(diff(F(x), x), x), x) = HFloat(-2.453032561333428)]

  

["data set a1", "data set a2", "data set a3", "data set a4"]

  

 
> 

 

 

Download eq1_only.mw

@shashi598 Let's try to synchro...

mmcdara 6149
November 02 2023
0 0

@shashi598 

Let's try to synchronize our exchanges.

What I propose is for you to read my previous reply (that I updated meanwhile), and for me to wait for your return.
Maybe, or maybe not, the procedure LaneEmdenSolution  my last reply contains will give you some answers.

If not, I propose you to explain again what you want to achieve . For instance:

  • Is it enough for you to compute the solution (n=1) from 0 to Pi?
    If not, look below to the attached file.
     
  • Is the problem for n <> 1 of the same nature (a condition, or boundary condition) at some specific value of xi?
    The LaneEmdenSolution procedure was written assuming this, but I can be cometely out of the way.

See you soon

BTW, here is a procedure that builds the numerical solution of the Lane-Emden equafion on any range 0..some_value:

> 

restart


Numerical solution of the Lane-Emden equation on arbitrary ranges

Procedure:

   1  Solve numerically the Emden equation  (result = EmdenN)

   2  For the Lane-Emden equation

        2.1   Use the known option option to instruct dsolve that it has to use the EmdenN solution
        2.2   Solve numerically a boundary value problem for the Lane-Emden equation (result = LaneEmden_bc)
                from 0 to a given abscissa xi = BC_IC_transition_point with conditions
       D(gamma)(0) = 0,
       gamma(BC_IC_transition_point)+BC_IC_transition_point*D(gamma)(BC_IC_transition_point) = 0

        2.3   Solve numerically an initial value problem for the Lane-Emden equation (result = LaneEmden_ic)
                from  xi = BC_IC_transition_point  to xi = EndPoint with conditions
       gamma(BC_IC_transition_point)    = eval(gamma(xi), LaneEmden_bc(BC_IC_transition_point)
       D(gamma)(BC_IC_transition_point) = eval(diff(gamma(xi), xi), LaneEmden_bc(BC_IC_transition_point)
                

> 

LaneEmden := proc(n, alpha, BC_IC_transition_point, EndPoint)
  local gamma:
  local ode1, ode2, EmdenN, LaneEmden_bc, LaneEmden_ic, f, bc, End_bc, ic:

  ode1   := diff(xi^2*(diff(theta__0(xi), xi)), xi)/xi^2 = -theta__0(xi)^n;
  EmdenN := dsolve({ode1, theta__0(0) = 1, (D(theta__0))(0) = 0}, theta__0(xi), numeric);

  f := proc(z)
    if not type(evalf(z),'numeric') then
      'procname'(z);
    else
      eval(theta__0(xi), EmdenN(z));
    end if;
  end proc;

  printf("Emden = 0 at xi = %a\n", fsolve(f(xi)=0, xi=0..EndPoint));
  ode2         := diff(gamma(xi), xi, xi)-2*gamma(xi)/xi^2+alpha^2*gamma(xi) = -f(xi)^n*xi^2:
  bc           := D(gamma)(0) = 0, gamma(BC_IC_transition_point)+BC_IC_transition_point*D(gamma)(BC_IC_transition_point) = 0;
  LaneEmden_bc := dsolve({ode2, bc}, known=f, numeric, method=bvp[midrich], maxmesh=1024):

  End_bc       := eval([gamma(xi), diff(gamma(xi), xi)], LaneEmden_bc(BC_IC_transition_point));
  ic           := gamma(BC_IC_transition_point) = End_bc[1], D(gamma)(BC_IC_transition_point) = End_bc[2];
  LaneEmden_ic := dsolve({ode2, ic}, known=f, numeric):
  
  DocumentTools:-Tabulate(
    [
      plots:-display(
        plots:-odeplot(
          EmdenN
          , [xi, theta__0(xi)]
          , xi=0..EndPoint
          , color=blue
          , gridlines=true
          , title = cat("Emden, n=", N, ", numerical solution")
        )
      ),
      
      plots:-display(
        plots:-odeplot(
          LaneEmden_bc
          , [xi, gamma(xi)]
          , xi=0..BC_IC_transition_point
          , color=blue
          , gridlines=true
          , title = cat("Lane-Emden, n=", N, ", numerical solution")
        ),
        plots:-odeplot(
          LaneEmden_ic
          , [xi, gamma(xi)]
          , xi=BC_IC_transition_point..EndPoint
          , color=red
          , gridlines=true
          , title = cat("Lane-Emden, n=", N, ", numerical solution")
        )
      ),
      
      plots:-display(
        plots:-odeplot(
          LaneEmden_bc
          , [xi, gamma(xi)+xi*D(gamma)(xi)]
          , xi=0..BC_IC_transition_point
          , color=blue
          , gridlines=true
          , title = typeset('gamma(xi)+xi*D(gamma)(xi)', `, n`=N)
        ),
        plots:-odeplot(
          LaneEmden_ic
          , [xi, gamma(xi)+xi*D(gamma)(xi)]
          , xi=BC_IC_transition_point..EndPoint
          , color=red
          , gridlines=true
          , title = typeset('gamma(xi)+xi*D(gamma)(xi)', `, n`=N)
        )
      )
    ], width=80
  ):
end proc:
> 

LaneEmden(1, 1.4, Pi, 10)

Emden = 0 at xi = 3.141592311

  
> 

LaneEmden(3, 1.4, Pi, 10)

Emden = 0 at xi = 6.896848025

  
> 

LaneEmden(3, 1.4, 6.896848025, 10)

Emden = 0 at xi = 6.896848025

  
> 

 

Download General_Emden_and_Lane-Emden_2.mw

First 7 8 9 10 11 12 13 Last Page 9 of 125