I thought Id' solve this problem by tightening up the accuracy on pdsolve(), reducing the time range and then using the returned solution module's 'value' method as input to a straightforwar numerical integration.

The value() method seems to correctly return numeric values for u(0.5,t) over the complete time range which I specified. However the Int() calculation returns 'unevaluated' when the integration range is much more than 0..0.1.

So in the following, why does

evalf(Int(g, 0..0.1)); #work, but
evalf(Int(g, 0..0.2)); #returns unevaluated

Download heatPDE.mw

for n positive, the term 1+1/n is a;ways more tham 1.

So (1+1/n)^n is always more than 1^n ------ and 1^n is always 1

So if I perform the sum(`1, n=0..infinity), then the amswer will be infinity

Therefore the sum( 1+1/n)^n , n=0..infinity) will be rather more than infinity

Is this response useful?

You seem to think that the Kronecker Delta function and the Dirac Delta function are synonymous - Errrrrrr No

The Kronecker Delta function is a function of two variables, defined as f(x,y)=piecewise( x=y,1,0).

So you can justifiably expect that f(x,0)=1 if x=0 (or f(0,y)=1 if y=0). However the integral of the Kronecker Delta function is always zero -because (sloppy mathematics!) it has no spatial extent. Think of its integral as a rectangle whose "height\" is 1 and whose "base" is 0 - what is the area? A similar (sloppy) argument will convince you that the product of a |Kronecker Delta and any "sensible" function will always be zero.

Now the (rather naive) definition of the Dirac Delta function is that it is infinite when x=0, and zero otherwise. With this simple-minded definition (sloppy mathematics again), its integral is a rectangle of infinite "height" and zero "base", so its area is zero*infinity - and you get to take your pick what value that should be.

A "more careful" definition of the Dirac Delta function is as the limit as eps->0 of the function which extends from -eps/2 to eps/2 and has height eps. Notice that the "area" of this function is always 1
 

For instance, in the example you provide

define U(x,t)=diff(C(x,t),x) and substitute in the original PDE

then essentially rinse and repeat

At the end you will have four first order PDEs in four dependent functions.. For your example, you could use, the attached. Obviously this code would have to be improvedextended to handle derivatives wrt all variables, mixed derivatives, etc, but so far as I know, this approach always works

Download pdeSub.mw

I'm pretty sure that the attached does what you want, but it wasn't exactly clear in what form you wanted the output - so I guessed.

This will generate an array of plots, and an array of 'integral sums' which yu can access in any appropriate way

loop.mw

so far as I can tell.

You just have to accept that just because you can write an equation does mean that it can be solved - at least analytically.

Given numeric values for all the parameters, it may be possible to solve it numerically

Too many syntax errors to list.

I have fixed all those I can find. The attached now runs and allows you to examine the values of the matrix u_new at each timestep

finElem.mw

It would seem that you cannot first generate a matrix, and the subsequently set its attributes using a convert command. You have to generate the matrix with the required attributes.

Must admit I didn't know that.
However, in your case, it is relatively easy to  generate the required matrix with shape=symmetric and datatype=float. In which case all generated eigenavlues will be real and determining min/max is pretty trivial.

realEV.mw

 I assumed that theta[1], theta[2], phi[1], phi[2] are all real. Defining 0<p<1 means that it is real also.

I couldn't persuade the built-in Optimization:-Maximize command to produce a solution, but I did get a solution using the add-on DirectSearch() package available free from the Maple application centre. This was

[1., [p = 1., `&phi;`[1] = 1.20193938785539, `&phi;`[2] = 1.12831141336330, `&theta;`[1] = 6.28318530006250, `&theta;`[2] = 6.28318530717958], 489]

Rather obviously, if the maximum is 1, obtained at p=1, then the term multiplied by sqrt(1-p) contributes nothing - so phi[1], phi[2] could actually be any value! In addition although DirectSearch returns 2*Pi for both theta[1]=0, theta[2]=0, p=1 would return the same maximum.

LifeLesson: always treat the output of optimisers with care!!!

I have uploaded the attached including Maple output, becuase you (probably) won't be able to re-run all of it (unless you do happen to have the DirecSearch() package loaded)

opt.mw

Add the option numpoints=1000 to the first plot command, as in

plot(result3, f = 100 .. 5000, numpoints = 1000)

The equation T=.... is quadratic in the variable 'p', so using solve will return two answers for p

It is however possible to generate a Taylor series for each of these possible solutions about the point b=0.

See the attached

Download taylorProb.mw

I've just tried this in 64-bit Maple 18.02 running on 64-bit Win7.

GIves exactly the same  answer as Maple2017.

No error messages, no issues, just the solution.

Try posting a worksheet (with output) which exhibits the problem using the big green up-arrow in the MaplePrimes toolbar

I think most of your confusion is arises because you seem uncertain about the distinction between defining a function and calling it.

See if the attached makes anything clearer, please read the comments carefully

pdeProb.mw

From the 'subs' help page

The action of substitution is not followed by evaluation by default.  There are two ways of achieving a full evaluation:

So either of the following will work in your case

restart;
Ex1 :=expand(diff(f(x+y)+x*f(x-y), x)^2)+expand(diff(f(x-y)-y*f(x+y), x)^2):
eval(subs(D(f) = (t -> 1-t), Ex1));
subs[eval](D(f) = (t -> 1-t), Ex1);

Note that there are subtle differences in the two methods of forcing the evaluation. Read the advice on the 'subs' help page, and check the examples carefully

 As written, you only have four independent variables - a,b,c,d (which may, or may not) depend on further parameters  - e, f.

So

solve( P(a,b,c,d,e,f)=0., [a, b, c, d] )

ought to work. Stated in this way, with a,b,c,d identified as the independent variables, then any/all of these may be returned explicitly, or as expressions involving 'e' and/or 'f''