This can be done using MathML coding/
Example

restart

 # example 1

`#msup(mo(u),mo(2))` := 3;
`#msub(msup(mo(u),mo(2)),mo(h))` := 4;
3*`#msup(mo(u),mo(2))` + a*`#msub(msup(mo(u),mo(2)),mo(h))`;

# example 2
 
`#msup(mo(u),mo(2))` := '`#msup(mo(u),mo(2))`':
`#msub(msup(mo(u),mo(2)),mo(h))` := ' `#msub(msup(mo(u),mo(2)),mo(h))`':

f := (x, y) -> 1/log(x)+y^3:
f(`#msup(mo(u),mo(2))` , `#msub(msup(mo(u),mo(2)),mo(h))`);

But manipulating such names is not that easy.
Using aliases can help simplifying all this:

restart

# alias definitions

alias(`#msup(mo(u),mo(2))`=__x, `#msub(msup(mo(u),mo(2)),mo(h))`=__y);
       #msup(mo(u),mo(2)), #msub(msup(mo(u),mo(2)),mo(h))

# example 1

3*__x+a*__y;

# example 2
​​​​​​​
f := (x, y) -> 1/log(x)+y^3:

f(__x, __y);  # smart output
lprint(%);    # detailed output

Here are a some improvements to what tomleslie already did:

1. Use of a simple strategy distribute the values of the contour according to the distribution of the 25x25 values
   (based on Statistics:-Quantile)
2. Use of CurveFitting:-ArrayInterpolation to smooth the contour plot.
3. Use of Explore to visualize how the rendering is affected by yhe changes of the number K of contours, the size N of the interpolation grid and the method T of interpolation.

Clever (?) strategies for choosing the contours could be used, while Maple offers (ImageTools package) other strategies to "smooth" the rendering of the plot.
contours.mw

Look at this
Student:-LinearAlgebra:-LinearSolveTutor( M );
Next click on the button "Next", look to the top right textfield what is done an observe the result in the main window.+

I don't understand how your Mi matrix could be useful to you.
If you really want to solve your system of equations using matrices (the LinearAlgebra package) you can do this

eqs := [ 5*x+4*y+3*z = 2, 7*x+6*y+4*z = 5, 7*x+3*y+5*z = -11];
vars := [indets(eqs)[]];
A, B := LinearAlgebra:-GenerateMatrix(eqs, vars)
sol := LinearAlgebra:-LinearSolve(A, B)
#check
eval(eqs, vars =~ convert(sol, list))

But solve(eqs) will do the job as the same

What do you mean exactly by multivariate polynomial regression here are a different possibilities
Given P regressors X₁.., X_P and Q dependent variables Y₁.., Y_Q :

1. some say multivariate in the case P > 1 and Q = 1
   (the univariate case being P=Q=1)
2. others say multivariate if Q > 1 and (generally P > 1 too) 

Case 1 can be handled by "native" Statistics:-LinearFit procedure
Case 2 can or can't depending the type of regression model you are looking for:

1. the simplest thing is to run Q independent "multivariate regressions" of type 1 (P > 1 and Q=1);
   this will give you Q independent regression models.
2. but you might be looking for a unique combination of X₁.., X_P (or f₁(X₁.., X_P), ...  f_R(X₁.., X_P) in case of polynomial regression) that globally explains the Y₁.., Y_Q :
   1. if the Y₁.., Y_Q are correlated ways are regression on principal component factors or redundancy analysis
   2. if the  X₁.., X_P are correlated too you must turn to partial least regression

Situation 1 can be done with Statistics:-LinearFit again; situation 2 require specific development.

Here is an example for P = 3 and Q = 2 wher Y₁ and Y₂  are assumed independent.
Let me know if you are interested with Y₁.., Y_Q correlated (a weaker notion than dependancy)

example.mw

Step 1: type this

restart: 
P1:=plot(sin(x), x=0..Pi);
dataMat:=plottools:-getdata(P1);
whattype(dataMat);

