Advanced Search

tomleslie

13821 Reputation

20 Badges

14 years, 290 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

Well...

tomleslie 13821
April 07 2023
0 1

I'd read the whole data set from EXxcel using ExcelTools:-Import() and then use Maple to select/organise dependent and independent variables for fitting.

A "toy" example is shown in the attached, where each fit outputs [rsquared, rsquaredadjusted, tprobability]. from the help page

?Solutions Returned by the Regression Commands

one has

  • rsquared -- The coefficient of determination, often used to indicate how well data fits a model.
  • rsquaredadjusted -- The adjusted coefficient of determination that accounts for the number of variables in the model.
  • tprobability -- The p-value for the hypothesis test for which the t value is the test statistic.

This ought(?) to be the answer to the second question which you have asked.

 

In this "toy" example the values for all of the dependent and independent variables are fictitious, as is the model to which they are being fitted so don't expect fits to be much good

By the way if you really have 30 columns of independent variables and you want to perform a fit to any pair of columns  - this will mean performing 870 fit operations -you may want to deal with a "smaller" example whilst writing/debugging your code

> 

  restart:
  with(LinearAlgebra):
  ExperimentalData := < v, 0.531, 0.341, 0.163, 0.641, 0.713, -0.040
                      | x, 1, 1, 1, 2, 2, 2
                      | y, 1, 2, 3, 1, 2, 3
                      | z, 1, 2, 3, 4, 5, 6
                      >;
  for j in combinat:-permute(op([1,2], ExperimentalData)-1, 2) do
      labs:=`<|>`(Column(ExperimentalData[1,..], [(j+~1)[], 1]));
      data:=`<|>`(Column(ExperimentalData[2..-1,..], [(j+~1)[], 1]));
      Statistics:-LinearFit( labs[1,1]*a + b*labs[1,1]^2/labs[1,2],
                             data,
                             convert(labs[1, 1..2], list),
                             output=[rsquared, rsquaredadjusted, tprobability]
                           );
  od;

Matrix(7, 4, {(1, 1) = v, (1, 2) = x, (1, 3) = y, (1, 4) = z, (2, 1) = .531, (2, 2) = 1, (2, 3) = 1, (2, 4) = 1, (3, 1) = .341, (3, 2) = 1, (3, 3) = 2, (3, 4) = 2, (4, 1) = .163, (4, 2) = 1, (4, 3) = 3, (4, 4) = 3, (5, 1) = .641, (5, 2) = 2, (5, 3) = 1, (5, 4) = 4, (6, 1) = .713, (6, 2) = 2, (6, 3) = 2, (6, 4) = 5, (7, 1) = -0.40e-1, (7, 2) = 2, (7, 3) = 3, (7, 4) = 6})

  

Vector[row](3, {(1) = x, (2) = y, (3) = v})

  

Matrix(6, 3, {(1, 1) = 1, (1, 2) = 1, (1, 3) = .531, (2, 1) = 1, (2, 2) = 2, (2, 3) = .341, (3, 1) = 1, (3, 2) = 3, (3, 3) = .163, (4, 1) = 2, (4, 2) = 1, (4, 3) = .641, (5, 1) = 2, (5, 2) = 2, (5, 3) = .713, (6, 1) = 2, (6, 2) = 3, (6, 3) = -0.40e-1})

  

[.730248825887470, .595373238831206, [.590040726626360, .380105705667960]]

  

Vector[row](3, {(1) = x, (2) = z, (3) = v})

  

Matrix(6, 3, {(1, 1) = 1, (1, 2) = 1, (1, 3) = .531, (2, 1) = 1, (2, 2) = 2, (2, 3) = .341, (3, 1) = 1, (3, 2) = 3, (3, 3) = .163, (4, 1) = 2, (4, 2) = 4, (4, 3) = .641, (5, 1) = 2, (5, 2) = 5, (5, 3) = .713, (6, 1) = 2, (6, 2) = 6, (6, 3) = -0.40e-1})

  

