## 5284 Reputation

7 years, 121 days

## @Kitonum  Thanks for your feedback...

@Kitonum

Concerning point 2 "But how will you solve a similar problem if the number  n  is large", this was the meaning of my my writting "regardless of its computational efficiency".
I was aware that that "my" solution will present serious problems of computational efficiency when n is large.

## @Kitonum No offense intended, but I...

@Kitonum

No offense intended, but I think your program has two shortcomings:

• firstly, it seems too complex,
• secondly, it's not a good idea (IMO) to specify a range as the second argument of Partition.

For instance you write

```n := 9:
Partition(n, 0..n, {2,3,6,9});
```

because you know that there is a one-element partition of 9, but why do you write

`Partition(8, 2..8, {2,3,6,9});`

Why?
In fact a correct range for partitionning any number n should have been

`Partition(n, 0..n, {2,3,6,9});`

, which proves the specification of this range is useless.
Unfortunately Partition requires the first operator of the range to be a strictly positive integer

```Partition(8, 0..8, L);
Error, (in Partition) invalid input: k_Partition expects its 2nd argument, k, to be of type posint, but received 0
```

Note that Partition does not always return the correct answer:

```Partition(8, 1..8, L);
[, [2, 6], [2, 3, 3], [2, 2, 2, 2]]
```

To answer the OP's question this simpler code below seems[*]  enough (regardless of its computational efficiency) and doesn't require the specification of any range.

```f := proc(n::posint, L::set)
local P := combinat:-partition(n):
select((x -> is ({x[]} subset L)), P, L)
end proc:

L := {2, 3, 6, 9}:
f(8, L);
[[2, 2, 2, 2], [2, 3, 3], [2, 6]]
f(9, L);
[[2, 2, 2, 3], [3, 3, 3], [3, 6], ]
f(3, L);
[]
f(1, L)
[]
```

[*]  we both implicitely assume then the number to be partionned was a positive integer. What does the OP want if this number is negative and the set of number contains negative numbers, for instance:

```n := -1:
L := {-2, 1}:
# then n writes
n := L + L```

Your procedure  Partition may be used in a lot of ways and is then more general than the simple code above.
Nevertheless I would suggest you to modify Partition in order to eliminate this unnecessary range argument  range argument and correct the "Partition(8, 1..8, L) error"

I hope I wasn't too blunt in my comment.

@dharr

I almost missed your reply because the contributor counter was stuck at 1.
That is a good tip indeed, thanks

## @dharr C <= 0: yes of course,...

C <= 0: yes of course, what was I thinking?
Only later did I think to check whether the matrix was symmetric.
Finally, I don't understand either why my approach doesn't work after correcting the first condition.

Even if I can't execute your code with Maple 2015, I carefully read it: you used a very beautiful approach and I vote up.

## @Rouben Rostamian   For incormatio...

@Rouben Rostamian

For incormation: I believe this question has connections to this one
https://www.mapleprimes.com/questions/236730-Solving-Schrodinger-Equation

## @somestudent It's very tiresome...

It's very tiresome that you spend your time modifying your initial question according to my answers and that you only divulge information on an ad hoc basis.
I feel like I'm wasting my time.

So this is my last answer.
Notice I rewrote your code in a much clever way, for doing copy-paste to define pole1, pole23re and pole23im is a bit stupid (no offense intended).

Your code ends with this command

`maximize(MyFilterABS, f = 0 .. .5, location = true)`

which is a nonsense for MyFilterABS depends on x, y, w, and f and maximize requires that all the quantities but those you want to maximize with respect to, are given a numerical value.
So I've taken the initiative to replace MyFilterABS by its evaluation when x, y, w given numerical values.

