## 787 Reputation

9 years, 177 days

## 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):`

## Another point of view: remove singularit...

Both the rhs have potentially 2 singularities: one at R=0, the other at R=1.
Computing limit(rhs(Sys[2]), R=1) returns alpha^2 and thus there is no problem at R=1 for the Phi-ode.

To resolve the 3 remaining singularities I propose to replace localy the rhs by its series expansion.
More precisely:

• T-ode:
• replace the rhs at R=0 by a series expansion T_SL in the interval [0, h],
• note T_CS the rhs in the interval [h, 1-h],
• replace the rhas at R=0 by a series expansion T_SR in the interval [1-h, L), where L is any number > 1,
(possibly L = 1+h and a fourth range [1+h, L) could be included where T_CS2 is equal to the rhs)
• dsolve exactly {diff(T(R), R) = T_LS, T(0)=0.5},
• dsolve numerically {diff(T(R), R) = T_CS, T(h)=U} where U is the value of T at point R=h for the previous solution,
• dsolve exactly {diff(T(R), R) = T_RS, T(1-h)=V} where U is the value of T at point R=h for the previous solution.
(possibly dsolve numerically in the fourth range)
• Phi-ode
• replace the rhs at R=0 by a series expansion Phi_SL in the interval [0, h],
• note T_CS the rhs in the interval [h, L), where L is any number > 1,
• solve exactly {diff(T(R), R) = T_LS, T(0)=0.5},
• solve numerically {diff(T(R), R) = T_CS, T(h)=U} where U is the value of T at point R=h for the previous solution.

Merge the pieces of solution.
I took arbitrary h=0.25, but this value can be adjusted; a comparison with @Preben Alsholm's solution is also presented
Finally, I've chosen alpha=1, T(0)=1/2 and Phi(0)=0.

 > restart
 > with(plots):
 > Sys := diff(T(R),R)=((1-1/R)*(sqrt(1-(alpha/R)^2*(1-1/R))))^(-1),diff(Phi(R),R)=(alpha/R)^2*(sqrt(1-(alpha/R)^2*(1-1/R)))^(-1);
 (1)
 > T_LS := eval(convert(series(rhs(Sys[1]), R=0, 6), polynom), alpha=1) assuming R > 0; T_RS := eval(convert(series(rhs(Sys[1]), R=1, 6), polynom), alpha=1); h := 0.25: T_CS := eval(rhs(Sys[1]), alpha=1);
 (2)
 > T_Lsol  := Re( dsolve({diff(T(R),R)=T_LS, T(0)=1/2}, T(R)) ): T_C1sol := dsolve({diff(T(R),R)=T_CS, T(h)=eval(rhs(T_Lsol), R=h)}, numeric): T_Rsol  := Re( dsolve({diff(T(R),R)=T_RS, T(1-h)=eval(T(R), T_C1sol(1-h))}, T(R)) ): T_C2sol := dsolve({diff(T(R),R)=T_CS, T(1+h)=eval(rhs(T_Rsol), R=1+h)}, numeric): T_plot := display(   plot(rhs(T_Lsol), R=0...h, color=red, legend=T),   odeplot(T_C1sol, [R, T(R)], R=h..1-h, color=red),   plot(Re(rhs(T_Rsol)), R=1-h..1+h, color=red),   odeplot(T_C2sol, [R, T(R)], R=1+h..2, color=red) ):
 > Phi_LS := eval(convert(series(rhs(Sys[2]), R=0, 6), polynom), alpha=1) assuming R > 0; Phi_RS := eval(rhs(Sys[2]), alpha=1);
 (3)
 > Phi_Lsol := dsolve({diff(Phi(R),R)=Phi_LS, Phi(0)=0}, Phi(R)); Phi_Rsol := dsolve({diff(Phi(R),R)=Phi_RS, Phi(h)=eval(rhs(Phi_Lsol), R=h)}, numeric): Phi_plot := display(   plot(rhs(Phi_Lsol), R=0...h, color=blue, legend=Phi),   odeplot(Phi_Rsol, [R, Phi(R)], R=h..2, color=blue) ):
 (4)
 > display(T_plot, Phi_plot)
 > NumSol := dsolve(eval({Sys}, alpha=1) union {T(1e-10)=1/2, Phi(1e-10)=0}, numeric, events=[[R=1-1e-10, halt]])
 (5)
 > NumSol_plot := odeplot(                  NumSol                  , [[R, T(R)], [R, Phi(R)]]                  , R=1e-10..1-1e-10                  , color=[red, blue]                  , style=[point, point]                  , symbol=[circle, circle]                  , legend=[T, Phi]                ):
 > display(T_plot, Phi_plot, NumSol_plot, view=[default, -10..3])
 >

