## 13579 Reputation

13 years, 102 days

## Probab;y(?)...

can only be done numerically, as in the attached

 > restart;   with(plots):   Pr:=1:   etamax := 10;   eq:= 3/4*F(eta)*diff(theta(eta),eta) = diff(theta(eta),eta,eta),        1/Pr*(1/2*diff(F(eta),eta)^2 - 3/4*F(eta)*diff(F(eta),eta,eta)) = -diff(F(eta),eta,eta,eta) + theta(eta);   bcs := F(0) = 0, theta(0)=1, D(F)(0)=0,theta(etamax)=0,D(F)(etamax)=0;   sol:=dsolve( [ eq, bcs ],                [ F(eta), theta(eta)],                numeric              );   odeplot( sol,            [ [ eta, F(eta)],              [eta, theta(eta) ]            ],            eta=0..etamax,            color=[ red, blue]          );
 >

## A lot better...

See the attached, where I have used arbitrary values for alpha, c, gamma and I am assuming that you are only interested in real roots

 >
 >
 (1)
 > # # Arbitrary values for alpha, c, gamma #   pvals:=[alpha=1, c=2, gamma=3];   poly:=eval( numer(diff(rel, k)), pvals); # # Get all roots of poly including complex ones #   [fsolve(poly, complex)]; # # Select only real roots #   crit_points:=select( type, %, realcons);
 (2)
 >
 (3)

## The basics...

If you examine the numerator of diff(rel, k), which has to be zero to get the critical points, you will find that it is 7-th order polynomial in k. Only polynomials of order 4 or less can be solved analytically - so your solve() commands will produce nothing useful.

If you have numerical values for {alpha, c, gamma, w}, tyhen you can use fsolve()

## I doubt...

that

most of papers or books give directly the formula without any proof.

because as https://en.wikipedia.org/wiki/Characteristic_equation_(calculus) states

In mathematics, the characteristic equation (or auxiliary equation[1]) is an algebraic equation of degree n upon which depends the solution of a given nth-order differential equation[2] or difference equation.[3][4] The characteristic equation can only be formed when the differential or difference equation is linear and homogeneous, and has constant coefficients.

Note that the characteristic equation can only be found when the differential equation is linear and you specifically state

I have a system of nonlinear differential equation

## Use the inert form of '/'...

ie '%/', as in the attached

 > restart; # # Evaluates the argument before generating latex #   latex((4*n-1)/9-7/16*n); # # Use the inert form of '/' #   latex((4*n-1)%/9-7%/16*n);
 \frac{n}{144}-\frac{1}{9} \frac{4 n -1}{9}-\frac{7}{16} n
 >

## Hmmm...

The " code" which you present has obviously been cut-and-pasted from a workseet using 2D input. It has multiple instances where the 2D input is interpreting spaces as implicit multiplication. Consider the snips from your (incorrect) code, with the highlighted multiplication signs below

display*([draw*[A(color = black, symbol = solidcircle, symbolsize = 6),
B(color = black, symbol = solidcircle, symbolsize = 6),
C(color = black, symbol = solidcircle, symbolsize = 6),
ABC(color = blue)],
textplot*([[coordinates(A)[], "A"],
[coordinates(B)[], "B"],
[coordinates(C)[], "C3"]],
align = [above, right])],
axes = none,
title = "Lemme du Chevron");

There are also multiple unbalanced parentheses.

I have fixed all of these in the attached.

 > restart; with(geometry): with(plots): _EnvHorizontalName = 'x': _EnvVerticalName = 'y': EQ := proc(M, N) local eq, sol; eq := simplify(expand((y - M[2])/(x - M[1]) - (N[2] - M[2])/(N[1] - M[1]))); sol := solve(eq, y); RETURN(y = sol); end proc: _local(D): point(A, [-2, 7]): point(B, [-5, -2]): point(C, [8, -7]): point(E, [1, 4]): EQ([-5, -2], [8, -7]): point(D, [1, subs(x = 1, rhs(%))]): dsegment(sgAD, [A, D]): BD := distance(B, D): DC := distance(C, D): triangle(ABC, [A, B, C]): area(ABC): triangle(ABD, [A, B, D]): area(ABD): triangle(ADC, [A, D, C]): area(ADC): is(area(ABD)/area(ADC) = BD/DC): triangle(EBD, [E, B, D]): area(EBD): triangle(EDC, [E, D, C]): area(EDC): triangle(AEC, [A, E, C]): area(AEC): triangle(ABE, [A, B, E]): area(ABE): is(area(ABE)/area(AEC) = BD/DC): display ( [ draw     ( [ A(color = black, symbol = solidcircle, symbolsize = 6),         B(color = black, symbol = solidcircle, symbolsize = 6),         C(color = black, symbol = solidcircle, symbolsize = 6),         ABC(color = blue)       ]     ),     textplot     ( [ [coordinates(A)[], "A"],         [coordinates(B)[], "B"],         [coordinates(C)[], "C3"]       ],       align = [above, right]     )   ],   axes = none,   title = "Lemme du Chevron" );
 >