Last point: your sentence "both points are inside the unit circle, both have a maximum at f = 0.5 but for the second one there is one maximum more" is unintelligible.

 > restart
 > with(ListTools): with(Optimization): with(plots):
 > MyFilter := (z-1)^3/(x*z^2+z^3+y*z-w); (1)
 > d := denom(MyFilter); alias(alpha = RootOf(d, z)); evala(Factor(d, alpha))   (2)
 > poles := [allvalues(alpha)]: sel   := select(`not`@has, poles, I); (3)
 > [indets(sel)[]]: select(type, %, `^`); ToReplace := op(1, %);  (4)
 > SimplificationRule := ToReplace=A; Simpler_Poles := eval(poles, SimplificationRule): ReIm := map(sp -> [selectremove(`not`@has, sp, I)] *~ [1, -I], Simpler_Poles);  (5)
 > PPZero := solve(numer(MyFilter) = 0, z); (6)
 > cs := [cos(2*Pi*f), sin(2*Pi*f)] (7)
 > MyFilterABS := mul(sqrt((cos(2*Pi*f)-z)^2+sin(2*Pi*f)^2), z in [PPZero]) / mul(sqrt(add((cs-~ri)^~2)), ri in ReIm) (8)
 > test_1 := eval(eval(MyFilterABS, (rhs=lhs)(SimplificationRule)), [w = .1, x = .3, y = .7]): plot(test_1, f=0..1, gridlines=true);  # I see a minimum at f=1/2, not a maximum! test_2 := eval(eval(MyFilterABS, (rhs=lhs)(SimplificationRule)), [w = 0.2e-1, x = -.2, y = .1]): plot(test_2, f=0..1, gridlines=true)  > # Two maxima found maximize(test_1, f=0..1, location) (9)
 > # One maximum found maximize(test_2, f=0..1, location) (10)

 >

To find the maxima in a larger f range:

```# Ten maxima found

maximize(test_1, f=0..5, location);
numelems(%);
16.32115142, {[{f = 0.2935930416}, 16.32115142],

[{f = 0.7064069584}, 16.32115142],

[{f = 1.293593042}, 16.32115142],

[{f = 1.706406958}, 16.32115142],

[{f = 2.293593042}, 16.32115142],

[{f = 2.706406958}, 16.32115142],

[{f = 3.293593042}, 16.32115142],

[{f = 3.706406958}, 16.32115142],

[{f = 4.293593042}, 16.32115142],

[{f = 4.706406958}, 16.32115142]}
10
# Five maximum found

maximize(test_2, f=0..5, location);
6.060606061, {[{f = _Z29 + 0.5000000000}, 6.060606061]}
{_Z29}
is used in the following assumed objects
[_Z29] assumed AndProp(integer,RealRange(0,4))
```

## @apronya So much the better...

So much the better

## @somestudent complexfunction is a f...

complexfunction is a function of x, y, f.
Thus the maximizer is a triple (x, y, f=1) and x --> -oo, y --> +oo (obvious result).
The value of complexfunction at this maximizer is equal to +oo (also obvious).

When you write the constraint

`{Maximize(complexfunction, {f <= 1}) = 1}`

for the "outer" Maximize, you write in fact the constraint

`{+oo <= 1}`

This is always false and the constraint will never be reached in the "outer" Maximize ... giving you no solution.

I advise you to reformulate your problem correctly before considering solving it with Maple.

Could_it_be_this.mw

On what base do you claim that  what @tomleslie presents is not based on a Finite Difference Method?
It is not because you do not see the code behind pdsolve that pdsolve doesn't use FDM!

Type this in some wotksheet

`help(pdsolve[numeric])`

In section Details section, you will read this
The default method uses a second order (in space and time) centered, implicit finite difference scheme to obtain the computed solution.

`help(pdsolve[numeric][method])`

In section Description it's written :
In addition to the default interface to pdsolve/numeric, which uses the default theta methods (the theta family of methods), additional functionality is available for use of numerical pdsolve for educational purposes.
...
additional functionality is available for use of numerical pdsolve for educational purposes.
[the underlining is mine]
This additional functionality includes the use of 11 standard methods, specification of numerical boundary conditions, specification of startup schemes for two-stage methods, and extensibility for your custom methods