[.821717463177158, .732576194765737, [.668002244228826, .133784018957921]]

  

Vector[row](3, {(1) = y, (2) = x, (3) = v})

  

Matrix(6, 3, {(1, 1) = 1, (1, 2) = 1, (1, 3) = .531, (2, 1) = 2, (2, 2) = 1, (2, 3) = .341, (3, 1) = 3, (3, 2) = 1, (3, 3) = .163, (4, 1) = 1, (4, 2) = 2, (4, 3) = .641, (5, 1) = 2, (5, 2) = 2, (5, 3) = .713, (6, 1) = 3, (6, 2) = 2, (6, 3) = -0.40e-1})

  

[.533423416834850, .300135125252275, [.142994934540483, .281965124452022]]

  

Vector[row](3, {(1) = y, (2) = z, (3) = v})

  

Matrix(6, 3, {(1, 1) = 1, (1, 2) = 1, (1, 3) = .531, (2, 1) = 2, (2, 2) = 2, (2, 3) = .341, (3, 1) = 3, (3, 2) = 3, (3, 3) = .163, (4, 1) = 1, (4, 2) = 4, (4, 3) = .641, (5, 1) = 2, (5, 2) = 5, (5, 3) = .713, (6, 1) = 3, (6, 2) = 6, (6, 3) = -0.40e-1})

  

[.380215523936068, 0.703232859041014e-1, [.405926630635241, .698700881346495]]

  

Vector[row](3, {(1) = z, (2) = x, (3) = v})

  

Matrix(6, 3, {(1, 1) = 1, (1, 2) = 1, (1, 3) = .531, (2, 1) = 2, (2, 2) = 1, (2, 3) = .341, (3, 1) = 3, (3, 2) = 1, (3, 3) = .163, (4, 1) = 4, (4, 2) = 2, (4, 3) = .641, (5, 1) = 5, (5, 2) = 2, (5, 3) = .713, (6, 1) = 6, (6, 2) = 2, (6, 3) = -0.40e-1})

  

[.917247152618675, .875870728928013, [0.583110431704181e-2, 0.968590476544140e-2]]

  

Vector[row](3, {(1) = z, (2) = y, (3) = v})

  

Matrix(%id = 36893488148279964604)

  

[HFloat(0.541823133651818), HFloat(0.312734700477727), [HFloat(0.995893258860949), HFloat(0.4743896868797166)]]

 (1)
> 

 

Download fitdata.mw

See ther attached...

tomleslie 13821
April 07 2023
2 1

where  some typos are fixed.

I don't think this code will *ever* produce a meaningful answer, because there is no constraint on x[4], so it can be increased arbitrarily, with corresponding increases in the values of f() and g()

restart

f := proc (x) options operator, arrow; sum(ln(x[i]), i = 1 .. 4) end proc

proc (x) options operator, arrow; sum(ln(x[i]), i = 1 .. 4) end proc

 (1)

g := proc (x) options operator, arrow; ln(x[1])+ln(x[2])+ln(x[3])+ln(x[4]) end proc

proc (x) options operator, arrow; ln(x[1])+ln(x[2])+ln(x[3])+ln(x[4]) end proc

 (2)

f([1, 1, 1, 1])

0

 (3)

cons := {x[1]+x[2] <= 1, x[1]+x[3] <= 2}

{x[1]+x[2] <= 1, x[1]+x[3] <= 2}

 (4)

with(Optimization); ip := [x[1] = .5, x[2] = .5, x[3] = 1, x[4] = 100]

[x[1] = .5, x[2] = .5, x[3] = 1, x[4] = 100]

 (5)

Maximize(f(x), cons, assume = nonnegative, initialpoint = ip)

[16.4237121534987303, [x[1] = HFloat(0.5028668384263846), x[2] = HFloat(0.49713316157361515), x[3] = HFloat(1.4971331615736156), x[4] = HFloat(3.62694858829851e7)]]

 (6)