makes a fairly trivial calculation look really complicated. In the attached, I deleted absolutely everything that you do not need. So much simpler to understand what is happening!

 > restart:   with(plots):   Params := [ fw = 0.2, M = 0.5, Q = 0.5, Pr = 6.2, phi = 0.05,               rf = 997.1, kf = 0.613, cpf = 4179, btf = 0.00003,               p1 = 0.01, p2 = 0.05, p3 = 0.05, rs1 = 5100, ks1 = 3007.4,               cps1 = 410, bs1 = 0.0002, rs2 = 2200, ks2 = 5000, cps2 = 790,               bs2 = 0.0005, rs3 = 3970, ks3 = 40, cps3 = 765, bs3 = 0.0004,               A1 = B1*p1 + B2*p2 + B3*p3, B1 = 1 + 2.5*phi + 6.2*phi^2,               B2 = 1 + 13.5*phi + 904.4*phi^2, B3 = 1 + 37.1*phi + 612.6*phi^2,               B4 = (ks1 + 2*kf - 2*phi*(kf - ks1))/(ks1 + 2*kf + phi*(kf - ks1)),               B5 = (ks2 + 3.9*kf - 3.9*phi*(kf - ks2))/(ks2 + 3.9*kf + phi*(kf - ks2)),               B6 = (ks3 + 4.7*kf - 4.7*phi*(kf - ks3))/(ks3 + 4.7*kf + phi*(kf - ks3)),               A2 = 1 - p1 - p2 - p3 + p1*rs1/rf + p2*rs2/rf + p3*rs3/rf,               A3 = B4*p1 + B5*p2 + B6*p3, A4 = 1 - p1 - p2 - p3 + p1*rs1*cps1/(rf*cpf) + p2*rs2*cps2/(rf*cpf) + p3*rs3*cps3/(rf*cpf)           ]:   ODEs:= [ A1*(diff(f(x), x, x, x))/(A2*phi)-(diff(f(x), x))^2-M^2*(f(x))+f(x)*(diff(f(x), x, x)),            A4*Pr*phi*(diff(Theta(x), x, x))/A3+f(x)*(diff(Theta(x), x))+Q*Theta(x)          ]:   ODEs:=eval[recurse](ODEs, Params):   BCs := [f(0) = fw, D(f)(0) = 1, Theta(0) = 1, D(f)(1) = 0, Theta(1) = 0]:   BCs := eval[recurse](BCs,Params):   sol:=dsolve([ODEs[], BCs[]], numeric);   odeplot( sol,            [ [ x, f(x)],              [x, Theta(x)]            ],            x=0..1,            color=[red, blue]          );
 >

## Check the help...

at ?plot3d,coords. Note in particular the interpretation of the expression being plotted, which is stated as

For spherical coordinates the interpretation is:

"plot3d(r(theta, phi), theta = a .. b, phi = c .. d, coords =spherical)"

So if you want to plot a sphere of radius 1, you would use something like the attached

(By the way, there are many other ways to plot a sphere!!)

 > restart;   plot3d( 1,           theta = 0..Pi,           phi = 0..2*Pi,           coords = spherical,           style =surface,           scaling=constrained         );
 >

## Nothing wrong...

with generating a conic from five points, although I would probably let the the command conic() in Maple's geometry() package do most of the "work" - as in the attached.

 > restart:   with(geometry):      Pn := [ (5*x + 12*y)*(-384 + 48*x + 55*y)/(225*x^2 + 225*y^2 - 1800*x - 602*y),           2*(12*x - 5*y)*(-384 + 48*x + 55*y)/(225*x^2 + 225*y^2 - 1800*x - 602*y)         ]:   line(l1, 2+3*x-y, [x,y]): # # Define the conic from five points #   conic( c,          [ seq            ( point              ( cat(P__, j),                eval                ( eval( Pn, isolate( Equation(l1), y ) ),                  x=j                )              ),              j=1..5            )          ],          [x,y]         ): # # Some properties of the conic #   form(c);   simplify(coordinates(center(c)));   simplify~(coordinates~(foci(c)));   simplify(Equation(c)); # # Draw the conic #   c_cen:=coordinates(center(c)):   vr:= z-> [ c_cen[1]-z..c_cen[1]+z,              c_cen[2]-z..c_cen[2]+z            ]:   draw( c(color=red),         view=vr(30)       );
 >

## This would be so much simpler to diagnos...

if you ha uploaded a worksheet using the big green up-arrow in the Mapleprimes toolbar,

Taking a wild guess, try replacing

YFtpt5 := TtT -> solution(TtT, 0.5)[3]

with

YFtpt5 := TtT -> rhs(solution(TtT, 0.5)[3])

## Why...

