sand15

Ask a Question

Create a Post

11 Badges

Member for: 8 years, 311 days

Contact sand15

MaplePrimes Activity

These are answers submitted by sand15

Have you any idea of the size of Bd?...

Answered: sand15 615

September 28 2023

0 6

restart
B := Matrix(9$2, symbol=m):
Bd:= CodeTools:-Usage( LinearAlgebra:-Determinant(B) ):
memory used=0.82GiB, alloc change=0.60GiB, cpu time=17.64s, real time=11.98s, gc time=9.07s

# approximately the number of symbols used to represent Bd

L := length(Bd)
                            40642561
# printing Bd on a 400 charaters per line display will require
# this number of lines

iquo(L, 400) + `if`(irem(L, 400)=0, 0, 1)
                             101607

# And thus a pocket book (50 characters per line and 60 lines
# per page) of about this number of pages

round(evalf(L / 50 / 60))
                             13548

So how can you even read this expression?

I would say no...

Answered: sand15 615

September 21 2023

1 1

On might think using the package ImageTools, providing the image has an appropriate format and is of good quality.
But it will be quite difficult to extract the contours of the points, take their centers, extract the correct scales (did I miss comething?)

My advice is to use a digitizer and do the things manually (some can be use online, for instance https://plotdigitizer.com/app)

Computation of the determinant of B and ...

Answered: sand15 615

September 21 2023

1 3

... is to compute the determinant X of some symbolic matrix A of elements x_{i, j}, and next to replace in the expression of X each x_{i, j} by the corresponding expression in B (x[i, j]=B[i, j])
The corresponding lines are in brown courier font.

Download determinant.mw

EDITED:
Following the (judicious) remark @acer made, here is a simple say to get the determinant of B1 (not the one of Omeg):

detH  := Determinant(H):
detB1 := Determinant(B1):
detB  := detH * detB1:

(1)

assume(x::real);

(2)

>	t0 := time(): detH := Determinant(H): detB1 := Determinant(B1): detB := detH * detB1: time()-t0;

(3)

>	simplify(detH, size); simplify(detB1, size)

(4)

Download determinant_edited.mw

As matrix Omeg := B + idn8 = H.B1 + idn8 (BTW idn8 := IdentityMatrix(8)) and because H:=Matrix(8, 8, symbol=h), (thus H is formally equivalent to the matrix X used in the "trick"), I think the "trick" can be applied to compute Determinant(Omeg).
Of course you must be aware of @acer's remark if you have to evaluate Determinant(Omeg) for some instance of H.

What loop are you talking about? A_...

Answered: sand15 615

September 11 2023

0 8

UPDATED

~~What loop are you talking about?~~ (shouldn't have been here)

Percentage of variation of peak heights

A_ref := eval(A, lambda = 1.3015):
x_ref := fsolve(diff(A_ref, x), x=-2..0);
y_ref := evalf(eval(A_ref, x=x_ref));

                         -0.7765307466
                          1.920953647

data_1 := [beta=0.1, Q=1.3015, lambda=0.9986];
B_1    := eval(B, data_1):
x_1    := fsolve(diff(B_1, x), x=-2..0);
y_1    := evalf(eval(B_1, x=x_1));

PercentageOfVariation := (y_1-y_ref)/y_ref*100

           [beta = 0.1, Q = 1.3015, lambda = 0.9986]
                         -0.8959052730
                          1.752467429
                          -8.770967392

data_2 := [beta=0.3, Q=1.3015, lambda=0.9986];
B_2    := eval(B, data_2):
dB_2   := diff(B_2, x):
x_2    := fsolve(diff(B_2, x), x=-2..0);
y_2    := evalf(eval(B_2, x=x_2));

PercentageOfVariation := (y_2-y_ref)/y_ref*100

           [beta = 0.3, Q = 1.3015, lambda = 0.9986]
                          -1.088722291
                          1.547856941
                          -19.42247313

The main problem with Word is that the c...

Answered: sand15 615

August 26 2023

0 0

The main problem with Word is that the character sizes will automatically decrease as the "depth" of the continued fraction increases (I don't know if it's possible to overcome this?) :

Why not use the following notation instead (see https://en.wikipedia.org/wiki/Continued_fraction, section Motivation and notation) and just copy-paste the result into your MS-Word document (next convert it as an equation if you prefer)?

phi := convert((1+sqrt(5))*(1/2), confrac)
               [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
# or
phi := [phi[], `#mo("&#x2026;")`]
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...]

See here https://www.wikihow.life/Start-Working-with-Continued-Fractions for some advices, in particular topic (4) "to avoid the cumbersome staircase notation"
(see also Pringsheim notation in the same reference).

PS: note that the conventional notation should be [1 ; 1, 1, ...] instead of [1 , 1, 1, ...].

You already found the solution...

Answered: sand15 615

August 18 2023

0 1

in your previous question
https://www.mapleprimes.com/questions/236897-Differentiate-Three-Regions-With-Implicitplot#comment296558

FirmModelPP := proc (alpha, delta) local p0, xi0, q0, Firmpf0, G0, Recpf0, Unsold0, Environ0; xi0 := 1; p0 := min(s+sqrt((v-s)*(c-s)), delta*v+sExp); q0 := u*(v-p0)/(v-s); f(N) := 1/u; F(N) := N/u; G0 := int(F(N), N = 0 .. q0); Firmpf0 := (p0-c)*q0-(p0-s)*G0; Recpf0 := (sExp-cr)*xi0*q0; Environ0 := q0+G0; Unsold0 := G0; return p0, q0, Firmpf0, Recpf0, Environ0, Unsold0 end proc:

pltPP3 := plot(['FirmModelPP(alpha, .2)[3]'], 'alpha' = 0.1e-1 .. .99, 'linestyle' = ['solid'], 'legend' = ["Current PP strategy with delta=0.2"], 'labels' = [alpha, "Firm profit"], 'labeldirections' = ["horizontal", "vertical"], 'color' = ['blue'], 'axes' = 'boxed'):

pltFC3 := plot(['FirmModelFC(alpha, .2, .2)[3]'], 'alpha' = 0.1e-1 .. .99, 'linestyle' = ['solid'], 'legend' = ["Current FC strategy with delta=0.2"], 'labels' = [alpha, "Firm profit"], 'labeldirections' = ["horizontal", "vertical"], 'color' = ['red'], 'axes' = 'boxed'):

$" pltHmax3:= plot([seq('FirmModelHmax(alpha,0.2,delta)[3]' ,delta=0.1..0.5,0.2)] , alpha=0..1,linestyle=[dot,dashdot,dash] ,legend=[seq('delta'=delta,delta=0.1..0.5,0.2)] , legendstyle=[location=left] ,labels=["\""`alpha"`,"SiS firm profit"] ,labeldirections = ["horizontal", "vertical"] ,legendstyle=[location=bottom]): "$

Error, Got internal error in Typesetting:-Parse:-Postprocess : "internal error: invalid object ""

Error, (in plots:-display) expecting plot structures but received: [pltHmax3]

WhyNot := proc (alpha, delta) if not [alpha, delta]::(list(numeric)) then return ('procname')(args) end if; FirmModelHmax(alpha, .2, delta)[3] end proc:

pltHmax3 := plot(
  [seq(WhyNot(alpha, delta),delta=0.1..0.5,0.2)]
  , alpha=0..1
  , linestyle=[dot,dashdot,dash]
  , legend=[seq('delta'=delta,delta=0.1..0.5,0.2)]
  , legendstyle=[location=left]
  , labels=["alpha","SiS firm profit"]
  , labeldirections =["horizontal", "vertical"]
  , legendstyle=[location=bottom]
):

display(pltHmax3);

>	display([pltPP3, pltFC3, pltHmax3])

Download WhyNot.mw

After diffr2 and diffr3 corrected...

Answered: sand15 615

August 17 2023

2 3

And after using some trick, here is the result you wanted (whose rendering can likely be improved)
Compare_three_regions_sand15.mw

A way...

Answered: sand15 615

August 15 2023

1 0

restart; with(LinearAlgebra); A := `<,>`(`<|>`(c, 0, 0), `<|>`(0, -1, a), `<|>`(0, b, -1))

(1)

(2)

Warning, solutions may have been lost

(3)

for u in [s] do `assuming`([`not`(is(z[2] >= 0))], [op(u)]); `assuming`([is(z[3] < 0)], [op(u)]); print() end do

(4)

for u in [s] do `assuming`([coulditbe(z[2] < 0)], [op(u)]); `assuming`([coulditbe(z[3] < 0)], [op(u)]); print() end do

(5)

Download A_way.mw

@Rana47 The solution I get (note th...

Answered: sand15 615

July 18 2023

1 6

@Rana47

The solution I get (note that either I did something wrong, either the reference you provide contains a typo, see eq. (10) in the attached file... but it looks like you made the same observation on yout own).

>	restart: # In the sequel [REF] refers to the paper in your question

>	PDEtools[declare](f(x), prime = x); PDEtools[declare](v(x), prime = x);

(1)

>	ansatz := N -> g(x) = add((p^(i))*f[i](x), i = 0 .. N);

(2)

>	HPMEq0 := n -> (1-p)diff(g(x), x$2)+p(diff(g(x), x$2)+nb(diff(g(x), x))^(n-1)*(diff(g(x), x$2))-k)

proc (n) options operator, arrow; (1-p)*(diff(g(x), `$`(x, 2)))+p*(diff(g(x), `$`(x, 2))+n*b*(diff(g(x), x))^(n-1)*(diff(g(x), `$`(x, 2)))-k) end proc

(3)

>	# For future comparison. # I assume that "b could be small" BandK := {b=0.2, k=1}: solnum := dsolve(eval({HPMEq0(5), g(0)=1, D(g)(1)=0}, [b=1, k=1, p=1]), numeric);

(4)

>	ansatz_order := 3; exponent_degree := 5; # The initial guess corresponding to relation [REF] eq (29). v[0] := k/2(x^2-2x)+1; Lv[0] := diff(v[0], x$2);

(5)

>	# Plug the ansatz into HPMEq0: HPMEq1 := collect(eval(HPMEq0(exponent_degree), ansatz(ansatz_order)), p):

>	# Derive the successive edos according to [REF] eq (23), (25), (27) and # generalize. equ[0] := coeff(HPMEq1, p, 0) - diff(v[0], x$2) = 0; # [REF] eq 23 for i from 1 to ansatz_order do equ[i] := coeff(HPMEq1, p, i) + diff(v[0], x$(i+1)) = 0 end do;

(6)

>	# BCs given by eq (21) in the paper you present. cond[0] := f[0](0) = 1, (D(f[0]))(1) = 0; # BCs on the successive corrections for j from 1 to ansatz_order do cond[j] := f[j](0) = 0, (D(f[j]))(1) = 0 end do;

(7)

>	answers := {dsolve({cond[0], equ[0]})}; for j from 1 to ansatz_order do ans := dsolve(eval({cond[j], equ[j]}, answers)): answers := answers union {ans} end do: #print~(answers):

${f[0](x) = (1/2)*k*x^2-k*x+1}$

(8)

>	# Check if f[1] is equal to the 3rd term in [REF] eq (33) # (OK) collect(eval(f[1](x), answers), [b, k]);

(9)

>	# Check if f[2] is equal to the 4th term in [REF] eq (32a) # (Not OK: coefficient "5" missing in the numerator, # could it be an error in [REF]?) collect(eval(f[2](x), answers), [b, k]);

(10)

>	# Check [REF] eq (32a) collect(eval(f[0](x), answers), k) + collect(eval(f[1](x), answers), [b, k]) + collect(eval(f[2](x), answers), [b, k]): eval(%, [5=n, 1/6=1/(n+1), 9=2*n-1])

(11)

>	# f[3] is not given in [REF], here is its expression collect(eval(f[3](x), answers), [b, k])

(12)

>	# Order 3 expansion of the solution collect(eval(f[0](x), answers), k) + collect(eval(f[1](x), answers), [b, k]) + collect(eval(f[2](x), answers), [b, k]) + collect(eval(f[3](x), answers), [b, k]);

(13)

ANSWERS := eval(answers, BandK):

plots:-display(
  plots:-odeplot(solnum, [x, g(x)], x=0..1)
  , plot(eval(eval(rhs(ansatz(0)), p=1), ANSWERS), x=0..1, color=blue, legend="0th order")
  , plot(eval(eval(rhs(ansatz(1)), p=1), ANSWERS), x=0..1, color=green, legend="1st order")
  , plot(eval(eval(rhs(ansatz(2)), p=1), ANSWERS), x=0..1, color=cyan, legend="2nd order")
  , plot(eval(eval(rhs(ansatz(3)), p=1), ANSWERS), x=0..1, color=gold, legend="3rd order")
)

>	plots:-display( plot(eval(f[1](x), ANSWERS), x=0..1, color=green, legend="1th order") , plot(eval(f[2](x), ANSWERS), x=0..1, color=cyan, legend="2th order") , plot(eval(f[3](x), ANSWERS), x=0..1, color=gold, legend="3th order") )

Download hpm_sand15_to_rana47.mw

Still need help?...

Answered: sand15 615

July 17 2023

0 0

If my previous reply didn't enlighten you, here is a way to proceed

restart:

>	N := 2: ansatz := f(x) = add((p^(i))*(f\|\|i)(x), i = 0 .. N) ;

(1)

>	HPMEq0 := (1-p)(diff(f(x), `$`(x, 3)))+p(diff(f(x), `$`(x, 3))+2alphaRE(diff(f(x), x))f(x)+(4-H)alpha^2(diff(f(x), x)))

(1-p)*(diff(diff(diff(f(x), x), x), x))+p*(diff(diff(diff(f(x), x), x), x)+2*alpha*RE*(diff(f(x), x))*f(x)+(4-H)*alpha^2*(diff(f(x), x)))

(2)

>	HPMEq1 := eval(HPMEq0, ansatz)

(1-p)*(diff(diff(diff(f0(x), x), x), x)+p*(diff(diff(diff(f1(x), x), x), x))+p^2*(diff(diff(diff(f2(x), x), x), x)))+p*(diff(diff(diff(f0(x), x), x), x)+p*(diff(diff(diff(f1(x), x), x), x))+p^2*(diff(diff(diff(f2(x), x), x), x))+2*alpha*RE*(diff(f0(x), x)+p*(diff(f1(x), x))+p^2*(diff(f2(x), x)))*(f0(x)+p*f1(x)+p^2*f2(x))+(4-H)*alpha^2*(diff(f0(x), x)+p*(diff(f1(x), x))+p^2*(diff(f2(x), x))))

(3)

>	for i from 0 to N do equ[1][i] := coeff(HPMEq1, p, i) = 0 end do;

diff(diff(diff(f0(x), x), x), x) = 0

diff(diff(diff(f1(x), x), x), x)+2*alpha*RE*(diff(f0(x), x))*f0(x)+(4-H)*alpha^2*(diff(f0(x), x)) = 0

diff(diff(diff(f2(x), x), x), x)+2*alpha*RE*(diff(f0(x), x))*f1(x)+2*alpha*RE*(diff(f1(x), x))*f0(x)+(4-H)*alpha^2*(diff(f1(x), x)) = 0

(4)

>	cond[1][0] := f0(0) = 1, (D(f0))(0) = 0, f0(1) = 0; for j from 1 to N do cond[1][j] := (f \|\| j)(0) = 0, (D(f \|\| j))(0) = 0, (f \|\| j)(1) = 0 end do;

(5)

>	answers := {dsolve({cond[1][0], equ[1][0]})}

${f0(x) = -x^2+1}$

(6)

>	H := 5; alpha := (1/180)(53.142); RE := 10

(7)

>	for j from 1 to N do ans := dsolve(eval({cond[1][j], equ[1][j]}, answers)): answers := answers union {ans} end do: print~(evalf(answers)):

(8)

>	plots:-display( plot(eval(f0(x), answers), x=0..1, color=red, legend="0th order"), plot(eval(f0(x)+f1(x), answers), x=0..1, color=green, legend="1st order"), plot(eval(f0(x)+f1(x)+f2(x), answers), x=0..1, color=blue, legend="2nd order"), title="Successive approximations" )

>	plots:-display( plot(eval(f1(x), answers), x=0..1, color=green, legend="1st order"), plot(eval(f2(x), answers), x=0..1, color=blue, legend="2nd order"), title="Successive corrections" )

Download hpm_error_sand15.mw

After having corrected the (likely) typo...

Answered: sand15 615

July 15 2023

1 1

After having corrected the (likely) typo (equal instead of equa) and changed a few little things:

Download HPM_1.mw

You willsee in the attached file that f0=1=...=0(solving directly the initial ode without even using the ansatz

f = x -> add(p^i*(f||i)(x), i = 0 .. L)

already gives f=0 given the bcs you write.

When on uses this kind of ansatz (small parameter expansion?), one generally set non identivally null bcs for the "trend" (f0) a,d null bcs for the perturbations of successive orders (to guarantee that successive refinments preserve the bcs on f0).

With bcs arbitrarily set to

f0(0)=1, D(f0)(0)=1, f0(1)=0

the attached file shows how to compute (in a step-by-step way) perturbations of increasing orders an thus a sequence of augmented solutions.
Once approved this method could be encapsulated into a specific procedure.

>	with(plots):

(1)

(2)

>	HO := (1-p)(diff(f(x), `$`(x, 3))+p(diff(f(x), `$`(x, 3)))+2alphaRE(diff(f(x), x))f(x)+(4-H)alpha^2(diff(f(x), x)))

(1-p)*(diff(diff(diff(f(x), x), x), x)+p*(diff(diff(diff(f(x), x), x), x))+1.745555556*(diff(f(x), x))*f(x)-0.7617410494e-2*(diff(f(x), x)))

(3)

>	ansatz := f = (x -> add(p^i*(f \|\| i)(x), i = 0 .. L));

(4)

>	HO1 := eval(HO, ansatz):

Zero order solution f0

>	HO1 := collect(expand(HO1), p): coeff(HO1, p, 0); sol0 := dsolve({coeff(HO1, p, 0), f0(0)=1, D(f0)(0)=1, f0(1)=0}, numeric);

-0.7617410494e-2*(diff(f0(x), x))+diff(diff(diff(f0(x), x), x), x)+1.745555556*(diff(f0(x), x))*f0(x)

$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) = .0, (2) = .11078278220502231, (3) = .22364174273643112, (4) = .336400857639379, (5) = .44797504433316043, (6) = .5582094408957031, (7) = .6677790172685235, (8) = .7779200458343468, (9) = .8895341149733956, (10) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(10, 3, {(1, 1) = 1.0, (1, 2) = 1.0, (1, 3) = -3.983477327766739, (2, 1) = 1.0859803004343302, (2, 2) = .549346065085729, (2, 3) = -4.1393578797975605, (3, 1) = 1.121419727242096, (3, 2) = 0.7740148632607177e-1, (3, 3) = -4.207364448156272, (4, 1) = 1.1034190552999128, (4, 2) = -.3960219069400372, (4, 3) = -4.172548043626565, (5, 1) = 1.0334856968721362, (5, 2) = -.8551317943920755, (5, 3) = -4.042652033602019, (6, 1) = .9150351107388788, (6, 2) = -1.290241093967279, (6, 3) = -3.8421142104577206, (7, 1) = .7510656023540799, (7, 2) = -1.6983897118813762, (7, 3) = -3.6049293155861535, (8, 1) = .5426235137367472, (8, 2) = -2.082307479123934, (8, 3) = -3.3711644866560833, (9, 1) = .2896295135416169, (9, 2) = -2.447663710557477, (9, 3) = -3.1893238955632652, (10, 1) = .0, (10, 2) = -2.7947930788176802, (10, 3) = -3.1183169372581445}, datatype = float[8], order = C_order); YP := Matrix(10, 3, {(1, 1) = 1.0, (1, 2) = -3.983477327766739, (1, 3) = -1.737938145506, (2, 1) = .549346065085729, (2, 2) = -4.1393578797975605, (2, 3) = -1.0371772019479302, (3, 1) = 0.7740148632607177e-1, (3, 2) = -4.207364448156272, (3, 3) = -.15092384429708883, (4, 1) = -.3960219069400372, (4, 2) = -4.172548043626565, (4, 3) = .7597529210528466, (5, 1) = -.8551317943920755, (5, 2) = -4.042652033602019, (5, 3) = 1.5361495967515448, (6, 1) = -1.290241093967279, (6, 2) = -3.8421142104577206, (6, 3) = 2.0510023517096023, (7, 1) = -1.6983897118813762, (7, 2) = -3.6049293155861535, (7, 3) = 2.2136969872974745, (8, 1) = -2.082307479123934, (8, 2) = -3.3711644866560833, (8, 3) = 1.9564571436311755, (9, 1) = -2.447663710557477, (9, 2) = -3.1893238955632652, (9, 3) = 1.2188067920130512, (10, 1) = -2.7947930788176802, (10, 2) = -3.1183169372581445, (10, 3) = -0.21289086127144367e-1}, 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) = .0, (2) = .11078278220502231, (3) = .22364174273643112, (4) = .336400857639379, (5) = .44797504433316043, (6) = .5582094408957031, (7) = .6677790172685235, (8) = .7779200458343468, (9) = .8895341149733956, (10) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(10, 3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = -0.13603335696444332e-7, (2, 1) = 0.7265581332725971e-9, (2, 2) = -0.9179594754749238e-9, (2, 3) = -0.14360711048599372e-7, (3, 1) = 0.935850951586146e-9, (3, 2) = -0.33022746441561657e-8, (3, 3) = -0.15358678429518584e-7, (4, 1) = 0.1523118379956452e-9, (4, 2) = -0.5942457533255831e-8, (4, 3) = -0.14810819342344207e-7, (5, 1) = -0.13595390189336288e-8, (5, 2) = -0.7496861395347225e-8, (5, 3) = -0.13899207188268168e-7, (6, 1) = -0.29580647058515214e-8, (6, 2) = -0.7708912520795418e-8, (6, 3) = -0.15201085345852005e-7, (7, 1) = -0.40645381116283325e-8, (7, 2) = -0.72449801484776965e-8, (7, 3) = -0.20236709048260123e-7, (8, 1) = -0.42352021900713655e-8, (8, 2) = -0.6908010054283183e-8, (8, 3) = -0.28476365215460827e-7, (9, 1) = -0.29553323497707087e-8, (9, 2) = -0.6885055401174663e-8, (9, 3) = -0.3650784387748412e-7, (10, 1) = .0, (10, 2) = -0.6407492020611782e-8, (10, 3) = -0.3661757260073899e-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[10] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(3.661757260073899e-8) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [3, 10, [f0(x), diff(f0(x), x), diff(diff(f0(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(3, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(10, 3, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(3, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(10, 3, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[f0(x), diff(f0(x), x), diff(diff(f0(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[10] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(3.661757260073899e-8) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [3, 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(3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(10, 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`(10, 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, f0(x), diff(f0(x), x), diff(diff(f0(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('[f0(x), diff(f0(x), x), diff(diff(f0(x), x), x)]'[i] = res[i+1], i = 1 .. 3)] catch: error end try end proc$

(5)

F0 := proc(t)
  if not type(evalf(t),'numeric') then
    'procname'(t);
  else
    evalf(eval(f0(x), sol0(t)));
  end if;
end proc:

dF0 := proc(t)
  if not type(evalf(t),'numeric') then
    'procname'(t);
  else
    evalf(eval(eval(diff(f0(x), x), x=t), sol0(t)));
  end if;
end proc:

ddF0 := proc(t)
  if not type(evalf(t),'numeric') then
    'procname'(t);
  else
    evalf(eval(eval(diff(f0(x), x$2), x=t), sol0(t)));
  end if;
end proc:

subs0 := {f0(x)=F0(x), diff(f0(x), x)=dF0(x), diff(f0(x), x$2)=ddF0(x)};

{diff(diff(f0(x), x), x) = ddF0(x), diff(f0(x), x) = dF0(x), f0(x) = F0(x)}

(6)

First order correction f1

>	coeff(HO1, p, 1): eval(%, subs0); sol1 := dsolve({%, f1(0)=0, D(f1)(0)=0, f1(1)=0}, f1(x), numeric);

1.745555556*dF0(x)*f1(x)+1.745555556*(diff(f1(x), x))*F0(x)-1.745555556*dF0(x)*F0(x)+diff(diff(diff(f1(x), x), x), x)-0.7617410494e-2*(diff(f1(x), x))+0.7617410494e-2*dF0(x)

$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(16, {(1) = .0, (2) = 0.6404561717292101e-1, (3) = .13548250894342245, (4) = .2122716707745767, (5) = .2884888958349662, (6) = .3620331204527494, (7) = .4316113205671453, (8) = .4977368087758108, (9) = .5608265805477214, (10) = .6220917756931348, (11) = .6832353851897301, (12) = .7445626983123486, (13) = .8078204843197364, (14) = .8733810552722883, (15) = .9378567463688994, (16) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(16, 3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = -0.3212187661126076e-2, (2, 1) = 0.6564679961247358e-4, (2, 2) = 0.3114080432866802e-2, (2, 3) = 0.963683331251039e-1, (3, 1) = 0.6087451195508886e-3, (3, 2) = 0.13026771944825986e-1, (3, 3) = .17486744344098132, (4, 1) = 0.21753193364114107e-2, (4, 2) = 0.28279271091437002e-1, (4, 3) = .21419262565092856, (5, 1) = 0.4954626157255434e-2, (5, 2) = 0.44517132006358755e-1, (5, 3) = .20359601897522392, (6, 1) = 0.8738299477267725e-2, (6, 2) = 0.5768688395048585e-1, (6, 3) = .14728647498744749, (7, 1) = 0.1304213815879093e-1, (7, 2) = 0.6497183975446072e-1, (7, 3) = 0.5654829594084956e-1, (8, 1) = 0.17383420415696352e-1, (8, 2) = 0.6507721183392393e-1, (8, 3) = -0.57225106548140114e-1, (9, 1) = 0.212936602529797e-1, (9, 2) = 0.575497509540539e-1, (9, 3) = -.18367691820637957, (10, 1) = 0.2439299425024134e-1, (10, 2) = 0.4227794437310392e-1, (10, 3) = -.3157553158356348, (11, 1) = 0.26304429992423214e-1, (11, 2) = 0.18885964739179e-1, (11, 3) = -.4489745129134839, (12, 1) = 0.2653705961586595e-1, (12, 2) = -0.1260033744207179e-1, (12, 3) = -.5761147465647417, (13, 1) = 0.24506713309364235e-1, (13, 2) = -0.52817579643494386e-1, (13, 3) = -.6921816649495488, (14, 1) = 0.19482793601557198e-1, (14, 2) = -.10148682704173242, (14, 3) = -.7875818778401714, (15, 1) = 0.11253899896994142e-1, (15, 2) = -.15442304708314117, (15, 3) = -.8483267423312889, (16, 1) = .0, (16, 2) = -.2079779482830367, (16, 3) = -.8683725554267051}, datatype = float[8], order = C_order); YP := Matrix(16, 3, {(1, 1) = .0, (1, 2) = -0.3212187661126076e-2, (1, 3) = 1.737938145506, (2, 1) = 0.3114080432866802e-2, (2, 2) = 0.963683331251039e-1, (2, 3) = 1.3552081584558717, (3, 1) = 0.13026771944825986e-1, (3, 2) = .17486744344098132, (3, 3) = .8278280889766552, (4, 1) = 0.28279271091437002e-1, (4, 2) = .21419262565092856, (4, 3) = .1883677603012989, (5, 1) = 0.44517132006358755e-1, (5, 2) = .20359601897522392, (5, 3) = -.46447601588938187, (6, 1) = 0.5768688395048585e-1, (6, 2) = .14728647498744749, (6, 3) = -1.056149999809692, (7, 1) = 0.6497183975446072e-1, (7, 2) = 0.5654829594084956e-1, (7, 3) = -1.535754685752323, (8, 1) = 0.6507721183392393e-1, (8, 2) = -0.57225106548140114e-1, (8, 3) = -1.8859545164871088, (9, 1) = 0.575497509540539e-1, (9, 2) = -.18367691820637957, (9, 3) = -2.102092604529736, (10, 1) = 0.4227794437310392e-1, (10, 2) = -.3157553158356348, (10, 3) = -2.188630853345391, (11, 1) = 0.18885964739179e-1, (11, 2) = -.4489745129134839, (11, 3) = -2.1475052546467324, (12, 1) = -0.1260033744207179e-1, (12, 2) = -.5761147465647417, (12, 3) = -1.9774487568282264, (13, 1) = -0.52817579643494386e-1, (13, 2) = -.6921816649495488, (13, 3) = -1.6703846603530887, (14, 1) = -.10148682704173242, (14, 2) = -.7875818778401714, (14, 3) = -1.2177366978634918, (15, 1) = -.15442304708314117, (15, 2) = -.8483267423312889, (15, 3) = -.6461131491510258, (16, 1) = -.2079779482830367, (16, 2) = -.8683725554267051, (16, 3) = 0.19704832721372575e-1}, 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(16, {(1) = .0, (2) = 0.6404561717292101e-1, (3) = .13548250894342245, (4) = .2122716707745767, (5) = .2884888958349662, (6) = .3620331204527494, (7) = .4316113205671453, (8) = .4977368087758108, (9) = .5608265805477214, (10) = .6220917756931348, (11) = .6832353851897301, (12) = .7445626983123486, (13) = .8078204843197364, (14) = .8733810552722883, (15) = .9378567463688994, (16) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(16, 3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = 0.7134055481678411e-9, (2, 1) = -0.21932230435644282e-10, (2, 2) = 0.4590245883148544e-11, (2, 3) = 0.7461923959253619e-9, (3, 1) = -0.6082398462942408e-10, (3, 2) = 0.3065326957403016e-10, (3, 3) = 0.8683866891000823e-9, (4, 1) = -0.8695065362407275e-10, (4, 2) = 0.11153256807576377e-9, (4, 3) = 0.9810666598014438e-9, (5, 1) = -0.6053555201500314e-10, (5, 2) = 0.20181250467793794e-9, (5, 3) = 0.9251916928967297e-9, (6, 1) = -0.41073477234674456e-12, (6, 2) = 0.2495361553785567e-9, (6, 3) = 0.7704948562363915e-9, (7, 1) = 0.6056818024074013e-10, (7, 2) = 0.2515251688390155e-9, (7, 3) = 0.6439500967633543e-9, (8, 1) = 0.10868684843511718e-9, (8, 2) = 0.2301904357179198e-9, (8, 3) = 0.5936739598428944e-9, (9, 1) = 0.14076986220351994e-9, (9, 2) = 0.20625579736177594e-9, (9, 3) = 0.623410717783828e-9, (10, 1) = 0.15987175494471696e-9, (10, 2) = 0.1900168529372196e-9, (10, 3) = 0.7165453504154215e-9, (11, 1) = 0.16825746824195036e-9, (11, 2) = 0.18477923694672967e-9, (11, 3) = 0.8650304385793945e-9, (12, 1) = 0.16514866680272493e-9, (12, 2) = 0.19808534216613374e-9, (12, 3) = 0.10595450219154226e-8, (13, 1) = 0.14657074566307816e-9, (13, 2) = 0.2315425152103885e-9, (13, 3) = 0.13008516723077958e-8, (14, 1) = 0.10422197575497192e-9, (14, 2) = 0.29065656833632424e-9, (14, 3) = 0.1564887441917633e-8, (15, 1) = 0.5088008014332878e-10, (15, 2) = 0.3817429595284989e-9, (15, 3) = 0.17589401418123942e-8, (16, 1) = .0, (16, 2) = 0.4863562698983413e-9, (16, 3) = 0.18476646867013011e-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[16] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.8476646867013011e-9) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [3, 16, [f1(x), diff(f1(x), x), diff(diff(f1(x), x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[16] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[16] 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`(16, 3, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(3, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(16, 3, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[f1(x), diff(f1(x), x), diff(diff(f1(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[16] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.8476646867013011e-9) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [3, 16, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[16] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[16] 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`(16, 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`(16, 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, f1(x), diff(f1(x), x), diff(diff(f1(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('[f1(x), diff(f1(x), x), diff(diff(f1(x), x), x)]'[i] = res[i+1], i = 1 .. 3)] catch: error end try end proc$

(7)

>	display( odeplot(sol0, [x, f0(x)], x=0..1, color=blue) , odeplot(sol1, [x, f1(x)+F0(x)], x=0..1, color=purple) )

F1 := proc(t)
  if not type(evalf(t),'numeric') then
    'procname'(t);
  else
    evalf(eval(f1(x), sol1(t)));
  end if;
end proc:

dF1 := proc(t)
  if not type(evalf(t),'numeric') then
    'procname'(t);
  else
    evalf(eval(eval(diff(f1(x), x), x=t), sol1(t)));
  end if;
end proc:

ddF1 := proc(t)
  if not type(evalf(t),'numeric') then
    'procname'(t);
  else
    evalf(eval(eval(diff(f1(x), x$2), x=t), sol1(t)));
  end if;
end proc:

subs1 := subs0 union {f1(x)=F1(x), diff(f1(x), x)=dF1(x), diff(f1(x), x$2)=ddF1(x)};

{diff(diff(f0(x), x), x) = ddF0(x), diff(diff(f1(x), x), x) = ddF1(x), diff(f0(x), x) = dF0(x), diff(f1(x), x) = dF1(x), f0(x) = F0(x), f1(x) = F1(x)}

(8)

Second order correction f2

>	coeff(HO1, p, 2): eval(%, isolate(coeff(HO1, p, 0), diff(f0(x), x$3))): eval(%, subs1); sol2 := dsolve({%, f2(0)=0, D(f2)(0)=0, f2(1)=0}, f2(x), numeric);

1.745555556*dF0(x)*f2(x)+1.745555556*dF1(x)*F1(x)+1.745555556*(diff(f2(x), x))*F0(x)-1.745555556*dF0(x)*F1(x)-1.745555556*dF1(x)*F0(x)+diff(diff(diff(f2(x), x), x), x)+0.7617410494e-2*dF1(x)-0.7617410494e-2*dF0(x)+1.745555556*dF0(x)*F0(x)-0.7617410494e-2*(diff(f2(x), x))

$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(16, {(1) = .0, (2) = 0.6401730673686455e-1, (3) = .13535434866443885, (4) = .21184479539811732, (5) = .2874106622391402, (6) = .36018298833690005, (7) = .4289747704390974, (8) = .4945158762055879, (9) = .5573544786848745, (10) = .6185541600774694, (11) = .6799679578806432, (12) = .7418691804331442, (13) = .805878954011701, (14) = .8724424641402733, (15) = .9376435075095818, (16) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(16, 3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = -0.10855106714429266e-1, (2, 1) = -0.9434324559718517e-4, (2, 2) = -0.4008921457300064e-2, (2, 3) = -.11022411295479144, (3, 1) = -0.7345226618178199e-3, (3, 2) = -0.14858406176389628e-1, (3, 3) = -.1874264034029753, (4, 1) = -0.24668572303698107e-2, (4, 2) = -0.30896215937434624e-1, (4, 3) = -.22337792269059026, (5, 1) = -0.5436523078143244e-2, (5, 2) = -0.4750302824001485e-1, (5, 3) = -.20761400745680145, (6, 1) = -0.939808779150627e-2, (6, 2) = -0.6061808070998779e-1, (6, 3) = -.14550383128999295, (7, 1) = -0.13843203706555673e-1, (7, 2) = -0.6751413435805735e-1, (7, 3) = -0.4954511550414553e-1, (8, 1) = -0.182941169975788e-1, (8, 2) = -0.6701748576587586e-1, (8, 3) = 0.6830228433029469e-1, (9, 1) = -0.22287523582588977e-1, (9, 2) = -0.5872929299610989e-1, (9, 3) = .19740923210996972, (10, 1) = -0.2542968960358908e-1, (10, 2) = -0.42601624432994124e-1, (10, 3) = .3300957365992241, (11, 1) = -0.27339849032211847e-1, (11, 2) = -0.18253967782892542e-1, (11, 3) = .46191475616382505, (12, 1) = -0.27503867881480085e-1, (12, 2) = 0.14230167062005505e-1, (12, 3) = .5854310388105278, (13, 1) = -0.25314865873590324e-1, (13, 2) = 0.5534105024997991e-1, (13, 3) = .6954860742794319, (14, 1) = -0.2001989833588037e-1, (14, 2) = .10473118863537298, (14, 3) = .783448101559649, (15, 1) = -0.11481766747797573e-1, (15, 2) = .15775651294698934, (15, 3) = .8372939488619414, (16, 1) = .0, (16, 2) = .21068302134443742, (16, 3) = .8543052611659095}, datatype = float[8], order = C_order); YP := Matrix(16, 3, {(1, 1) = .0, (1, 2) = -0.10855106714429266e-1, (1, 3) = -1.737938145506, (2, 1) = -0.4008921457300064e-2, (2, 2) = -.11022411295479144, (2, 3) = -1.3479189615379026, (3, 1) = -0.14858406176389628e-1, (3, 2) = -.1874264034029753, (3, 3) = -.7999164710320333, (4, 1) = -0.30896215937434624e-1, (4, 2) = -.22337792269059026, (4, 3) = -.13140940246672603, (5, 1) = -0.4750302824001485e-1, (5, 2) = -.20761400745680145, (5, 3) = .5455239983517919, (6, 1) = -0.6061808070998779e-1, (6, 2) = -.14550383128999295, (6, 3) = 1.1483341496355746, (7, 1) = -0.6751413435805735e-1, (7, 2) = -0.4954511550414553e-1, (7, 3) = 1.6223192555783386, (8, 1) = -0.6701748576587586e-1, (8, 2) = 0.6830228433029469e-1, (8, 3) = 1.951686014460034, (9, 1) = -0.5872929299610989e-1, (9, 2) = .19740923210996972, (9, 3) = 2.13471956047877, (10, 1) = -0.42601624432994124e-1, (10, 2) = .3300957365992241, (10, 3) = 2.17926647087871, (11, 1) = -0.18253967782892542e-1, (11, 2) = .46191475616382505, (11, 3) = 2.09200960450085, (12, 1) = 0.14230167062005505e-1, (12, 2) = .5854310388105278, (12, 3) = 1.8788522008739734, (13, 1) = 0.5534105024997991e-1, (13, 2) = .6954860742794319, (13, 3) = 1.5415600873402, (14, 1) = .10473118863537298, (14, 2) = .783448101559649, (14, 3) = 1.0852525314768466, (15, 1) = .15775651294698934, (15, 2) = .8372939488619414, (15, 3) = .5540197824891389, (16, 1) = .21068302134443742, (16, 2) = .8543052611659095, (16, 3) = -0.1809997366367583e-1}, 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(16, {(1) = .0, (2) = 0.6401730673686455e-1, (3) = .13535434866443885, (4) = .21184479539811732, (5) = .2874106622391402, (6) = .36018298833690005, (7) = .4289747704390974, (8) = .4945158762055879, (9) = .5573544786848745, (10) = .6185541600774694, (11) = .6799679578806432, (12) = .7418691804331442, (13) = .805878954011701, (14) = .8724424641402733, (15) = .9376435075095818, (16) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(16, 3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = -0.1109561503003624e-8, (2, 1) = 0.26699812435622908e-10, (2, 2) = -0.19918851259478768e-10, (2, 3) = -0.11812785942632816e-8, (3, 1) = 0.7290545265482255e-10, (3, 2) = -0.435518863276401e-10, (3, 3) = -0.1370038408939214e-8, (4, 1) = 0.10327208677955919e-9, (4, 2) = -0.9099775069557878e-10, (4, 3) = -0.15286310683797769e-8, (5, 1) = 0.738457274020005e-10, (5, 2) = -0.14594628936276227e-9, (5, 3) = -0.14424793442283092e-8, (6, 1) = 0.9057897962753132e-11, (6, 2) = -0.17240836061534957e-9, (6, 3) = -0.12205996244564129e-8, (7, 1) = -0.54425291935501215e-10, (7, 2) = -0.17098292347539843e-9, (7, 3) = -0.10323622296221041e-8, (8, 1) = -0.10351756071377093e-9, (8, 2) = -0.15663723131857096e-9, (8, 3) = -0.9356588405429689e-9, (9, 1) = -0.13549220262582118e-9, (9, 2) = -0.1452319944728306e-9, (9, 3) = -0.9358672788267563e-9, (10, 1) = -0.15331699329356284e-9, (10, 2) = -0.14588951928284822e-9, (10, 3) = -0.10145347686655287e-8, (11, 1) = -0.15964518339726942e-9, (11, 2) = -0.16289058555461386e-9, (11, 3) = -0.1162024809582045e-8, (12, 1) = -0.15341111964548522e-9, (12, 2) = -0.20676190160146303e-9, (12, 3) = -0.13701166123373674e-8, (13, 1) = -0.1306663509042881e-9, (13, 2) = -0.28766725617611577e-9, (13, 3) = -0.16404736621639975e-8, (14, 1) = -0.839918600087215e-10, (14, 2) = -0.4214750711551729e-9, (14, 3) = -0.1947795993242128e-8, (15, 1) = -0.3463981824448189e-10, (15, 2) = -0.5968250300218129e-9, (15, 3) = -0.21655810926463017e-8, (16, 1) = .0, (16, 2) = -0.780446270962002e-9, (16, 3) = -0.2262710129722444e-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[16] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(2.262710129722444e-9) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [3, 16, [f2(x), diff(f2(x), x), diff(diff(f2(x), x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[16] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[16] 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`(16, 3, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(3, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(16, 3, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[f2(x), diff(f2(x), x), diff(diff(f2(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[16] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(2.262710129722444e-9) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [3, 16, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[16] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[16] 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`(16, 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`(16, 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, f2(x), diff(f2(x), x), diff(diff(f2(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('[f2(x), diff(f2(x), x), diff(diff(f2(x), x), x)]'[i] = res[i+1], i = 1 .. 3)] catch: error end try end proc$

(9)

>	display( odeplot(sol0, [x, f0(x)], x=0..1, color=blue) , odeplot(sol1, [x, f1(x)+F0(x)], x=0..1, color=purple) , odeplot(sol2, [x, f2(x)+F0(x)+F1(x)], x=0..1, color=red) )

F2 := proc(t)
  if not type(evalf(t),'numeric') then
    'procname'(t);
  else
    evalf(eval(f2(x), sol2(t)));
  end if;
end proc:

dF2 := proc(t)
  if not type(evalf(t),'numeric') then
    'procname'(t);
  else
    evalf(eval(eval(diff(f2(x), x), x=t), sol2(t)));
  end if;
end proc:

ddF2 := proc(t)
  if not type(evalf(t),'numeric') then
    'procname'(t);
  else
    evalf(eval(eval(diff(f2(x), x$2), x=t), sol2(t)));
  end if;
end proc:

subs2 := subs1 union {f2(x)=F2(x), diff(f2(x), x)=dF2(x), diff(f2(x), x$2)=ddF2(x)};

{diff(diff(f0(x), x), x) = ddF0(x), diff(diff(f1(x), x), x) = ddF1(x), diff(diff(f2(x), x), x) = ddF2(x), diff(f0(x), x) = dF0(x), diff(f1(x), x) = dF1(x), diff(f2(x), x) = dF2(x), f0(x) = F0(x), f1(x) = F1(x), f2(x) = F2(x)}

(10)

Third order correction f3

>	coeff(HO1, p, 3): eval(%, isolate(coeff(HO1, p, 0), diff(f0(x), x$3))): eval(%, isolate(coeff(HO1, p, 1), diff(f1(x), x$3))); eval(%, subs2); sol3 := dsolve({%, f3(0)=0, D(f3)(0)=0, f3(1)=0}, f3(x), numeric);

1.745555556*dF0(x)*f3(x)+1.745555556*dF1(x)*F2(x)+1.745555556*dF2(x)*F1(x)+1.745555556*(diff(f3(x), x))*F0(x)-1.745555556*dF0(x)*F2(x)-1.745555556*dF1(x)*F1(x)-1.745555556*dF2(x)*F0(x)+diff(diff(diff(f3(x), x), x), x)+0.7617410494e-2*dF2(x)+1.745555556*dF0(x)*F1(x)+1.745555556*dF1(x)*F0(x)-1.745555556*dF0(x)*F0(x)-0.7617410494e-2*dF1(x)+0.7617410494e-2*dF0(x)-0.7617410494e-2*(diff(f3(x), x))

$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(15, {(1) = .0, (2) = 0.6868100819699492e-1, (3) = .14595638898984406, (4) = .22740500925807855, (5) = .30717748303792447, (6) = .38288700335472736, (7) = .4543480349020496, (8) = .5227453148711261, (9) = .5886887155389848, (10) = .6545230333992893, (11) = .7211386562635443, (12) = .7898630878316497, (13) = .8618912771704792, (14) = .9329501707442724, (15) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(15, 3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = 0.25694532785562645e-1, (2, 1) = 0.1491788703088783e-3, (2, 2) = 0.5551372935008765e-2, (2, 3) = .13099940113459047, (3, 1) = 0.10582537588995147e-2, (3, 2) = 0.1897376951598749e-1, (3, 3) = .2082481619796001, (4, 1) = 0.3338407184419461e-2, (4, 2) = 0.37348640764275555e-1, (4, 3) = .23267324550810226, (5, 1) = 0.7035325223644292e-2, (5, 2) = 0.54854993962537416e-1, (5, 3) = .19642015019561962, (6, 1) = 0.11680425308889414e-1, (6, 2) = 0.66768337313347e-1, (6, 3) = .11057958569318455, (7, 1) = 0.16639899614274103e-1, (7, 2) = 0.7061961064737002e-1, (7, 3) = -0.797137778600824e-2, (8, 1) = 0.21348091436839207e-1, (8, 2) = 0.6548908832301815e-1, (8, 3) = -.14490233927517804, (9, 1) = 0.25249173230482803e-1, (9, 2) = 0.5125732753073486e-1, (9, 3) = -.2875375207494506, (10, 1) = 0.27897092799920892e-1, (10, 2) = 0.27626301543563176e-1, (10, 3) = -.42938749383968755, (11, 1) = 0.2868373948753689e-1, (11, 2) = -0.5483361420010717e-2, (11, 3) = -.5620492132282962, (12, 1) = 0.26883334437606946e-1, (12, 2) = -0.4824239819134363e-1, (12, 3) = -.6781167048832408, (13, 1) = 0.21563292898862964e-1, (13, 2) = -.10057345513347787, (13, 3) = -.7692213545779051, (14, 1) = 0.12421270049959098e-1, (14, 2) = -.1573760517562036, (14, 3) = -.8233191152894311, (15, 1) = .0, (15, 2) = -.21331312092771168, (15, 3) = -.8394658359534058}, datatype = float[8], order = C_order); YP := Matrix(15, 3, {(1, 1) = .0, (1, 2) = 0.25694532785562645e-1, (1, 3) = 1.737938145506, (2, 1) = 0.5551372935008765e-2, (2, 2) = .13099940113459047, (2, 3) = 1.304932394663431, (3, 1) = 0.1897376951598749e-1, (3, 2) = .2082481619796001, (3, 3) = .6737940929200968, (4, 1) = 0.37348640764275555e-1, (4, 2) = .23267324550810226, (4, 3) = -0.8182620950393095e-1, (5, 1) = 0.54854993962537416e-1, (5, 2) = .19642015019561962, (5, 3) = -.8183665285133928, (6, 1) = 0.66768337313347e-1, (6, 2) = .11057958569318455, (6, 3) = -1.4284164631328513, (7, 1) = 0.7061961064737002e-1, (7, 2) = -0.797137778600824e-2, (7, 3) = -1.8625610701017996, (8, 1) = 0.6548908832301815e-1, (8, 2) = -.14490233927517804, (8, 3) = -2.112540126710141, (9, 1) = 0.5125732753073486e-1, (9, 2) = -.2875375207494506, (9, 3) = -2.185773121994767, (10, 1) = 0.27626301543563176e-1, (10, 2) = -.42938749383968755, (10, 3) = -2.0977739473107095, (11, 1) = -0.5483361420010717e-2, (11, 2) = -.5620492132282962, (11, 3) = -1.8627989688402546, (12, 1) = -0.4824239819134363e-1, (12, 2) = -.6781167048832408, (12, 3) = -1.497097500557855, (13, 1) = -.10057345513347787, (13, 2) = -.7692213545779051, (13, 3) = -1.0202162547215374, (14, 1) = -.1573760517562036, (14, 2) = -.8233191152894311, (14, 3) = -.49670159106503564, (15, 1) = -.21331312092771168, (15, 2) = -.8394658359534058, (15, 3) = 0.16475080057813184e-1}, 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(15, {(1) = .0, (2) = 0.6868100819699492e-1, (3) = .14595638898984406, (4) = .22740500925807855, (5) = .30717748303792447, (6) = .38288700335472736, (7) = .4543480349020496, (8) = .5227453148711261, (9) = .5886887155389848, (10) = .6545230333992893, (11) = .7211386562635443, (12) = .7898630878316497, (13) = .8618912771704792, (14) = .9329501707442724, (15) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(15, 3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = 0.22389404042291576e-8, (2, 1) = -0.5168218350433026e-10, (2, 2) = 0.57468118190898715e-10, (2, 3) = 0.24344806834023884e-8, (3, 1) = -0.1417976448449938e-9, (3, 2) = 0.8404174699561922e-10, (3, 3) = 0.28914120499479665e-8, (4, 1) = -0.17738113188552843e-9, (4, 2) = 0.10287376876120677e-9, (4, 3) = 0.3134948950104668e-8, (5, 1) = -0.10327528482632357e-9, (5, 2) = 0.1263970472149691e-9, (5, 3) = 0.28325343046744072e-8, (6, 1) = 0.10202960090514321e-10, (6, 2) = 0.1334652921831637e-9, (6, 3) = 0.23559724843213524e-8, (7, 1) = 0.1066953549222101e-9, (7, 2) = 0.12847849228984948e-9, (7, 3) = 0.2010184323492844e-8, (8, 1) = 0.17396091548402068e-9, (8, 2) = 0.12586648554756268e-9, (8, 3) = 0.18599761019367305e-8, (9, 1) = 0.21124400605476557e-9, (9, 2) = 0.14579545969855338e-9, (9, 3) = 0.18947527456398552e-8, (10, 1) = 0.22654582423044557e-9, (10, 2) = 0.197432459382499e-9, (10, 3) = 0.20815927918919977e-8, (11, 1) = 0.2185432085017301e-9, (11, 2) = 0.3023159295450315e-9, (11, 3) = 0.24129148737532237e-8, (12, 1) = 0.18226111479919648e-9, (12, 2) = 0.4902954681979294e-9, (12, 3) = 0.28885241734514984e-8, (13, 1) = 0.10489511429189739e-9, (13, 2) = 0.8083964868359335e-9, (13, 3) = 0.3489441941514877e-8, (14, 1) = 0.28291286689841743e-10, (14, 2) = 0.1218372130067739e-8, (14, 3) = 0.3933864731719928e-8, (15, 1) = .0, (15, 2) = 0.16233326858211381e-8, (15, 3) = 0.41201977518441975e-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[15] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(4.1201977518441975e-9) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [3, 15, [f3(x), diff(f3(x), x), diff(diff(f3(x), x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[15] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[15] 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`(15, 3, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(3, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(15, 3, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[f3(x), diff(f3(x), x), diff(diff(f3(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[15] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(4.1201977518441975e-9) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [3, 15, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[15] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[15] 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`(15, 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`(15, 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, f3(x), diff(f3(x), x), diff(diff(f3(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('[f3(x), diff(f3(x), x), diff(diff(f3(x), x), x)]'[i] = res[i+1], i = 1 .. 3)] catch: error end try end proc$

(11)

>	display( odeplot(sol0, [x, f0(x)], x=0..1, color=blue) , odeplot(sol1, [x, f1(x)+F0(x)], x=0..1, color=purple) , odeplot(sol2, [x, f2(x)+F0(x)+F1(x)], x=0..1, color=red) , odeplot(sol3, [x, f3(x)+F0(x)+F1(x)+F2(x)], x=0..1, color=gold) )

and so on...

Plot corrections:

(note that the process doesn't seem to converge)

>	display( odeplot(sol1, [x, f1(x)], x=0..1, color=purple) , odeplot(sol2, [x, f2(x)], x=0..1, color=red) , odeplot(sol3, [x, f3(x)], x=0..1, color=gold) )

Download HPM_1_iterative.mw

A few remarks...

Answered: sand15 615

July 15 2023

0 2

Remark 1
Your initial pde must contain partial derivatices wrt x.
When applying the Fourier transform the imaginary root I should appear elsewhere, unless all the partial derivatives wrt x. are of the even orders.
See here

edo := M*diff(u(x), x$2)+C*diff(u(x), x)+K*u(x)=f(x);
FT := inttrans:-fourier(edo, x, k):

ToReplace := remove(has, convert(select(has, indets(%, function), fourier), list), diff):
ReplaceBy := map(i -> parse(StringTools:-Capitalize(op(0, i)))(k), map2(op, 1, ToReplace)):

Fode := eval(FT, ToReplace =~ ReplaceBy);

M*(diff(diff(u(x), x), x))+C*(diff(u(x), x))+K*u(x) = f(x)

(1)

pde := C*diff(u(x, t), t)+lambda*diff(u(x, t), x$2) = f(x, t);
FT := inttrans:-fourier(pde, x, k);

ToReplace := remove(has, convert(select(has, indets(%, function), fourier), list), diff):
ReplaceBy := map(i -> parse(StringTools:-Capitalize(op(0, i)))(k, t), map2(op, 1, ToReplace)):

Fpde := eval(FT, ToReplace =~ ReplaceBy):

(2)

pde := C*diff(u(x, t), t)+v*diff(u(x, t), x) = f(x, t);
FT := inttrans:-fourier(pde, x, k):

ToReplace := remove(has, convert(select(has, indets(%, function), fourier), list), diff):
ReplaceBy := map(i -> parse(StringTools:-Capitalize(op(0, i)))(k, t), map2(op, 1, ToReplace)):

Fpde := eval(FT, ToReplace =~ ReplaceBy);

(3)

Download Remark_1.mw

Did you omit I or is there really no I?

Remark 2
Given the complexity of the ode in the Fourier space, I guess the initial pde contains terms zeta(x,t)^r where r is some integer.
Applying the Fourier transform to a term like this should introduce a convolution product (... unless if the initial pde contains itself convolution products which are transformed into ordinary products by the Fourier transform?)

Download Remark_2.mw

Remark 3
As your ode is extremely complex, are you suresolving directly the initial pde (in a numeric way of course) is that bad an idea?

Remark 4
Fourier transform has the advantage to transform a linear pde into a simple ode, but you must not omit the boundary condions.
This help page explain how to solve formally an ode with some integral transformations

help(dsolve[inttrans])

What is the best way to solve this rrobl...

Answered: sand15 615

July 10 2023

1 2

What is the best way to solve this rroblem? (apart from the obvious choice to fit the data and solve the equation with the fitted potential)

Why do you say that? Do you think this "obvious choice" doesn't deserve to be taken into account?

I guess that V(x) comes from the numerical solution of an ODE, which doesn't contain y(x) nor any of its derivative?

If it's not the case just solve a coupled system of odes.

If it is the case I advice you to look this help page

help(dsolve[numeric,efficiency])

and go to the section Solutions depending on other solutions.
You will find an example where a first ode is solved numerically to get a procedure

                  xf:=proc(t) ... end proc

(typically xf could be V) and then this procedure is plugged into a second ode system (typically the -y''(x) = (E-V(x))*y(x) ode) :

sys2 := {diff(y(t),t,t)=xf(t), y(0)=0, D(y)(0)=1};
sol2 := dsolve(sys2, numeric, known={xf});
                    sol2:=proc(x_rkf45) ... end proc

So you could write something like this

sys2 := {-diff(y(x), x$2) = (E-V(x))*y(x), +IC/BC};
sol2 := dsolve(sys2, numeric, known={V});

This avoids the "fitting step" which you do not like.

Depending of the complexity of the ode V(x) verifies and on the complexity of V(x), the "fitting" approach can perform better, or be worse, than the metof describe in the help page.

Expert's answer is required to confirm m...

Answered: sand15 615

June 06 2023

3 2

f1 can be integrated formally, but it seems f2 cannot.
(the assumptions I use come from the kc and kh ranges you give for f1).

restart:

>	# # for numerical integration, assume sigma=10, kc=1.6..3, kh=3..7, deltac=1-deltah , deltah=0.65 # f1:=deltac(1-x/(kc-3/2))^(-kc+1/2)+deltah(1-x/(sigma*(kh-3/2)))^(-kh+1/2);

deltac*(1-x/(kc-3/2))^(-kc+1/2)+deltah*(1-x/(sigma*(kh-3/2)))^(-kh+1/2)

(1)

There is no need to use a numerical integration for f1

# For some unknown reason
# int(f1, x=0..X) assuming kc > 1/2, kh > 1/2, X > 0;
# vives no answer with Maple 2015.2 after a reasonable amount of time.
#
# So this trick:

int(f1, x) assuming kc > 1/2, kh > 1/2;
If1 := unapply(eval(%, x=X)-eval(%, x=0), [X, sigma, kc, kh, deltac, deltah]):

plot3d(If1(0.15, 10, kc, kh, 1-0.65, 0.65), kc=1.6..3, kh=3..7, labels=["kc", "kh", "int"], style=surface)

-deltac*(1-x/(kc-3/2))^(-kc+3/2)*(kc-3/2)/(-kc+3/2)-deltah*(1-x/(sigma*(kh-3/2)))^(-kh+3/2)*sigma*(kh-3/2)/(-kh+3/2)

I guess numerical integration is needed for f2

>	f2 := (deltac(1-x/(kc-3/2))^(-kc+1/2)+deltah(1-x/(sigma(kh-3/2)))^(-kh+1/2)) (deltac(1-x/(kc-3/2))^(-kc+3/2)+deltahsigma(1-x/(sigma(kh-3/2)))^(-kh+3/2))

(deltac*(1-x/(kc-3/2))^(-kc+1/2)+deltah*(1-x/(sigma*(kh-3/2)))^(-kh+1/2))*(deltac*(1-x/(kc-3/2))^(-kc+3/2)+deltah*sigma*(1-x/(sigma*(kh-3/2)))^(-kh+3/2))

(2)

>	with(IntegrationTools):

>	If2 := Expand(Int(f2, x)) assuming kc > 1/2, kh > 3/2: if2 := value(If2);

(3)

>	select(has, {op(if2)}, int)

${deltah*deltac*(int((1-x/(sigma*(kh-3/2)))^(1/2)*(1-x/(kc-3/2))^(3/2)/((1-x/(sigma*(kh-3/2)))^kh*(1-x/(kc-3/2))^kc), x)), deltac*deltah*sigma*(int((1-x/(kc-3/2))^(1/2)*(1-x/(sigma*(kh-3/2)))^(3/2)/((1-x/(kc-3/2))^kc*(1-x/(sigma*(kh-3/2)))^kh), x))}$

(4)

>	(1-x/b)^(1/2-kh); # and int(expr, x)

(1-x/a)^(-kc+1/2)*(1-x/b)^(-kh+1/2)

int((1-x/a)^(-kc+1/2)*(1-x/b)^(-kh+1/2), x)

(5)

Download f1_and_f2.mw

The simple command res := eval(<vars&...

Answered: sand15 615

May 20 2023

0 0

The simple command

res := eval(<vars>, fsolve(eval({e||(1..8)}, [y1[n]=iny1, y2[n]=iny2]), {vars})):

works perfectly well and your code produces something which seems correct.
Did you think that tolerance is an option of eval? In this cas the answer is NO: eval has only 1 or 2 arguments (look to the corresponding help page)

Had it been the case, the lines

tolerance := 1e-6:
res := eval(<vars>, fsolve(eval({e||(1..8)}, [y1[n]=iny1, y2[n]=iny2]), {vars}), tolerance = tolerance):

would have produced an error; for instance:

restart
dsys := {diff(x(t),t)=y(t),diff(y(t),t)=-x(t),x(0)=1,y(0)=0}:
abserr := 1e-6;
                            0.000001
dsol := dsolve(dsys, numeric, abserr=abserr):
Error, (in dsolve/numeric/an_args/SC) keyword was '0.1e-5', optional keyword must be one of 'abserr', 
'delaymax', 'delaypts', 'differential', 'event_doublecross', 'event_initial', 'event_iterate', 
'event_maxiter', 'event_pre', 'event_project', 'event_relrange', 'event_stepreduction', 'events', 
'implicit', 'initstep', 'interpolate', 'interr', 'maxfun', 'maxstep', 'minstep', 'optimize', 
'output', 'projection', 'range', 'relerr', 'startinit', 'steppast'

To avaoid this error the correct syntax is (note the simple quotes arround the first occurence of abserr)

dsol := dsolve(dsys, numeric, 'abserr'=abserr):

1 2 3 4

Page 1 of 4

Share via:

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

sand15

615 Reputation

11 Badges

MaplePrimes Activity

These are answers submitted by sand15

Have you any idea of the size of Bd?...

I would say no...

Computation of the determinant of B and ...

What loop are you talking about? A_...

The main problem with Word is that the c...

You already found the solution...

After diffr2 and diffr3 corrected...

A way...

@Rana47 The solution I get (note th...

Still need help?...

After having corrected the (likely) typo...

A few remarks...

What is the best way to solve this rrobl...

Expert's answer is required to confirm m...

The simple command res := eval(<vars&...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again