Maximize(f(x), {x[1]+x[2] <= 1, x[1]+x[3] <= 2}, assume = nonnegative, initialpoint = ip)

[16.4237121534987303, [x[1] = HFloat(0.5028668384263846), x[2] = HFloat(0.49713316157361515), x[3] = HFloat(1.4971331615736156), x[4] = HFloat(3.62694858829851e7)]]

 (7)

Maximize(g(x), cons, assume = nonnegative, initialpoint = ip)

[16.4237121534987303, [x[1] = HFloat(0.5028668384263846), x[2] = HFloat(0.49713316157361515), x[3] = HFloat(1.4971331615736156), x[4] = HFloat(3.62694858829851e7)]]

 (8)

Maximize(g(x), {x[1]+x[2] <= 1, x[1]+x[3] <= 2}, assume = nonnegative, initialpoint = ip)

[16.4237121534987303, [x[1] = HFloat(0.5028668384263846), x[2] = HFloat(0.49713316157361515), x[3] = HFloat(1.4971331615736156), x[4] = HFloat(3.62694858829851e7)]]

 (9)

 

Maximize(sum(ln(x[i]), i = 1 .. 4), {x[1]+x[2] <= 1, x[1]+x[3] <= 2}, assume = nonnegative, initialpoint = ip)

[16.4237121534987303, [x[1] = HFloat(0.5028668384263846), x[2] = HFloat(0.49713316157361515), x[3] = HFloat(1.4971331615736156), x[4] = HFloat(3.62694858829851e7)]]

 (10)

``

NULL

Download opt.mw

Problem with plot...

tomleslie 13821
April 06 2023
1 2

ranges and scaling?? See the attached.

BTW the gridlines do not appear in an actual Maple worksheet, so the plot *looks* much nicer

> 

restart;
TWeq:=2.96736996560705*10^(-12)*p^2+1.31319840299485*10^(-13)*t^2-8.89549693662593*10^(-7)*p+8.53128393394231*10^(-7)*t-3.65558815509970*10^(-30)*p*t-1 = 0;
plots:-implicitplot(TWeq, p = -8*10^6..2*10^6, t=-8*10^6..2*10^6, scaling=constrained)

0.2967369966e-11*p^2+0.1313198403e-12*t^2-0.8895496937e-6*p+0.8531283934e-6*t-0.3655588155e-29*p*t-1 = 0

  

 
> 

 

> 

 

Download ell.mw

Do it numerically...

tomleslie 13821
April 05 2023
0 1

as in the attached

> 

restart; K1 := 10^(-2.12); K2 := 10^(-7.21); K3 := 10^(-12.32); Kw := 10^(-14); NaCl := .5004/(58.44); NaH2PO4 := .1092/(119.98); Na2HPO4 := 1.214/(141.96); Why := Kw = H*OH; Do := K1 = H*H2PO4/H3PO4; It := K2 = H*HPO4/H2PO4; Today := K3 = H*PO4/HPO4; When := Cl = NaCl; Its := Na = NaH2PO4+2*Na2HPO4+NaCl; Due := Na+H = Cl+3*PO4+2*HPO4+H2PO4+OH; Tomorrow := Na2HPO4+NaH2PO4 = H3PO4+H2PO4+HPO4+PO4; z := fsolve([Why, Do, It, Today, When, Its, Due, Tomorrow]); eval(H2PO4, z)

{Cl = 0.8562628337e-2, H = -0.2274751800e-1, H2PO4 = -0.4733982564e-2, H3PO4 = 0.1419582613e-1, HPO4 = 0.1283194935e-7, Na = 0.2657618944e-1, OH = -0.4396084003e-12, PO4 = -0.2699968014e-18}

  

-0.4733982564e-2

 (1)
> 

 

Download ChemNum.mw

Something *odd* about the implcitplot3d...

tomleslie 13821
April 04 2023
0 0

and I'm not sure what it is!