extract data from a plot, when it is simpler to generate the data used by the plot?

See the final execution group in the attached where I have arrabged the the data matrix is in the same "orientation" as the contour plot, just to make comparing the two a bit simpler.

data.mw

## Hmmm...

As well as fixing syntax problems, I had to make many wild guesses about what you are actually trying to achieve.

Maybe(??!) the attached will be helpful

 > restart:   N := 4: M := .1: Kp := .1: Gr := 0.1e-1: Gc := 0.1e-1: Pr := 1: S := 0.1e-1: Sc := .78: Kc := 0.1e-1: La := 1:   f__N:=  sum((p^(i))*f[i](x), i = 0 .. N) :   Theta__N:=  sum((p^(i))*Theta[i](x), i = 0 .. N) :   Phi__N:= sum((p^(i))*Phi [i] (x), i = 0 .. N):   HPMEq1 := (1-p)*(diff(f__N, x, x, x))+p*(diff(f__N, x, x, x)+(1/2)*(diff(f__N, x, x))*f__N-(M^2+Kp)*(diff(f__N, x)-La)+Gr*Theta__N+Gc*Phi__N):   HPMEq2 := (1-p)*(diff(Theta__N, x, x))/Pr+p*((diff(Theta__N, x, x))/Pr+(1/2)*(diff(Theta__N, x))*f__N+S*Theta__N):   HPMEq3 := (1-p)*(diff(Phi__N, x, x))/Sc+p*((diff(Phi__N, x, x))/Sc+(1/2)*(diff(Phi__N, x))*f__N+Kc*Phi__N):   for i from 0 to N do       equ[1][i] := coeff(HPMEq1, p, i) = 0:   end do:   for i from 0 to N do       equ[2][i] := coeff(HPMEq2, p, i) = 0:   end do:   for i from 0 to N do       equ[3][i] := coeff(HPMEq3, p, i) = 0:   end do:   cond[1][0] := f[0](0) = 0, (D(f[0]))(0) = 0, Theta[0](0) = 1, Phi[0](0) = 1, Theta[0](5) = 0, Phi[0](5) = 0, (D(f[0]))(5) = 1:   for j to N do        cond[1][j] := f[j](0) = 0, (D(f[j]))(0) = 0, Theta[j](0) = 0, Phi[j](0) = 0, Theta[j](5) = 0, Phi[j](5) = 0, (D(f[j]))(5) = 0:   end do:   ans:=dsolve({equ[1][0], equ[2][0], equ[3][0], cond[1][0]}):   for k from 1 by 1 to 4 do       ans:=union( ans, dsolve( { eval({equ[1][k], equ[2][k], equ[3][k]}, ans)[], cond[1][k]})):   od:   Phi:= evalf(eval(Phi[k-1](x), ans));   Theta:= evalf(eval(Theta[k-1](x), ans));   f:= evalf(eval(f[k-1](x), ans));   plot([Phi, Theta, f], x=0..5, color=[red, blue, green]);
 >

## Fairly obvious problem...

Consider the command

x:=x+1

If 'x' has previously been assigned a numeric value, the this will work as expected. But if 'x' has not been assigned a numeric value, then the above will evaluate (recursively) to

(x+1)+1;
((x+1)+1)+1
(((x+1)+1)+1)+1....

and so on. After a sufficient number of such recursive evaluations, Maple will halt and generate an error.

Almost any assignment can produce this. Unless you post a worksheet, using the big green up-arrow in the Mapleprimes toolbar, no-one will be able to tell "which command created it" - although it is probaly related to an assignment to the quantity 'PD'

## Problem...

with RealDomain() - and a version issue.

Maple 2022.2, returns x=-4, x=2, as solutions for your equation, whether or not you invoke RealDomain(). When RealDomain() is not specified, both sides of your equation evaluate to the same (complex) number. However, if one attempts to evaluate your equation at x=2, after invoking RealDomain(), then Maple returns undefined - see the attached.

 > # # No real domain assumption #   restart;   eqn:= log[2](x^2 - 6*x) = 3 + log[2](1 - x);   ans:=solve(eqn, x); # # x=2 gives the function complex values! #   simplify(eval(eqn, x=ans[1]));   simplify(eval(eqn, x=ans[2]));
 (1)
 > # # Use real domain #   with(RealDomain):   ans:=solve(eqn, x);   simplify(eval(eqn, x=ans[1]));   simplify(eval(eqn, x=ans[2]));
 (2)

## Clarify...

Apart from the restricted viewing ranges in the required plots (which I have removed in the attached), I don't see any problem with this worksheet

 >
 >
 >
 (1)
 >
 (2)
 >
 (3)
 >
 >
 >
 >
 (4)
 >
 >
 (5)
 >
 (6)
 >
 (7)
 >
 (8)
 >
 (9)
 >
 (10)
 >
 >
 >
 >
 >
 >
 >
 >