Advanced Search

tomleslie

13836 Reputation

20 Badges

14 years, 342 days
Contact tomleslie

MaplePrimes Activity


These are replies submitted by tomleslie

Just assign...

tomleslie 13836
October 17 2022
0 0

@romanrieme 

the result to something as in the last execution group of the attached, where I have used Acer's worksheet, because it is a better solution than mine

> 

maxPt := Optimization:-Maximize(180*arccos(temp(t))/Pi, t = 0 .. 2); maxPt := [rhs(maxPt[2][]), maxPt[1]]; 2*maxPt[2]

[HFloat(22.00633969122702), [t = HFloat(0.6180259177597804)]]

[HFloat(0.6180259177597804), HFloat(22.00633969122702)]

HFloat(44.01267938245404)

Download maxPt.mw

Actually...

tomleslie 13836
October 16 2022
0 0

@Carl Love 

I think it may be a "units" issue - the attached, which is OP's original (with all references to units removed) appears to work much better for both plot() with option axis[1] = [mode = log] and semilogplot()

Impedance Graph

 

restart

Example

A noise voltage specification of
 `ΔV__noise` := 50*10^(-3)

Worst case transient current
 I__wc_transient := 1

Z__target := `ΔV__noise`/I__wc_transient

1/20

 (1)

For a clock frequency
f__clk := 2*10^9 

the clock period  `τ__clk` := 1/f__clk

1/2000000000

 (2)

NULL

The maximum number of switching gate cycles
 N__sw_gates := 20

The transient time

t__transient := N__sw_gates*`τ__clk`

1/100000000

 (3)

 

The roll-off frequency for the maximum transient current         f__roll_off := .35/t__transient

35000000.00

 (4)

 Therefore the target impedance for this example is

Z__target_ac := proc (f) local X1, Z__target_1; Z__target_1 := Z__target; X1 := proc (f) options operator, arrow; evalf(20*log[10](Z__target_1))+20*log[10](sqrt(1+f^2/f__roll_off^2)) end proc; 10^((1/20)*X1(f)) end proc

NULL

Z__target_ac(3.82*10^6)

0.5029692240e-1

 (5)

``

f__range := 10^6 .. 10^10

NULL

Target Impedance plot

NULL

plot(Z__target_ac, f__range, gridlines, axis[1] = [mode = log], axes = boxed, title = "Target Impedance", color = red, labels = ["Frequency (Hz)", "Impedance (Ω)"]); plots:-semilogplot(Z__target_ac, f__range, gridlines, axes = boxed, title = "Target Impedance", color = red, labels = ["Frequency (Hz)", "Impedance (Ω)"])

  

This graph should have a horizontal line all the way along to the roll-off/corner frequency and then rise.

 

NULL


 

Download noUnits.mw

A disappointing thing...

tomleslie 13836
October 16 2022
0 0

@Carl Love 

is that, while I agree with your reasoning, the command semilogplot() also uses "linear"  selection for the values of the independent variable - so has the same problem as plot() with the option axis=[mode=log].

I don't know why but I would have expected semlogplot() to do better? It didn't, which is why I came up wiyj the multiple plot approach

I have just run a chich check...

tomleslie 13836
October 14 2022
0 0

on mmcdara's response, and the failure problem is a Maple version issue.

Maple version         Outcome

        18                     works

       2015                  works

       2016                  works

       2017                  works

       2018                  fails with "Error"

       2019                  fails with "Error"

       2020                  fails with "Error"

       2021                  fails with "Error"

       2022                  fails with "Error"

One of the reasons...

tomleslie 13836
October 14 2022
0 0

@mmcdara 

I posted a" pale copy" was that when I tried to run your original worksheet, it produced an "Error" quite early which affected all subsequent calculations, so didn't even produce an answer for the OP's original case of sa=2, sd=1. See the attached.