The theta-method is a family of numerical scheme controlled by a parameter "theta".
For 1D-space heta equation the choice theta=1/2 corresponds to the Crank-Nicolson scheme/

• https://folk.ntnu.no/leifh/teaching/tkt4140/._main065.html
• https://www.math.pku.edu.cn/teachers/lizp/courses/Numerical_PDE/Numerical%20PDE%20Lecture%20Slides/numpde_lecture_5_c2.pdf
• https://itp.uni-frankfurt.de/~rezzolla/lecture_notes/2007/ECT_hyperbolic_PDEs_0807.pdf

In @tomleslie's code you can force the Crank-Nicolson scheme by setting

``` PrVals:=[0.71, 1.00, 3.00, 7.00]:
colors:=[red, green, blue, black]:
for j from 1 to numelems(PrVals) do
pars1:=`union`( pars, {Pr=PrVals[j]}):
pdSol:= pdsolve( eval([pdes], pars1),
eval([conds], pars1),
numeric,
method=CrankNicholson,
);
plt[j]:=pdSol:-plot( u(y,t), t=2, y=0..inf, numpoints=200, color=colors[j]);
od:
plots:-display( [seq(plt[j], j=1..numelems(PrVals))]);
Error, (in pdsolve/numeric/par_hyp) the CrankNicholson method can only be used for a single PDE

```

The error message is clear.
I guess (see above) that the development team didn't bother to implement the Crank-Nicolson scheme for pde systems for two reasons:

• Crank-Nicolson scheme is mainly for educational purposes.
• and theta-method is more general.
So why do the job twice?

Note that if you really want to implement the Crank-Nicolson scheme for pde systems, you can do it yourself, see here

`help(pdsolve[numeric][pdemethods])`

As your question was "How to get same graph from maple with finite difference method for differential equations " I consider that @tomleslie has perfectly answered it (and I vote up).

## Could you give the expressions of g and ...

Could you provide us the expressions of g and f?
It is absolutely necessary for one can build extremely simple examples where no solution does exist:

 > restart:
 > with(Optimization):
 > # Your expression takes this form in 1D # Maximize(f(x), {Maximize(g(x,f),{0 < f, f < 1}) < Thres},initialpoint = {x=1},iterationlimit = 9999);
 > # A simple example for which no solution exist f := x -> x^2 (1)
 > g := unapply(f(x)/x, x) (2)
 > max_1 := Maximize(g(x),{0 <= x^2, x^2 <= 1}, initialpoint={x=2}) (3)
 > Thres := 1/2 (4)
 > Maximize(f(x), {max_1 <= Thres}, initialpoint = {x=1});
 >

## @AHSAN You can customize "my&q...

@AHSAN

You can customize "my" code the way you want to get the rendering which pleases you.
But writting a procedure which accounts for all your desires is, IMO, a total waste of time.