Luckily, you are only using it to produce the vertical plane x=y, which is trivial to do uing the plottools:polygon() command. Everything now works in the attached. I did change some colors and radii on the tubeplots, to make clearer (to me) what was going on. You may want to change these back, or reset to other values. Same issue with "scaling", Some of the plots are on rather different scales, which makes combining them look a bit odd - you may want to consider some different plot ranges, or even scaling=constrained.

The actual plots look a lot better in the worksheet than they render on this site - honest!
 

Relud surface(2)

 

 

> 

restart;
with(plots):
with(plottools):
> 

G1 := tubeplot([4*t, -1, t], t = -10 .. 10, radius = 2, color = blue):
G2 := tubeplot([1, t, t^2], t = -15 .. 15, radius = 2, color = cyan):
G3 := plottools:-polygon([[-5, -5, -4], [-5, -5, 4], [5, 5, 4], [5, 5, -4]], color = red):
bk := display([G1, G2, G3]);

 
> 

Line:=(q,t)->[4*q+(1-4*q)*t,-1+(1-4*q)*t,q+(16*q^2-q)*t]:
F:= Q-> display(  line(Line(Q,-1.5),Line(Q,+2.5),color=red,thickness=5),
                  plot3d([4*q+y+1,y,q+(y+1)*(16*q^2-q)/(1-4*q)], q=-4..Q, y=(-2.5+6*q)..(1.5-10*q),
                         style=patchcontour, lightmodel=light4, axes=framed)
                        ):
> 

animate(F, [q], q=-4..4, frames=50, lightmodel=light4,background=bk);

 

NULL

Download oddPlot.mw

Maybe...

tomleslie 13821
April 03 2023
0 5

something like this?