## Assuming both alpha and k are strictly p...

Assuming both alpha and k are strictly positive:

```Jn := (alpha, k) -> k*Int((BesselK(0, k*Zeta)/Zeta^3)*Zeta*BesselK(1, alpha*Zeta), Zeta = 1 .. infinity):
r  := rand(0. .. 10):
ak := r(), r();
evalf(Jn(ak))
9.631107879, 5.713953505
-8
1.463786783 10
```

Multiply this value by 4*h1*lambda to get the final result

## What is cot?...

You can define cot(x) as 1/tan(x) or cot(x) as tan(Pi/2-x) and even cot(x) as -tan(x+Pi/2) (and while not say that cot(x) = cos(x)/sin(x) ?).
In your case (where does the value 90 come from?) you can obtain evry different simplifications: they all depend on the rules you apply and the order you apply them.

Result y is, IMO, even simpler and easier to read than MMA's result. It uses the subjective simplification tan = 1/cot.

Result z is closer to MMA's but replaces Pi/2 by Phi in order to avoid earlier (and autimatic) simplifications. It uses the subjective simplification cot = 1/tan (which is just one of the many definitions of cot).

The last result starts from z and uses nother subjective simplification cot(Pi/2+alpha/2)=-cot(Pi/2-alpha/2).

```restart:
x := -sin(alpha)*(L + (e + r)/tan(-Pi/2 + alpha/2))/2 + sin(alpha)*L:
eval(%, cot=(u -> 1/tan(u))):
y := collect(%, sin);
/1     1            /1      \\
|- L + - (e + r) tan|- alpha|| sin(alpha)
\2     2            \2      //

x := -sin(alpha)*(L + (e + r)/tan(Phi + alpha/2))/2 + sin(alpha)*L:

eval(%, tan=(u -> 1/cot(u))):
z := collect(%, sin);
/1     1            /      1      \\
|- L - - (e + r) cot|Phi + - alpha|| sin(alpha)
\2     2            \      2      //
eval(z, cot(Phi+alpha/2)=-cot(Phi-alpha/2));
/1     1            /      1      \\
|- L + - (e + r) cot|Phi - - alpha|| sin(alpha)
\2     2            \      2      //
```

simplify.mw

## s := singular(1/sin(x), x): var := (ind...

```s := singular(1/sin(x), x):
var := (indets(%) minus {x})[]:
eval(s, var = `#mo("&#x2124;")`);
```

## As @tomleslie said, one BC is missi...

As @tomleslie said, one BC is missing.

An example with a periodic BC for D[1](u) ( = diff(u(x, t), x) )

```restart:
PDE  := diff(u(x, t), t) = -u(x, t)*diff(u(x, t), x) + 0.1 * diff(u(x, t), x\$2):
IC   := u(x, 0) = sin(x):
BC   := u(0, t)=u(2*Pi, t), D[1](u)(0, t)=D[1](u)(2*Pi, t):
IBC  := {IC, BC}:

pds := pdsolve(PDE, IBC, numeric, time=t, range=0..2*Pi):  # just do as explained in the help page
plots:-display(
seq(
pds:-plot(
t=tau
, numpoints=50
, color=ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12])
, legend=('t'=nprintf("%1.2e", tau))
)
, tau in [seq](0..2, 0.25)
)
)
```

## ExplanationP := plot3d( [ b, ( b^2 ) / 4...

Explanation

`P := plot3d( [ b, ( b^2 ) / 4, -b / 2 ], b = -4 .. 4, c = -4 .. 4, thickness = 5, color = red  ); `

P is made of two elements

`op(P) = MESH(...., COLOR(....)), THICKNESS(...)`

When a MESH structure is plotted the "construction lines" of the mesh are too (default choice).
To avoid plotting them (see PLOT3D help page) you have to use the option STYLE(HIDDEN).

Then (no need to the THICKNESS option)

`PLOT3D( MESH( op(op(P)), STYLE(HIDDEN) ) )`

does the job.

Not that simple but this explains why the curve you got is black: the black boundaries of each element of the MESH "overwrite" the red on their interiors.

 1 2 3 4 5 Page 3 of 5
﻿