## 615 Reputation

8 years, 311 days

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

```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

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

So how can you even read this expression?

## I would say no...

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 ...

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

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)
 >
 >
 (2)
 >
 >
 >

 > t0 := time(): detH  := Determinant(H): detB1 := Determinant(B1): detB  := detH * detB1: time()-t0;
 (3)
 > simplify(detH, size); simplify(detB1, size)
 (4)
 >
 >

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_...

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...

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"

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

## You already found the solution...

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

 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 > 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])
 >

## After diffr2 and diffr3 corrected...

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

 >
 >
 (1)
 >
 (2)
 >
 (3)
 >
 >
 >
 (4)
 >
 (5)

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

@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)+n*b*(diff(g(x), x))^(n-1)*(diff(g(x), x\$2))-k)
 (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-2*x)+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)
 (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") )
 >

## Still need help?...

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))+2*alpha*RE*(diff(f(x), x))*f(x)+(4-H)*alpha^2*(diff(f(x), x)))
 (2)
 > HPMEq1 := eval(HPMEq0, ansatz)
 (3)
 > for i from 0 to N do   equ[1][i] := coeff(HPMEq1, p, i) = 0 end do;
 (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)
 (6)
 > H := 5; alpha := (1/180)*(5*3.142); RE := 10
 (7)
 (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" )
 >

## After having corrected the (likely) typo...

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

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)))+2*alpha*RE*(diff(f(x), x))*f(x)+(4-H)*alpha^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);
 (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)};
 (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);
 (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)};
 (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);
 (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)};
 (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);
 (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) )
 >

## A few remarks...

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);
 (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)
 >

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?)

 > inttrans:-fourier(u(x)*u(x), x, k) = int(U(kappa)*U(k-kappa), kappa=-infinity..+infinity)
 (1)
 >

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...

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...

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);
 (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)

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)+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)
 (4)
 > (1-x/b)^(1/2-kh); # and int(expr, x)
 (5)
 >

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

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
﻿