dataMat is ia list of 3 elements:

1. "curve"
2. a list of ranges
3. a matrix

Converting this to a Matrix as you do won't help.
What you want is to "extract" this matrix alone.

# do this to see what P1 is made of
op(P1);
# and now this to extract the matrix of points
dataMat := op([1, 1], P1);
# then export
ExportMatrix( "test.mat", dataMat, target=Matlab);

For your 3D example:

restart:
z=sin(x)*cos(y);
P2:=plot3d(z);
op(P2);
dataMat := convert(op([1, 3], P2), Matrix);  # as op(...) is an Array, not a Matrix
ExportMatrix( "test.mat", dataMat, target=Matlab);

Done with Maple 2015

Ar the end of your worksheet type this

lprint(B[..,3])
lprint(B[..,4])

You will see that the third column of B has elements ot type Hfloat while the fourth is of "simple" type float.
This means that some computation in your worksheet is done with Hfloats (from the Hfloat help page : The presence of a hardware floating-point number in an expression generally implies that the computation will use hardware floating-point evaluation, unless the settings of Digits and UseHardwareFloats specify otherwise (see UseHardwareFloats).):

The simplest (?) trick to avoid this is to write this at the top of your worksheet (it's better to write the restart in a separate block) 

UseHardwareFloats:=false

I didn't try to seek where Hfloats are generated (I had some difficulties accommodating my eyesight to the colors you use :-) )

You write:
In an academic assignment, the load function acting on top of a beam contains dirac function and its derivatives.

No, the load function acting on the top of a beam contains  "concentrated forces" ( the one you call "Dirac forces") and density of forces ("distributed forces), never derivatives or Dirac functions.
Euler–Bernoulli_beam_theory

The first thing to notice is that "Dirac forces" do not exist: they are just a useful mathematical representation of the reality (a force always act on a domain of strictly positive measure otherwise the local pressure would be infinite). More of this a concentrated force has a finite magnitude while the amplitude of a dirac is not defined.

If you really want (but it is useless) to define a loading in terms of a set of Dirac(s), here is a simple trick

# a force of magnitude M at location x=a

F := A * diff(Heaviside(x-a), a)

Another trick would be to define a uniform density of force (A)  in the interval [a-eps, a+eps] such that 2*epsilon*A=M.

Why is this useless?
Because the shear force, which is a key element in the ode that governs the flexion of the beam, is the antiderivative of the loading (concentrated and distributed forces)... and thus integrating Dirac(s) gives Heaviside(s).


The attached file gives an illustration
Beam_for_fun.mw

If I'm not mistaken, here is the "first order" aymptotic expansion.
(This could be done by hand using the Laplace's Method)

The result is of the form

G(x) * N*exp(-N*sqrt(x^2+1)+N*arcsinh(1/x))

(G is given below and in the attached file).

Note: this asymptotic expansion goes to 0 if x is larger than a critical value equl to 0.6627 and goes to +oo is x < 0.6627.
This means this asymtotic expansion is wrong (at least) for x < 0.6627 (the integral vanishes for any positive x as N --> +oo)... either one should push the expansion to larger orders or I did something wrong

# the critical value of x verifies
fsolve(-sqrt(x^2+1)+arcsinh(1/x))
                     0.66274341934918158097

I do not know if Maple can do this without help
 

restart
f := (n, x) -> Int(exp(-n*(x*cosh(t)+t)), t = 0 .. infinity);

# the previous relation is of the form

f__x(epsilon) = Int(g(t)*exp(phi__x(t)/epsilon), t = 0 .. infinity);

# with 
#   g(t)=1, 
#   phi__x(t)=-(x*cosh(t)+t)   (x is consider as a known parameter)
#   epsilon=1/n

# Rewritten this way one can apply the Laplace's Method:
#     g(t) and phi__x(t) are smooth functions
#     epsilon is a small positive parameter
# Doing this one must assume:
#     1/ that phi__x has a global maximum at some value t=c
#        (cf below)
#     2/ diff(phi__x(t), t$2) < 0 at t=c
#        The positiveness of x implies this condition is verified too 

phi__x := t -> -(x*cosh(t)+t) ;
c      := solve(diff(phi__x(t), t)=0, t);

diff(phi__x(t), t$2);
eval(%, t=c);         # negativeness of diff(phi__x(t), t$2) at t=c


# Laplace's Method:
# replace phi__x(t) by its taylor expansion up to order 2
# to get 

f__x(epsilon) = Int(exp(('phi__x(c)'+1/2*D[2]('phi__x(c)')*(t-'c')^2)/epsilon), t = 0 .. infinity);

h :=convert(taylor(phi__x(t), t=c, 3), polynom);

# verify this is identical to 

phi__x(c) +1/2*eval(diff(phi__x(t), t$2), t=c)*(t-c)^2;

simplify(%-h)
         
                               0

# note that phi__x(c) doesn't depend on t, thus 

f__x(epsilon) = exp('phi__x(c)'/epsilon)*Int(exp(1/2*D[2]('phi__x(c)')*(t-'c')^2/epsilon), t = 0 .. infinity);


# The integral could easily be computed by hand (a ERF function)
# but I use Maple
#    I set A = D[2]('phi__x(c)') (which is negative)
#    and C=c to avoid premature evaluation and keep the expression simple.


K := int(exp(1/2*A*(t-C)^2)/epsilon, t=0..+infinity) assuming A < 0, epsilon > 0;

# Finally

J := simplify( 
       exp(phi__x(c)/epsilon) 
       * 
       eval(K, [C=c, A=eval(diff(phi__x(t), t$2), t=c)]) 
     ) assuming x > 0;


# and thus the asymtotic expansion with regard to N:

Asympt_J := combine(asympt(eval(J, epsilon=1/N), N), exp);

AsymptoticExpansion.mw

In the plotsetup help page it is written :
          A list of colors cannot be provided to the plots[pointplot] or plots[pointplot3d] command.

Here is a workaround (Maple 2015)
 

restart:
List := [[[0, 0], 1], [[1, 2], 0], [[2, 3], 0], [[3, 1], 1]]:
C := n -> COLOR(RGB, op(`if`(List[n][2]=1, [1, 0, 0], [0, 0, 1])));

p := PLOT(seq( POINTS([List[i][1]], C(i)), i=1..nops(List) )):
# verify
plots:-display(p);

plotsetup(ps, plotoutput=cat(currentdir(), "/Desktop/test"), plotoptions=`color=cmyk`);
plots:-display(p);
plotsetup(default):

Here is the content of the file test.eps (the extension is automatically added).

As Mapleprimes doesn't accept to load eps files I opened the eps file with a ps wiever and exported it to a png file... but believe me, test.eps contains colors.

Hi,

I understand your dismay if you are not familiar with signal processing, DFT (Discrete Fourier Transform), and signal sampling.

You will find a few explanations in the attached file but I strongly advice you to read some book or course on the topic (I have a few that I used when I was a student but they are getting old and are in French).

Basically the DFT differs from the FFT on at least two points (I put aside considerations related to the sampling [=discretization] of the signal):

• signals must be causal (because they are assumed to be real signals!)
• the signal being observed on a finite period of time (or space), one must decide what this signal is outside of this observation window:
  • either it is considered as perioc
  • either it is considered as aperiodic

These considerations implie the FFT and the DFT geberally give different "raw" results.
But knowing the assumptions used by the DFT one can recover the result a FFT gives.

In your case, knowing these assumptions
(cf formulas used in DiscreteTransforms[FourierTransform] help page)
helps to explain the discrepancy betwen the exact and the numerical solutions, and also enable "adjusting" the numerical one.

DFT.mw

Hi, 
Here is a solution

restart:
interface(version)
Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895

eq__1 := W__mn^2*n__1+L__11*U__mn+L__12*V__mn+L__13*W__mn = C__1*`ΔT`:
eq__2 := W__mn^2*n__2+L__21*U__mn+L__22*V__mn+L__23*W__mn = 0:
 
with(LinearAlgebra):
A, b := GenerateMatrix([eq__1, eq__2], [U__mn, V__mn]);
sol := collect~(LinearSolve(A, b), W__mn):

# print solution
< [U__mn, V__mn] > = sol;

# print coefficients
[seq(a__||k, k=0..2)] =~ coeff~(sol[1], W__mn, [0, 1, 2]);
[seq(b__||k, k=0..2)] =~ coeff~(sol[1], W__mn, [0, 1, 2]);

solution.mw

I don't know what you expect from Maple  (what do you doubt of as the name of your file seems to mean)?
Given the abstract definitions of RV x and beta, you have absolutely no chance of finding (if that's what you were hoping for)  any closed form expression for q__ER nor Pi__ER.

Here is the farthest you can go with Maple (and probably with any other CAS or even by hand)

restart:
with(Statistics):

#---------------------------------------------------x

d := Distribution(
       PDF = (proc (t) options operator, arrow; f__D(t) end proc), 
       Quantile = (proc (q) options operator, arrow; int(f__D(t), t = -infinity .. q) end proc)
    ):

x := RandomVariable(d);

F := t -> Quantile(x, t);
F__c := t -> 1- F(t);


t -> Statistics:-Quantile(x, t);
t -> 1 - F(t);


#---------------------------------------------------beta

d1 := Distribution(
        PDF = (t -> g__beta(t)), 
        Quantile = (q -> int(g__beta(t), t=-infinity..q))
      );

beta := RandomVariable(d1);

G := t -> Quantile(beta, t);
G__c := t -> 1- G(t);


#---------------------------------------------------Pi__ER

Pi__ER := beta__R -> piecewise(
            q/(G__c(beta__R)+m*G(beta__R))-x >= 0, 
               ((p-s)*(G__c(beta__R)+m*G(beta__R))-(c-s)*q-(p-s)*(1-m)*G__c(beta__R))*x, 
            `and`(x-q/(G__c(beta__R)+m*G(beta__R)) > 0, q/G__c(beta__R)-x > 0), 
               (p-s)*q-(c-s)*q-(p-s)*(1-m)*G__c(beta__R)*x, 
            x-q/G__c(beta__R) >= 0, 
               (p-s)*q-(c-s)*q-(p-s)*(1-m)*q);



# Differentiate first and take the mean after

dp := diff(Pi__ER(beta__R), q); 

# write the conditions under simpler forms

UV := [ op(1, rhs(op(1, dp))) = U, -op(1, rhs(op(-2, dp))) = V ];

dpUV := eval(dp, UV);
  
# simplify the writing of dpUV

dpUV := piecewise(_R <= U, op(2, dpUV), _R < V, op(4, dpUV), op(6, dpUV))
       

with(Student[Calculus1]):

mean_dp := Mean(dpUV) assuming U < V:
mean_dp := eval(%, int=Int):
mean_dp := add(rhs~(map(Rule[`c*`], [op(mean_dp)])));
             

mean_dp := (-c+s)*(Int(f__D(_t)*_t, _t = -infinity .. U))+(p-c)*(Int(f__D(_t), _t = U .. V))+(m*p-m*s-c+s)*(Int(f__D(_t), _t = V .. infinity))

              
# let make "q" to appear (if you want it)

Krule  := op(1, denom(lhs(UV[2]))) = K;
UVrule := eval((rhs=lhs)~(UV), Krule);
mean_dp_with_q := eval(mean_dp, UVrule);

     
# It is obvious that q__ER defined by solve(mean_dp_with_q=0, q) cannot be computed 
# for f__D has no algebraic expression.
#
# So let's assume q__ER is known from any other mean


# Proceed as previously to compute Mean_p = Mean(Pi__ER):

pUV := eval(Pi__ER(beta__R), UV):
pUV := piecewise(_R <= U, op(2, pUV), _R < V, op(4, pUV), op(6, pUV)):

mean_p := Mean(pUV) assuming U < V:
mean_p := eval(%, int=Int):
mean_p := add(rhs~(map(Rule[`c*`], [op(mean_p)])));

 

mean_p__ER := eval(mean_p, q=q__ER);

ABCDE.mw

I've just added a few "partially instanciated" examples at the end of the attached file to (try to) convince you that there is no hope in finding a formal expression  of q__ER and Pi__ER  until f__D and g__beta are abstrat

 

Hi,

As long as you have not defined the PDF or CDF of D (that you must declare as local as D is a protected name) you will not be able to find the expression of q.

In the attached file you will see how far you can go FORMALLY with an abstract random variable D (that is a random variable for wich the PDF is just f__D(t) [I assume  D is continuous and f__D exists] ).
At the end you will find 4 examples where solve (..., q) doesn't fail.
 

restart:
with(Statistics):
local D:
d := Distribution(PDF = (t -> f__D(t)), Quantile = (z -> int(f__D(t), t=-infinity..q)));
D := RandomVariable(d);

TC := piecewise(q-D >= 0, c__o*(q-D), q-D < 0, c__u*(-q+D));

M := Mean(TC);
diff(%, q):
expand(%);
eval(%, op([2, -1], %)=1-op([1,-1], %));
IntegralTerm := op([1,-1], %);
eval(%%, IntegralTerm = J);
isolate(%, J);
expr := eval(%, J=IntegralTerm);


'Quantile(Q, q)' = Quantile(D, q);

# here is the rule to rewrite the lhs of expr as the value of the CDF of D at quantile q

MyRule := q -> subs(t=lhs(op(-1, lhs(expr))), Quantile(D, q)) = F__D(q):

Expr := eval(expr, MyRule(q));

solve(Expr, q);


# Example 1:

eval(Expr, F__D(q) = Quantile(Uniform(0, 1), q));
solve(%, q);
                             c__u    
                          -----------
                          c__o + c__u


# Example 2:

eval(Expr, F__D(q) = Quantile(Exponential(lambda), q));
solve(%, q);
                    /          c__u        \    
                -exp|- --------------------| + 1
                    \  lambda (c__o + c__u)/    


# Example 3:

eval(Expr, F__D(q) = Quantile(BetaDistribution(a, b), q));
solve(%, q);


                a                                            
   /   c__u    \           /                        c__u    \
   |-----------|  hypergeom|[a, 1 - b], [1 + a], -----------|
   \c__o + c__u/           \                     c__o + c__u/
   ----------------------------------------------------------
                          Beta(a, b) a                       


# Example 4:

eval(Expr, F__D(q) = Quantile(Normal(mu, sigma), q));
solve(%, q);

                /                            (1/2)\    
           1    |(c__o mu + c__u mu - c__u) 2     |   1
         - - erf|---------------------------------| + -
           2    \      2 sigma (c__o + c__u)      /   2

Solve_is_impossible.mw
 

A few errors :

• plot
  • the stntax is not correct: a range is needed for t
  • M is still undefined: you need to set it to a numerical value before plotting ; for instance

    plot(eval( dyr(t), M = 1), t=0..1);
    

• dyr(5) cannot work because the expression of dyr(t) contains derivatives with respect to t: what would it mean to "derive with respect to 5"?

  eval(dyr(t), t=5)
            /               (1/4)                    (3/4)
          M \-0.6609839668 4      + 0.0005882438542 4     
  
                                (1/2)\
             - 0.0001114078453 4     /
  
  
  

• This is basically the sme problem for z(1) ; do this instead

  z := unapply(diff(dyr(t), t), t):
  z(5)
           /                (1/4)                    (3/4)
         M \-0.07155390824 4      + 0.0001061325461 4     
  
                                (1/2)\
            - 0.00002412060546 4     /
  

   

dyr.mw