Nevertheless here is an evolution of my previous code.

 > > > ind := [indets(Nu, name)[]] (1)
 > # You can't assign beta here because beta is no longer a name but a list in the expression of Nu # Display Nu to understand what happens # beta := [0.1, 0.3, 0.5, 0.7, 0.9]; # Nu
 > # Define a function which maps the indeterminates onto Nu: nu := unapply(Nu, ind): # Define a procedure which draws a Univariate BarGraph (only one parameter change) UnivariateBarGraph := proc(                             param::name                             , values::list                             , others::list(`=`)                             , col::{name, string}                             , {                                 dx::positive:=0                                 , shadow::boolean:=false                                 , Halign::anything:=NULL                                 , Hpos::numeric:=0                               }                           )   local instances, N, r, h, nuval, sy, sx, Tr:   if `not`(is(param in ind)) then     error cat("param ", sprintf("%s", parameter), " is not an indeterminate of Nu"):   end if:       instances := i -> [others[], param=values[i]]:   N         := numelems(values):   r         := proc(x, y, H, c, {pol::boolean:=false})                  plottools:-rectangle([x-H, 0], [x+H, y], color=c, `if`(pol, style=polygon, NULL))                end proc:   if dx = 0 then     h := (values-values)/4:   else     h := dx   end if:      nuval := seq(nu(eval(ind, instances(i))[]), i=1..N):   sy    := 1/4:   sx    := sy/max(abs~([nuval])):   Tr    := plottools:-transform((x, y) -> [x+y*sx, y*sy]):   plots:-display(     seq(       r(values[i], nuval[i], h, col)       , i=1..N     ),     seq(       `if`(shadow, Tr(r(values[i], nuval[i], h, "LightGray", pol=true)), NULL)       , i=1..N     ),     seq(       plots:-textplot(         [           values[i]+Hpos           , nuval[i]           , nprintf("%1.3e", nuval[i])         ]         , align=`if`(                    Halign = NULL                    , `if`(nuval[i] > 0, above, below)                    , `if`(nuval[i] > 0, {above, Halign}, {below, Halign})                  )       )       , i=1..N     )     , labels=[typeset(param), "Nu"]     , axes=boxed     , axis=[tickmarks=values]     , view=[         min(values)-3*h..max(values)+3*h         , (1.1*min-0.1*max)(nuval).. (1.1*max-0.1*min)(nuval)       ]   ) end proc:
 > UnivariateBarGraph(beta, [seq(-0.9..-0.1, 0.2)], [Q = .4290, lambda = -0.388e-1, k = .3, x = 0, Br = .1], "Chartreuse") > UnivariateBarGraph(beta, [seq(-0.9..-0.1, 0.2)], [Q = .4290, lambda = -0.388e-1, k = .3, x = 0, Br = .1], "Chartreuse", shadow=true) > res_1 := UnivariateBarGraph(            beta            , [seq(0.1..0.9, 0.2)]            , [Q = .4290, lambda = -0.388e-1, k = .3, x = 0, Br = .1]            , "Gold"            , dx=0.025            , Halign=left            , Hpos=0.02          ): res_2 := UnivariateBarGraph(            beta            , [seq(0.1..0.9, 0.2)]            , [Q = .4290, lambda = -0.388e-1, k = .3, x = 0, Br = .2]            , "Brown"            , dx=0.025            , Halign=right            , Hpos=-0.02          ): plots:-display(   plottools:-translate(res_1, -0.025, 0)   , plottools:-translate(res_2, +0.025, 0)   , title = cat("Q = .4290, lambda = -0.388e-1, k = .3, x = 0 \n Gold: Br = .1, Brown: Br = .2")   , labels=[beta, "Nu"]   , axis=[gridlines=[seq(0.1..0.9, 0.2)]]   , size=[800, 400]   , view=[0..1, 0..-0.025] ) >

If tou wnat to change the indication at he end of a bar, you must add an optional parameter to procedure UnivariateBarGraph (all these optional parameter are gathered in curly braces in the sequence of parameters of the procedure).
For instance:

` , HeadTitle::string:=""`

and change

```plots:-textplot(
[
values[i]+Hpos
, nuval[i]
, nprintf("%1.3e", nuval[i])
]
, align=`if`(
Halign = NULL
, `if`(nuval[i] > 0, above, below)
, `if`(nuval[i] > 0, {above, Halign}, {below, Halign})
)
)
```

by

```plots:-textplot(
[
values[i]+Hpos
, nuval[i]
]
, align=`if`(
Halign = NULL
, `if`(nuval[i] > 0, above, below)
, `if`(nuval[i] > 0, {above, Halign}, {below, Halign})
)
)
```

## Some mistakes in your code.....

@AHSAN