(which won't display "inline", becauase this site is broken AGAIN

fit.mw

So many ways to do this...

tomleslie 13821
April 03 2023
0 0

The attached shows three of them

  1. using the Student[MultivariateCalculus]() package;
  2. using the LinearAlgebra() package
  3. using the geom3d() package

(And this site is failing to produce an inline display of wroksheets -AGAIN!)

ptline.mw

 

You are using...

tomleslie 13821
March 22 2023
0 1

the procedures f1() and f2() in generating equations in e1.. e6, however you are calling these with "square" brackets, eg f1[n]. These should be "round" brackets. In Maple "square" brackets are only used to access entries in an indexable data container

I have fixed this problem in the attached, where I have also restricted the extent of the main "for" loop, just to avoid generating too much output for this page,. Also, just for display purposes, I slightly modified the format strings in the main printf() statement, just for clarity.

For some reason this worksheet won't display inline here, so will you have to download to check

fsol.mw

Maybe...

tomleslie 13821
March 22 2023
1 0

using PDEtools:-dchange(), as shown in the attached?

> 

restart
> 

"T(x,eta):=Te+Te*(1-g(x,eta))*phi(x); PDE:=diff(T(x,eta),eta$2)+ f(x,eta)*diff(T(x,eta),eta)=x*(diff(f(x,eta),eta)*diff(T(x,eta),x)-diff(T(x,eta),eta)*diff(f(x,eta),x)) "

proc (x, eta) options operator, arrow, function_assign; Te+Te*(1-g(x, eta))*phi(x) end proc

  

-Te*(diff(diff(g(x, eta), eta), eta))*phi(x)-f(x, eta)*Te*(diff(g(x, eta), eta))*phi(x) = x*((diff(f(x, eta), eta))*(-Te*(diff(g(x, eta), x))*phi(x)+Te*(1-g(x, eta))*(diff(phi(x), x)))+Te*(diff(g(x, eta), eta))*phi(x)*(diff(f(x, eta), x)))

 (1)
> 

> 

NULL
> 

NULL
> 

NULL
> 

PDEtools:-dchange(eta = (ue/(c*x))^(1/2)*y, PDE, [y], params = [ue, c])

-Te*(diff(diff(g(y, x), y), y))*c*x*phi(x)/ue-f(y, x)*Te*(diff(g(y, x), y))*phi(x)/(ue/(c*x))^(1/2) = x*((diff(f(y, x), y))*(-Te*(diff(g(y, x), x))*phi(x)+Te*(1-g(y, x))*(diff(phi(x), x)))/(ue/(c*x))^(1/2)+Te*(diff(g(y, x), y))*phi(x)*(diff(f(y, x), x))/(ue/(c*x))^(1/2))

 (2)

NULL

Download chngeVar.mw

Or...

tomleslie 13821
March 21 2023
0 0

you could use complexplot(), as in the attached

> 

restart
> 

with(plottools)
> 

with(plots)
> 

with(CurveFitting)
> 

Digits := 10
> 

with(GaussInt)
> 

w := GInearest(0+I)

I

 (1)
> 

NULL
> 

"f(t):=7.0*(e)^((-(t-13180)^(2))/(2000000))+4.7*(e)^((-(t-16000)^(2))/(3200000)):"
> 

p1 := plot(f(t), t = 0 .. 20000, color = green); plots[display]({p1})

 
> 

NULL
> 

D1 := 15
> 

epsilon := 200000
> 

L := 6500
> 

v := .7

.7

 (2)
> 

t := 10000

10000

 (3)
> 

i = sqrt(-1)

i = I

 (4)
> 

"k(n) := evalf((2 *Pi*n)/(L))"

proc (n) options operator, arrow, function_assign; evalf(2*Pi*n/L) end proc

 (5)
> 

f(n) = (int(f(t)*exp(-w*k(n)*x), x = 0 .. L))/L
> 

``
> 

"C(x, t) := (&sum;) exp(-v* t*k(n)- D1 *t*(k(n)^())^(2)- epsilon *t*(k(n))^(4)) *f(n)* exp(w*k(n)* x)"

proc (x, t) options operator, arrow, function_assign; sum(exp(-v*t*k(n)-D1*t*k(n)^2-varepsilon*t*k(n)^4)*f(n)*exp(w*k(n)*x), n = 1 .. 10) end proc

 (6)
> 

uu10000 := [seq(evalf(C(L-j, t)), j = 0 .. 6500, 100)]; complexplot(uu10000, x = min(`~`[Re](uu10000)) .. max(`~`[Re](uu10000)), scaling = constrained, style = point)

 
> 

``

Download cplxplt.mw

The OP's primary objective...

tomleslie 13821
March 19 2023
0 0

is to solve the ODE given as f'''+ff''-f'^2-Mf'=0, with f(0)=0, f'(0)=1, and f'(5)=0 with M=0.5

So let's just solve it, as shown in the attached.

> 

  restart:
  with(plots):
> 

  M:=0.5;
  ode:= diff(f(x), x$3)+f(x)*diff(f(x), x$2)-diff(f(x),x)^2-M*diff(f(x),x)=0;
  conds:= f(0)=0, D(f)(0)=1, D(f)(5 )-0;
  ans:=dsolve( [ode, conds], numeric);
  odeplot(ans, [x, f(x)], x=0..5);
  

.5

  

diff(diff(diff(f(x), x), x), x)+f(x)*(diff(diff(f(x), x), x))-(diff(f(x), x))^2-.5*(diff(f(x), x)) = 0

  

f(0) = 0, (D(f))(0) = 1, (D(f))(5)

  

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) = .0, (2) = .3648083980208248, (3) = .739453123146836, (4) = 1.1294231031240023, (5) = 1.5404396951722952, (6) = 1.9674194416734547, (7) = 2.4067451277606597, (8) = 2.850771504517652, (9) = 3.2971678108063687, (10) = 3.7442923732839364, (11) = 4.191663594165653, (12) = 4.620140941040035, (13) = 5.0}, datatype = float[8], order = C_order); Y := Matrix(13, 3, {(1, 1) = .0, (1, 2) = 1.0, (1, 3) = -1.2249064035636923, (2, 1) = .294194453432655, (2, 2) = .6396128765313516, (2, 3) = -.7836115514965457, (3, 1) = .48635488180500724, (3, 2) = .40415085049570865, (3, 3) = -.49534834686074514, (4, 1) = .6116367675260377, (4, 2) = .250541601628673, (4, 3) = -.3073631245893844, (5, 1) = .6925060924519302, (5, 2) = .15125045647457197, (5, 3) = -.18593456512252354, (6, 1) = .7427378377645504, (6, 2) = 0.893905555320538e-1, (6, 3) = -.110380533553483, (7, 1) = .7730294979143654, (7, 2) = 0.51843487043024994e-1, (7, 3) = -0.6463752743281678e-1, (8, 1) = .7906832855559992, (8, 2) = 0.29654107059807183e-1, (8, 3) = -0.3774106284132804e-1, (9, 1) = .8007523087394047, (9, 2) = 0.1661644456196504e-1, (9, 3) = -0.22098378262475696e-1, (10, 1) = .8063158313666053, (10, 2) = 0.8938405451087109e-2, (10, 3) = -0.1307531449846363e-1, (11, 1) = .8092030889605896, (11, 2) = 0.4352616339821983e-2, (11, 3) = -0.7907830531921609e-2, (12, 1) = .810440946160643, (12, 2) = 0.1626507257535621e-2, (12, 3) = -0.5077258781403503e-2, (13, 1) = .8107323262402987, (13, 2) = .0, (13, 3) = -0.3612766608982084e-2}, datatype = float[8], order = C_order); YP := Matrix(13, 3, {(1, 1) = 1.0, (1, 2) = -1.2249064035636923, (1, 3) = 1.5, (2, 1) = .6396128765313516, (2, 2) = -.7836115514965457, (2, 3) = .9594452421864268, (3, 1) = .40415085049570865, (3, 2) = -.49534834686074514, (3, 3) = .6063284218940224, (4, 1) = .250541601628673, (4, 2) = -.3073631245893844, (4, 3) = .37603648294155106, (5, 1) = .15125045647457197, (5, 2) = -.18593456512252354, (5, 3) = .2272627479658001, (6, 1) = 0.893905555320538e-1, (6, 2) = -.110380533553483, (6, 3) = .13466974800716747, (7, 1) = 0.51843487043024994e-1, (7, 2) = -0.6463752743281678e-1, (7, 3) = 0.7857620604810917e-1, (8, 1) = 0.29654107059807183e-1, (8, 2) = -0.3774106284132804e-1, (8, 3) = 0.4554764716317479e-1, (9, 1) = 0.1661644456196504e-1, (9, 2) = -0.22098378262475696e-1, (9, 3) = 0.26279655923937467e-1, (10, 1) = 0.8938405451087109e-2, (10, 2) = -0.1307531449846363e-1, (10, 3) = 0.1509193089776011e-1, (11, 1) = 0.4352616339821983e-2, (11, 2) = -0.7907830531921609e-2, (11, 3) = 0.8594294332320503e-2, (12, 1) = 0.1626507257535621e-2, (12, 2) = -0.5077258781403503e-2, (12, 3) = 0.4930717565329715e-2, (13, 1) = .0, (13, 2) = -0.3612766608982084e-2, (13, 3) = 0.29289866770633205e-2}, 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) = .0, (2) = .3648083980208248, (3) = .739453123146836, (4) = 1.1294231031240023, (5) = 1.5404396951722952, (6) = 1.9674194416734547, (7) = 2.4067451277606597, (8) = 2.850771504517652, (9) = 3.2971678108063687, (10) = 3.7442923732839364, (11) = 4.191663594165653, (12) = 4.620140941040035, (13) = 5.0}, datatype = float[8], order = C_order); Y := Matrix(13, 3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = 0.1273975608644069e-9, (2, 1) = -0.9699958654485496e-7, (2, 2) = 0.11898009711795105e-6, (2, 3) = -0.14531470840848373e-6, (3, 1) = -0.10320306463506545e-6, (3, 2) = 0.12691811465663893e-6, (3, 3) = -0.154678061248875e-6, (4, 1) = -0.7871067918956964e-7, (4, 2) = 0.9744291713856059e-7, (4, 3) = -0.11818172410834043e-6, (5, 1) = -0.4865437588049475e-7, (5, 2) = 0.613589205139654e-7, (5, 3) = -0.7358044074399424e-7, (6, 1) = -0.21030326071357806e-7, (6, 2) = 0.28453569331208245e-7, (6, 3) = -0.3288511507213022e-7, (7, 1) = -0.12272437930054753e-8, (7, 2) = 0.5323548042522246e-8, (7, 3) = -0.4173055146179112e-8, (8, 1) = 0.9931557233651253e-8, (8, 2) = -0.7035041293820884e-8, (8, 3) = 0.11351363226824974e-7, (9, 1) = 0.14250962440952772e-7, (9, 2) = -0.10817614549677532e-7, (9, 3) = 0.1640579916342742e-7, (10, 1) = 0.14719058046583237e-7, (10, 2) = -0.965367179587185e-8, (10, 3) = 0.15469412664652983e-7, (11, 1) = 0.13753746410766736e-7, (11, 2) = -0.6455085789179595e-8, (11, 3) = 0.12146772829413411e-7, (12, 1) = 0.1268333027160663e-7, (12, 2) = -0.2912324487441674e-8, (12, 3) = 0.8508819311976027e-8, (13, 1) = 0.12136283725856026e-7, (13, 2) = .0, (13, 3) = 0.5688098437791177e-8}, 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(1.54678061248875e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [3, 13, [f(x), diff(f(x), x), diff(diff(f(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(3, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(13, 3, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(3, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(13, 3, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[f(x), diff(f(x), x), diff(diff(f(x), x), x)]'[i] = yout[i], i = 1 .. 3)] 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(1.54678061248875e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [3, 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(3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(13, 3, 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(3, {(1) = 0., (2) = 0., (3) = 0.}); `dsolve/numeric/hermite`(13, 3, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 3)] end proc, (2) = Array(0..0, {}), (3) = [x, f(x), diff(f(x), x), diff(diff(f(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)]'[i] = res[i+1], i = 1 .. 3)] catch: error end try end proc

  

 
> 

 

Download solODE.mw

8-th order polynomial...

tomleslie 13821
March 17 2023
0 0

which if you differentiate it, will give a 7-th order polynomial. Setting the latter equal to zero will rsult in seven roots (ie seven stationary points), although as mmcdara has pointed out three of these occur at x=-1. \Whether you consider this a single stationary point (with a multiplicity of 3), or three coincident stationary points is a quibble about language

The other four roots are complex

Second vertical axis is not logarithmic...

tomleslie 13821
March 17 2023
1 1

Check the tickmarks!

The problem appears to be that you have two sets of horizontal gridlines - one for each of the vertical axes in the dualaxis plot. Since the left and right tickmarks do not align, these two sets of gridlines are interleaved, giving an "illusion" of a weird vertical scale.

In a dualaxis plot, I can't think of any way to make the left and right tickmarks "align", since they are likely to be on very different scales.

Even if one could come up with a method of having horizontal gridlines only for the left axis (and I don't know if this can be done!), it means that these gridlines would not coincide with the tickmarks on the right axis - which would make the plot look a bit weird

anything...

tomleslie 13821
March 16 2023
0 0

wrong with the fsolve() command using the 'complex' option?

You couldread this response for some sug...

tomleslie 13821
March 14 2023
0 0

https://www.mapleprimes.com/questions/235892-How-Do-I-Solve-A-Differential-Equation-In-Maple

3 4 5 6 7 8 9 Last Page 5 of 207