I am assuming that this is a Maple version issue.  (I'm using Maple 2022.) However I didn't have the time/patience to work through the last few versions of Maple to discover the one in which your worksheet actually "functioned".

> 

restart:
interface(version);

`Standard Worksheet Interface, Maple 2022.1, Windows 7, May 26 2022 Build ID 1619613`

 (1)
> 

f := (x, x__0, S, C) -> -C*(exp(S*(x - x__0)^2*(-1/2)) - 1)*Heaviside(x - x__0)

proc (x, x__0, S, C) options operator, arrow; -C*(exp(-(1/2)*S*(x-x__0)^2)-1)*Heaviside(x-x__0) end proc

 (2)
> 

X0   := 0:
fDOG := unapply(f(x, X0, S__a, 1) - f(x, X0, S__d, 1), (x, S__a, S__d))

proc (x, S__a, S__d) options operator, arrow; -(exp(-(1/2)*S__a*x^2)-1)*Heaviside(x)+(exp(-(1/2)*S__d*x^2)-1)*Heaviside(x) end proc

 (3)
> 

DfDOG := diff(fDOG(x, S__a, S__d), x) assuming x > X0

S__a*x*exp(-(1/2)*S__a*x^2)-S__d*x*exp(-(1/2)*S__d*x^2)

 (4)
> 

solve({DfDOG=0, x > X0}, x) assuming S__a > S__d;
tpeak := eval(x, %[1])

piecewise(0 < sqrt((S__a-S__d)*ln(S__a/S__d))/(S__a-S__d), [{x = sqrt(2)*sqrt(ln(S__a)-ln(S__d))/sqrt(S__a-S__d)}], sqrt((S__a-S__d)*ln(S__a/S__d))/(S__a-S__d) < 0, [{x = -sqrt(2)*sqrt(ln(S__a)-ln(S__d))/sqrt(S__a-S__d)}], [])

  

Error, invalid input: eval received piecewise(0 < 1/(S__a-S__d)*((S__a-S__d)*ln(S__a/S__d))^(1/2),[{x = 2^(1/2)/(S__a-S__d)^(1/2)*(ln(S__a)-ln(S__d))^(1/2)}],1/(S__a-S__d)*((S__a-S__d)*ln(S__a/S__d))^(1/2) < 0,[{x = -2^(1/2)/(S__a-S__d)^(1/2)*(ln(S__a)-ln(S__d))^(1/2)}],[])[1], which is not valid for its 2nd argument, eqns

  
> 

hpeak := fDOG(tpeak, S__a, S__d) assuming S__a > S__d, S__d > 0

-(exp(-(1/2)*S__a*tpeak^2)-1)*Heaviside(tpeak)+(exp(-(1/2)*S__d*tpeak^2)-1)*Heaviside(tpeak)

 (5)
> 

half_hpeak := hpeak/2;

-(1/2)*(exp(-(1/2)*S__a*tpeak^2)-1)*Heaviside(tpeak)+(1/2)*(exp(-(1/2)*S__d*tpeak^2)-1)*Heaviside(tpeak)

 (6)
> 

eq := fDOG(x, S__a, S__d) = half_hpeak assuming x > X0;

MyChange := ln(S__a/S__d) = (u/sqrt(2))^2*(S__a-S__d); # u is > 0
eval(eq, MyChange);
eval(%, [S__a=2, S__d=1]);
 

-exp(-(1/2)*S__a*x^2)+exp(-(1/2)*S__d*x^2) = -(1/2)*(exp(-(1/2)*S__a*tpeak^2)-1)*Heaviside(tpeak)+(1/2)*(exp(-(1/2)*S__d*tpeak^2)-1)*Heaviside(tpeak)

  

ln(S__a/S__d) = (1/2)*u^2*(S__a-S__d)

  

-exp(-(1/2)*S__a*x^2)+exp(-(1/2)*S__d*x^2) = -(1/2)*(exp(-(1/2)*S__a*tpeak^2)-1)*Heaviside(tpeak)+(1/2)*(exp(-(1/2)*S__d*tpeak^2)-1)*Heaviside(tpeak)

  

-exp(-x^2)+exp(-(1/2)*x^2) = -(1/2)*(exp(-tpeak^2)-1)*Heaviside(tpeak)+(1/2)*(exp(-(1/2)*tpeak^2)-1)*Heaviside(tpeak)

 (7)
> 

xofu := [solve(%, x)];
 

[(-2*ln(1/2+(1/2)*(1+2*Heaviside(tpeak)*(exp(-(1/2)*tpeak^2))^2-2*Heaviside(tpeak)*exp(-(1/2)*tpeak^2))^(1/2)))^(1/2), -(-2*ln(1/2+(1/2)*(1+2*Heaviside(tpeak)*(exp(-(1/2)*tpeak^2))^2-2*Heaviside(tpeak)*exp(-(1/2)*tpeak^2))^(1/2)))^(1/2), (-2*ln(1/2-(1/2)*(1+2*Heaviside(tpeak)*(exp(-(1/2)*tpeak^2))^2-2*Heaviside(tpeak)*exp(-(1/2)*tpeak^2))^(1/2)))^(1/2), -(-2*ln(1/2-(1/2)*(1+2*Heaviside(tpeak)*(exp(-(1/2)*tpeak^2))^2-2*Heaviside(tpeak)*exp(-(1/2)*tpeak^2))^(1/2)))^(1/2)]

 (8)
> 

solve(MyChange, u);
x_half_peak := eval(xofu, u=%[1]);

2^(1/2)*((S__a-S__d)*ln(S__a/S__d))^(1/2)/(S__a-S__d), -2^(1/2)*((S__a-S__d)*ln(S__a/S__d))^(1/2)/(S__a-S__d)

  

[(-2*ln(1/2+(1/2)*(1+2*Heaviside(tpeak)*(exp(-(1/2)*tpeak^2))^2-2*Heaviside(tpeak)*exp(-(1/2)*tpeak^2))^(1/2)))^(1/2), -(-2*ln(1/2+(1/2)*(1+2*Heaviside(tpeak)*(exp(-(1/2)*tpeak^2))^2-2*Heaviside(tpeak)*exp(-(1/2)*tpeak^2))^(1/2)))^(1/2), (-2*ln(1/2-(1/2)*(1+2*Heaviside(tpeak)*(exp(-(1/2)*tpeak^2))^2-2*Heaviside(tpeak)*exp(-(1/2)*tpeak^2))^(1/2)))^(1/2), -(-2*ln(1/2-(1/2)*(1+2*Heaviside(tpeak)*(exp(-(1/2)*tpeak^2))^2-2*Heaviside(tpeak)*exp(-(1/2)*tpeak^2))^(1/2)))^(1/2)]

 (9)
> 

{abs~(x_half_peak)[]}

{2^(1/2)*abs(ln(1/2-(1/2)*(1+2*Heaviside(tpeak)*(exp(-(1/2)*tpeak^2))^2-2*Heaviside(tpeak)*exp(-(1/2)*tpeak^2))^(1/2)))^(1/2), 2^(1/2)*abs(ln(1/2+(1/2)*(1+2*Heaviside(tpeak)*(exp(-(1/2)*tpeak^2))^2-2*Heaviside(tpeak)*exp(-(1/2)*tpeak^2))^(1/2)))^(1/2)}

 (10)
> 

sa := 2:
sd := 1:

eval(x_half_peak, [S__a=sa, S__d=sd]):
select(is@`>`, %, 0);
evalf(%);

[]

  

[]

 (11)
> 

# Applying "select" directly on x_half_peak doesn't work.
# I believe the following is correct but this must be checked

`correct?` := {abs~(x_half_peak)[]};
eval(`correct?`, [S__a=sa, S__d=sd]):
evalf(%);

{2^(1/2)*abs(ln(1/2-(1/2)*(1+2*Heaviside(tpeak)*(exp(-(1/2)*tpeak^2))^2-2*Heaviside(tpeak)*exp(-(1/2)*tpeak^2))^(1/2)))^(1/2), 2^(1/2)*abs(ln(1/2+(1/2)*(1+2*Heaviside(tpeak)*(exp(-(1/2)*tpeak^2))^2-2*Heaviside(tpeak)*exp(-(1/2)*tpeak^2))^(1/2)))^(1/2)}

  

{1.414213562*abs(ln(.5000000000-.5000000000*(1.+2.*Heaviside(tpeak)*(exp(-.5000000000*tpeak^2))^2-2.*Heaviside(tpeak)*exp(-.5000000000*tpeak^2))^(1/2)))^(1/2), 1.414213562*abs(ln(.5000000000+.5000000000*(1.+2.*Heaviside(tpeak)*(exp(-.5000000000*tpeak^2))^2-2.*Heaviside(tpeak)*exp(-.5000000000*tpeak^2))^(1/2)))^(1/2)}

 (12)
> 

 

Download oddError.mw

 

I doubt it...

tomleslie 13836
October 13 2022
0 0

Some of the expressions get really ugly. About the best I thinnk I can do is to put the (slightlly modified) previou calculation in a procedure wrapper, so that you can call it with any (numeric) values for Sa, Sb, as shown in the attached

> 

  restart;
  with(plots):

  doSums:= proc(Sa, Sb)
                local
                    #
                    # OP's definition of a "gaussian"
                    #
                    f:=(x, x0, S, Gmax)-> -Gmax*(exp(S*(x - x0)^2*(-1/2)) - 1),
                    sv:=sort([Sa,Sb]),
                    ans, pt1, pt2, pt3, sols, p1, fwhm;
                #
                # Get the coordinates of the maximum of
                # the difference between two gaussians
                # Select the one for x>0, Store these
                # coordinates as pt1.
                #
                ans:= maximize
                      ( f(x, 0, sv[2], 1)-f(x, 0, sv[1], 1),
                        location=true
                      )[2]:
                ans:=select( j-> `if`(evalf(rhs(j[1][]))>0,true,false), ans)[];
                pt1:=[ rhs(ans[1][]), ans[2]];
                #
                # Find the x-values of the gaussian difference
                # function, where the y-values are half the
                # maximum obtained above.
                #
                # Store the relevant points as pt2 and pt3
                #
                sols:= [ solve
                         ( [ x>0,
                             f(x, 0, sv[2], 1)-f(x, 0, sv[1], 1 )=pt1[2]/2
                           ],
                           explicit
                         )
                       ]:
                pt2:= [ rhs(sols[1][]), pt1[2]/2];
                pt3:= [ rhs(sols[2][]), pt1[2]/2];
                #
                # Display the two gaussians, the difference between
                # them and the three points of interest
                #
                p1:= display([ plot
                               ( [ f(x, 0, sv[2], 1),
                                   f(x, 0, sv[1], 1),
                                   f(x, 0, sv[2], 1)-f(x, 0, sv[1], 1)
                                 ],
                                 x=-10..10,
                                 color=[red, blue, green]
                               ),
                               plot
                               ( [ pt1,
                                   pt2,
                                   pt3
                                 ],
                                 symbol=solidcircle,
                                 symbolsize=12,
                                 color=black,
                                 style=point
                               )
                             ]);
                #
                # Full width at half maximum of the gaussian difference
                # function is the difference between the x-coordinates of
                # pt2 and pt3.
                #
                fwhm:= pt2[1]-pt3[1];
                return [pt1, pt2, pt3, p1, fwhm ];
          end proc:
> 

 ans:=doSums(1,2):
 ans[1..3];
 ans[4];

[[2^(1/2)*ln(2)^(1/2), 1/4], [(-2*ln(1/2-(1/4)*2^(1/2)))^(1/2), 1/8], [(-2*ln(1/2+(1/4)*2^(1/2)))^(1/2), 1/8]]

  

 
> 

 ans:=doSums(1,10):
 ans[1..3];
 ans[4];

[[(1/3)*2^(1/2)*(ln(2)+ln(5))^(1/2), -exp(-(10/9)*ln(2)-(10/9)*ln(5))+exp(-(1/9)*ln(2)-(1/9)*ln(5))], [(-2*ln(RootOf(200*_Z^10+9*10^(8/9)-200*_Z, index = 1)))^(1/2), -(1/2)*exp(-(10/9)*ln(2)-(10/9)*ln(5))+(1/2)*exp(-(1/9)*ln(2)-(1/9)*ln(5))], [(-2*ln(RootOf(200*_Z^10+9*10^(8/9)-200*_Z, index = 2)))^(1/2), -(1/2)*exp(-(10/9)*ln(2)-(10/9)*ln(5))+(1/2)*exp(-(1/9)*ln(2)-(1/9)*ln(5))]]

  

 
> 

 

Download gmaxProc.mw

I suppose...

tomleslie 13836
October 13 2022
0 0

it depends on which 8 points you want to throw way!

Seriously, if you want help, get into the habit of posting an executable worksheet using the big green up-arrow in the Mapleprimes toolbar - because I suspect that the answer to your problem would be to make all the matrices the same size in the first place - but since I have no idea how you produced the relevant graph(s) in the first place, I really can't help much

Hmmmm...

tomleslie 13836
October 11 2022
0 0

Well since you are now asking for the quantity -log(H), I am going to assume that you are only interested in solutions for 'H' which are positive -otherwise log(H) would be complex. THis is pretty trivial to arrange.

\You ask for tw0 possible sequences of the independent variable V__b, either

  1. 5 to 15 in steps of 0.2, or
  2. 0-5 x 1.0, 5-9 x 0.5, 9-11 x 0.1, 11-15 x 0.5, 15-20 x 1.0

again pretty trivial to produce (both options are shown in the attached)

Plots of -log(H) versus V__b: plots for both of the above sets of V__b values are shown in the attached. As you woulld expect, these differ in minor details

Export to Excel - (why?), Standard export to Excel is also shown in the attached for both results matrizes, although you will have to change the file path to something appropriate for your machine - I just used my own (Windows) desktop for test purposes.

(The graphs in the attached "render" much better in the Maple worksheet than they do on this site - take my word for it!)

> 

  restart;
  with(plots):
  interface(rtablesize=60):
#
# Define desired expression
#
  expr1:= V__b=V__a*(( C__a/(1+H/K__a)-H+K__w/H)/(C__b/(1+K__w/H/K__b)+H-K__w/H));
#
# Define parameter values
#
  params:=[C__a=0.1, C__b=0.1, K__a=1.74*10^(-5), K__b=10^20, K__w=10^(-14), V__a=10]:
#
# Define a function which will return the value of 'H'
# for the set of parameters above, and a supplied value for
# 'V__b'
#
  getSol:= z-> -log10(rhs(solve([eval(expr1,[params[], V__b=z]), H>0])[])):
#
# Apply the above function for a aequence of values for V__b
#
# Results at OP's request are for V__b=5..15 in steps of 0.2
#  
  M1:= Matrix( [ [ 'V__b', 'H1'],
                   seq( [j, getSol(j)], j=5..15, 0.2)
               ]
             );
  pointplot( M1[2..-1,1..2],
             style=pointline,
             labels=[V__b, "-log10(H)"]
           );
  
 

expr1 := V__b = V__a*(C__a/(1+H/K__a)-H+K__w/H)/(C__b/(1+K__w/(H*K__b))+H-K__w/H)

  

Matrix(%id = 36893488148087723300)

  

 
> 

#
# Results at OP's request are for
# V__b=0-5 x 1.0, 5-9 x 0.5, 9-11 x 0.1, 11-15 x 0.5, 15-20 x 1.0
#
  VbVals:=[ seq(j, j=0..5,1),
            seq(j,j=5.5..9,0.5),
            seq(j, j=9..11, 0.1),
            seq(j, j=11.5..15, 0.5),
            seq( j, j=16..20, 1)
          ]:
  M2:= Matrix( [ [ 'V__b', 'H1'],
                   seq( [j, getSol(j)], j in VbVals)
               ]
             );
  pointplot( M2[2..-1,1..2],
             style=pointline,
             labels=[V__b, "-log10(H)"]
           );

Matrix(49, 2, {(1, 1) = `#msub(mi("V"),mi("b"))`, (1, 2) = H1, (2, 1) = 0, (2, 2) = 2.882589723, (3, 1) = 1, (3, 2) = 3.813306454, (4, 1) = 2, (4, 2) = 4.159642546, (5, 1) = 3, (5, 2) = 4.392562210, (6, 1) = 4, (6, 2) = 4.584019592, (7, 1) = 5, (7, 2) = 4.759903668, (8, 1) = 5.5, (8, 2) = 4.846987785, (9, 1) = 6.0, (9, 2) = 4.935877607, (10, 1) = 6.5, (10, 2) = 5.028590975, (11, 1) = 7.0, (11, 2) = 5.127689535, (12, 1) = 7.5, (12, 2) = 5.236806939, (13, 1) = 8.0, (13, 2) = 5.361723091, (14, 1) = 8.5, (14, 2) = 5.512971651, (15, 1) = 9.0, (15, 2) = 5.713870000, (16, 1) = 9, (16, 2) = 5.713870000, (17, 1) = 9.1, (17, 2) = 5.764423299, (18, 1) = 9.2, (18, 2) = 5.820319223, (19, 1) = 9.3, (19, 2) = 5.883003266, (20, 1) = 9.4, (20, 2) = 5.954591887, (21, 1) = 9.5, (21, 2) = 6.038365647, (22, 1) = 9.6, (22, 2) = 6.139819614, (23, 1) = 9.7, (23, 2) = 6.269253954, (24, 1) = 9.8, (24, 2) = 6.449790237, (25, 1) = 9.9, (25, 2) = 6.755189673, (26, 1) = 10.0, (26, 2) = 8.729262664, (27, 1) = 10.1, (27, 2) = 10.69685412, (28, 1) = 10.2, (28, 2) = 10.99569123, (29, 1) = 10.3, (29, 2) = 11.16963085, (30, 1) = 10.4, (30, 2) = 11.29243301, (31, 1) = 10.5, (31, 2) = 11.38721819, (32, 1) = 10.6, (32, 2) = 11.46428546, (33, 1) = 10.7, (33, 2) = 11.52912875, (34, 1) = 10.8, (34, 2) = 11.58502746, (35, 1) = 10.9, (35, 2) = 11.63409687, (36, 1) = 11.0, (36, 2) = 11.67778123, (37, 1) = 11.5, (37, 2) = 11.84365304, (38, 1) = 12.0, (38, 2) = 11.95860745, (39, 1) = 12.5, (39, 2) = 12.04575758, (40, 1) = 13.0, (40, 2) = 12.11539348, (41, 1) = 13.5, (41, 2) = 12.17300023, (42, 1) = 14.0, (42, 2) = 12.22184879, (43, 1) = 14.5, (43, 2) = 12.26404646, (44, 1) = 15.0, (44, 2) = 12.30103002, (45, 1) = 16, (45, 2) = 12.36317792, (46, 1) = 17, (46, 2) = 12.41373429, (47, 1) = 18, (47, 2) = 12.45593197, (48, 1) = 19, (48, 2) = 12.49184452, (49, 1) = 20, (49, 2) = 12.52287875})

  

 
> 

#
# Export both matrixes to Excel (why?)
#
# NB. OP will have to change the file path to something
# appropriate for his/her machine
#
  ExcelTools:-Export( M1, "C:/Users/TomLeslie/Desktop/M1.xlsx", 1, "B2");
  ExcelTools:-Export( M2, "C:/Users/TomLeslie/Desktop/M2.xlsx", 1, "B2")
> 

 

Download doSum3.mw

 

You should be aware...

tomleslie 13836
October 11 2022
0 0

that the equation you want to solve produces four possible solutions for each value of V__b - and (obviously) I have no idea which one you want, although he fact that three of these are negative may be a clue!

The attached is also using V__b=1..20, again because I have no idea which values you actually want

> 

  restart;
  interface(rtablesize=30):
#
# Define desired expression
#
  expr1:= V__b=V__a*(( C__a/(1+H/K__a)-H+K__w/H)/(C__b/(1+K__w/H/K__b)+H-K__w/H));
#
# Define parameter values
#
  params:=[C__a=0.1, C__b=0.1, K__a=1.74*10^(-5), K__b=10^20, K__w=10^(-14), V__a=10]:
#
# Define a function which will return the value of 'H'
# for the set of parameters above, and a supplied value for
# 'V__b'
#
  getSol:= z-> solve(eval(expr1,[params[], V__b=z]), H):
#
# Apply the above function for a aequence of values for V__b
#
  getSol:= z-> solve(eval(expr1,[params[], V__b=z]), H):
  Matrix( [ [ 'V__b', 'H1', 'H2', 'H3', 'H4'],
            seq( [j, getSol(j)], j=1..20)
          ]
        );  
 

[10, 10]

  

expr1 := V__b = V__a*(C__a/(1+H/K__a)-H+K__w/H)/(C__b/(1+K__w/(H*K__b))+H-K__w/H)

  

Matrix(%id = 36893488148162589628)

 (1)
> 

 

Download doSum2.mw

I have no idea...

tomleslie 13836
October 10 2022
0 0

@delvin 

what you are trying to achieve. The attached fixes the problem with the diff() command and generates a lot more solutions. Are any of these of interest?

Download solEq3.mw

Of course...

tomleslie 13836
October 06 2022
0 0

@Ahmed111 

you can - shown in the revised worksheet attached (The plot looks better in a Maple worksheet than it renders on this site  - honest!)

> 

  restart;
> 

#
# Make the asssumptions c>0, g>0, and suppress
# the assumption-indicator in Maple output
#
  assume(c>0, g>0);
  interface(showassumed=0):
> 

  q := diff(w(x), x)^2/2 + w(x)^4/8 - c*w(x)^2/2 - g = 0;
#
# If one assumes for them moment that w(x)=constant
# then the differential term in q above will be
# identically zero,and one essentially has a quartic
# equation in w(x). In general suxh a quartic equation
# will have four roots, and one can solve for produce
# these roots using a simple solve() command
#
  solve(eval( q, diff(w(x),x)=0), w(x));

(1/2)*(diff(w(x), x))^2+(1/8)*w(x)^4-(1/2)*c*w(x)^2-g = 0

  

(2*c+2*(c^2+2*g)^(1/2))^(1/2), -(2*c+2*(c^2+2*g)^(1/2))^(1/2), I*(2*(c^2+2*g)^(1/2)-2*c)^(1/2), -I*(2*(c^2+2*g)^(1/2)-2*c)^(1/2)

 (1)
> 

#
# The above four "roots" will also be produced by the
# dsolve() command, since they are all valid solutions
# for the ODE
#
  [ dsolve( q ) ];

[w(x) = (2*c+2*(c^2+2*g)^(1/2))^(1/2), w(x) = -(2*c+2*(c^2+2*g)^(1/2))^(1/2), w(x) = -I*(2*(c^2+2*g)^(1/2)-2*c)^(1/2), w(x) = I*(2*(c^2+2*g)^(1/2)-2*c)^(1/2), w(x) = 2*JacobiSN((1/2)*(2*(c^2+2*g)^(1/2)-2*c)^(1/2)*x+_C1, I*(-((c^2+2*g)^(1/2)*c-c^2-g)*g)^(1/2)/((c^2+2*g)^(1/2)*c-c^2-g))*g/(g*((c^2+2*g)^(1/2)-c))^(1/2)]

 (2)
> 

#
# One could select the constant solutions provided by
# dsolve, using
#
  select(j->not has(rhs(j), x), [dsolve(q)]);

[w(x) = (2*c+2*(c^2+2*g)^(1/2))^(1/2), w(x) = -(2*c+2*(c^2+2*g)^(1/2))^(1/2), w(x) = -I*(2*(c^2+2*g)^(1/2)-2*c)^(1/2), w(x) = I*(2*(c^2+2*g)^(1/2)-2*c)^(1/2)]

 (3)
> 

#
# Alternatively if one wants a solution where w(x)
# actually varies with the independent variable 'x'
# then one coud use
#
  select(j->has(rhs(j), x), [dsolve(q)]);

[w(x) = 2*JacobiSN((1/2)*(2*(c^2+2*g)^(1/2)-2*c)^(1/2)*x+_C1, I*(-((c^2+2*g)^(1/2)*c-c^2-g)*g)^(1/2)/((c^2+2*g)^(1/2)*c-c^2-g))*g/(g*((c^2+2*g)^(1/2)-c))^(1/2)]

 (4)
> 

  plot( eval(rhs(%[]), [_C1=1, c=1, g=1.5]), x=-10..10);

 
> 

 

NULL

Download odeProb2.mw

Not clear exactly sure...

tomleslie 13836
September 29 2022
0 0

@mmcdara 

 what point you are trying to make here.

What do you mean by  "eq1 and eq2 can be considered as "coupled", although "simultaneous" might be a better designation"?
Either the solution f1 of eq1 is a solution of eq2 and the solution of eq2 is also a solution of eq1, in what case eq1 anq eq2 could be said "compatible"; instead they are not and eq1 and eq2 are "uncompatible".

I don't like the use of the term "coupled" to describe these equations, Since they can be solved "simultaneously", with the same (formal - so fairly useless!)  solution satisfying both equations, I'd prefer to refer to them as "simultaneous". If you want to use the term "compatible", that too is absolutely fine with me.

I am only surprised that you claim these equations are "uncompatible". A reduced version of the original worksheet, showing that the (formal - so fairly useless!) solution returned by pdsolve([eq1, eq2]) does solve both equations is attached

> 

restart
> 

eq1 := (2*(r^2+a^2*cos(theta)^2))*(M*r-(1/2)*a^2-(1/2)*r^2)*(diff(f(r, theta), r, theta))+(2*(a^2*(M-r)*cos(theta)^2-M*r^2+a^2*r))*(diff(f(r, theta), theta))

2*(r^2+a^2*cos(theta)^2)*(M*r-(1/2)*a^2-(1/2)*r^2)*(diff(diff(f(r, theta), r), theta))+2*(a^2*(M-r)*cos(theta)^2-M*r^2+a^2*r)*(diff(f(r, theta), theta))

 (1)
> 

eq2 := sin(theta)*(r^2+a^2*cos(theta)^2)*(diff(f(r, theta), theta, theta))-cos(theta)*(diff(f(r, theta), theta))*(a^2*cos(theta)^2-2*a^2-r^2)

sin(theta)*(r^2+a^2*cos(theta)^2)*(diff(diff(f(r, theta), theta), theta))-cos(theta)*(diff(f(r, theta), theta))*(a^2*cos(theta)^2-2*a^2-r^2)

 (2)
> 

fsol := pdsolve([eq1, eq2])

{f(r, theta) = _F1(r)+(Int((r^2+a^2*cos(theta)^2)/((cos(theta)+1)^(1/2)*(cos(theta)-1)^(1/2)), theta))*_C1/(2*M*r-a^2-r^2)}

 (3)
> 

pdetest(fsol, eq1); pdetest(fsol, eq2); pdetest(fsol, [eq1, eq2]); simplify(eval(eq1, fsol)); simplify(eval(eq2, fsol))

0

  

0

  

[0, 0]

  

0

  

0

 (4)
> 

 

Download pdIssue2.mw

Oh dear...

tomleslie 13836
September 28 2022
0 0

2-D input again!?!

Two functional versions using 2D input are given in the attached.

The first is 2D-input code wrtten more-or-less "from scratch". I really don't see the point of doing this, because it looks like nothing you will see in a math textbook. So what exactly is the point of doing it?

If you feel that you have to use 2D-input then never, ever copy/paste 1D-code! This is pretty much guaranteed never to work. A better approach is to select all of the 1D code, then, from the menus, use Format->Convert To->2-D Math Input.

The second execution group in the attached is the result of this conversion on my original 1D code. Less readable, than writing the 2D code from scratch above - but it works. Again does it look like anything you will see in a math textbook?

You also ask

What does the " 'nolist' " do?

Try reading the help for the command entries() with the option 'nolist'

> 

restart

with(LinearAlgebra)

A := Matrix(4, 4, {(1, 1) = 0, (1, 2) = 1, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 1, (2, 4) = 0, (3, 1) = 1, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = -1})

IMat := IdentityMatrix(4)

solve({entries(`~`[`=`](MatrixPower(A, n), IMat[]), 'nolist')}, allsolutions = true)[]

about(_Z1)

Matrix(%id = 36893488147927702092)

  

n = 6*_Z1

  

Originally _Z1, renamed _Z1~:
  is assumed to be: integer
 

  
> 

restart
with(LinearAlgebra)A := `<,>`(`<|>`(0, 1, 0, 0), `<|>`(0, 0, 1, 0), `<|>`(1, 0, 0, 0), `<|>`(0, 0, 0, -1))IMat := IdentityMatrix(4)solve({entries(`~`[`=`](MatrixPower(A, n), ` $`, IMat[]), 'nolist')}, allsolutions = true)[]about(_Z1)

n = 6*_Z1

  

Originally _Z1, renamed _Z1~:
  is assumed to be: integer
 

  
> 

 

Download solMat2.mw

You mean...

tomleslie 13836
September 24 2022
0 0

@Syed Asadullah Shah 

you want something like the attached?

> 

  restart:
  with(plots):
  odes:=  diff( f(eta), eta$3)+f(eta)*diff(f(eta), eta$2)-diff(f(eta),eta)^2-M*diff(f(eta),eta)=0,
          diff(g(eta),eta)-2*diff(f(eta),eta$2)-2*f(eta)*diff(f(eta),eta)=0;
  Mvals:=[0, 0.5, 1, 1.5, 2]:
  conds:= f(0)=0, g(0)=0, D(f)(0)=1, D(f)(inf)=0;
  for j from 1 by 1 to numelems(Mvals) do
       params:= [M=Mvals[j], inf=10]:
       sol[j]:= dsolve(eval([odes, conds], params), numeric);
       plt[j]:= odeplot( sol[j],
                         [ [eta, f(eta)],
                           [eta, g(eta)],
                           [eta, diff(f(eta), eta)]
                         ],
                         eta=0..10,
                         color=[black, red, blue]
                       );
   od:
   display( [seq( plt[j], j=1..numelems(Mvals))]);

diff(diff(diff(f(eta), eta), eta), eta)+f(eta)*(diff(diff(f(eta), eta), eta))-(diff(f(eta), eta))^2-M*(diff(f(eta), eta)) = 0, diff(g(eta), eta)-2*(diff(diff(f(eta), eta), eta))-2*f(eta)*(diff(f(eta), eta)) = 0

  

f(0) = 0, g(0) = 0, (D(f))(0) = 1, (D(f))(inf) = 0

  

 
> 

 
   
 

 

Download odeProb2.mw

Try...

tomleslie 13836
September 22 2022
0 0

reading the help page at ?tabulate

First 10 11 12 13 14 15 16 Last Page 12 of 207