... the most dramatic is that you must not assign beta to a value or a list where tou did that.

Look to the attached file

 > > > ind := [indets(Nu, name)[]] (1)
 > # You can't assign beta here because beta is no longer a name but a list in the expression of Nu # Display Nu to understand what happens # beta := [0.1, 0.3, 0.5, 0.7, 0.9]; # Nu
 > # Define a function which maps the indeterminates onto Nu: nu := unapply(Nu, ind): # Define a procedure which draws a Univariate BarGraph (only one parameter change) UnivariateBarGraph := proc(param::name, values::list, others::list(`=`), col::{name, string})   local instances, N, r, h, nuval:   if `not`(is(param in ind)) then     error cat("param ", sprintf("%s", parameter), " is not an indeterminate of Nu"):   end if:   instances := i -> [others[], param=values[i]]:   N := numelems(values):   r := (x, y, H, c) -> plottools:-rectangle([x-H, 0], [x+H, y], color=c):      nuval := seq(nu(eval(ind, instances(i))[]), i=1..N):   h := (values-values)/4:   plots:-display(     seq(       r(values[i], nuval[i], h, col)       , i=1..N     ),     seq(       plots:-textplot(         [           values[i]           , nuval[i]           , nprintf("%1.3e", nuval[i])         ]         , align=`if`(nuval[i] > 0, above, below)       )       , i=1..N     )     , labels=[typeset(param), "Nu"]     , axes=boxed     , axis=[tickmarks=values]     , view=[         min(values)-3*h..max(values)+3*h         , (1.1*min-0.1*max)(nuval).. (1.1*max-0.1*min)(nuval)       ]   ) end proc:
 > UnivariateBarGraph(beta, [seq(0.1..0.9, 0.2)], [Q = .4290, lambda = -0.388e-1, k = .3, x = 0, Br = .1], "Chartreuse") > UnivariateBarGraph(Br, [seq(1..9, 2)], [Q = .4290, lambda = -0.388e-1, k = .3, x = 0, beta = .1], "Chartreuse") > # Why are the two graphs identical? # # When x=0 and beta=0.1 as Br is variable, or Br=0.1 as beta is variable, # the coefficients of Br, respectively beta, are the same: Nu_Br   := coeff(eval(Nu, [x=0, beta=0.1]), Br): Nu_beta := coeff(eval(Nu, [x=0, Br=0.1]), beta): simplify(Nu_Br-Nu_beta); (2)
 > # Thus the two plots are the same. # # Draw this one to vince you there is no error in procedure UnivariateBarGraph: UnivariateBarGraph(Br, [seq(1..9, 2)], [Q = .4290, lambda = -0.388e-1, k = .3, x = 0, beta = .2], "Chartreuse") >

## @AHSAN Nu is a function of Q, beta ...

Nu is a function of Q, beta and lambda:

```indets(Nu, name);
{Q, beta, lambda}
```

You write: "I want to draw bar chat for different value of beta [0.1,0.3,0.5,0.7,0.9) between Nu and y ( Nu is on vertical axies and y is on horizontal axis) and crosspedning to the values of beta the value of Lambda and Q also changes".

What is y (which is not part of the definition of Nu) ?

One doesn't draw (a) bar chart for different value of beta  between Nu and y... one draws (a) bar chart of N (and y if any) for different value of beta : I still do not understand what you want.

It would be helpfull if you coukld provide a sketch of what you want (use PowerPoint or do this drawing by hand and scan it).

## The content of your mw file has nothing ...

The content of your mw file has nothing to do with data that can be presented under BarChart form.
Basically it contains an expression N which depends upon beta, lambda, Q and that's all.
See here

`help(BarChart);`

Maybe you should explain clearly what your desired BarChart is aimed to represent?

## @arashghgood  Could you provide yo...

@arashghgood

Could you provide your Maple code?
Just click on the big green up-arrow in the menu bar.

﻿