Advanced Search

Rouben Rostamian

7711 Reputation

22 Badges

19 years, 39 days
Contact Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

Figures...

Rouben Rostamian 7711
March 14 2018
1 7

Save the figure in the Encapsulated PostScript (EPS) format.  Then in latex do:

\includegraphics[width=0.4\textwidth]{myfigure.eps}

 

Answer...

Rouben Rostamian 7711
March 14 2018
3 6

At the fixed end the displacement and slope are zero.  At the free end the moment and shear are zero.

So, if the beam goes from x=0 to x=L, we have:

u(0)=0,  u'(0)=0,  u''(L)=0,  u'''(L)=0.

 

Will this do?...

Rouben Rostamian 7711
March 13 2018
0 0

Worksheet: mw.mw

RootOf...

Rouben Rostamian 7711
March 12 2018
0 13

The _Z that appears in RootOf is not a global variable at all.  On the contrary, it is a local variable inside the RootOf.

Let's see.  RootOf(x^2−4) evaluates to +2, −2.  Not surprisingly, RootOf(y^2−4) also evaluates to +2, −2.  What I am saying is that in a RootOf expression the name of the variable is immaterial., that is, it does not affect the expression's value.  Such a variable is commonly known as a dummy variable.  Maple uses the generic dummy variable name _Z in the argument of RootOf to distinguish that variable from anything else that you may have in your worksheet.  That's all.

Exact solution...

Rouben Rostamian 7711
March 09 2018
1 0

Maple can solve your equation explicitly:

> 

restart;
> 

eq := eta*(2 + 9*exp(-eta*y)) - 22*eta*y = 0;

eta*(2+9*exp(-eta*y))-22*eta*y = 0

> 

sol := isolate(eq, y);

y = (9/22)*exp(-LambertW((9/22)*eta*exp(-(1/11)*eta))-(1/11)*eta)+1/11

> 

plot(rhs(sol), eta=0..1, view=0..0.6, color=red);

Download mw.mw

Solution...

Rouben Rostamian 7711
March 08 2018
0 0
> 

restart;
> 

with(LinearAlgebra):

Proc to generate the nth order Caley Omega differential operator.

> 

Caley_Omega := proc(n::posint)
  local L, q, x, term_to_D;
  Matrix(n,n, [x[i] $i=1..n^2]);
  Determinant(%);
  L := [op(1..-1, %)];
  term_to_D := proc(q)
    `if`(nops(q)=n,1,-1)*D[seq(op([-i,1],q), i=1..n)];
  end proc;
  add(term_to_D(q), q in L);
end proc:
> 

 

Illustration: Select n:

> 

n := 2;

2

Here is the Caley Omega operator:

> 

C := Caley_Omega(n);

D[4, 1]-D[3, 2]

Let's call the variables x__ij, i, j = 1 .. n:

> 

X := seq(seq(x[i,j], j=1..n), i=1..n);

x[1, 1], x[1, 2], x[2, 1], x[2, 2]

Apply C to F(X)

> 

C(F)(X);

(D[1, 4](F))(x[1, 1], x[1, 2], x[2, 1], x[2, 2])-(D[2, 3](F))(x[1, 1], x[1, 2], x[2, 1], x[2, 2])

Another example with variables named x_i, i = 1 .. n^2

> 

X := x[i] $i=1..n^2;

x[1], x[2], x[3], x[4]

> 

C(F)(X);

(D[1, 4](F))(x[1], x[2], x[3], x[4])-(D[2, 3](F))(x[1], x[2], x[3], x[4])

> 

 

Larger n:

> 

n := 3;

3

> 

C := Caley_Omega(n);

D[9, 5, 1]-D[8, 6, 1]-D[9, 4, 2]+D[7, 6, 2]+D[8, 4, 3]-D[7, 5, 3]

> 

X := x[i] $i=1..n^2;

x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9]

> 

C(F)(X);

(D[1, 5, 9](F))(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9])-(D[1, 6, 8](F))(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9])-(D[2, 4, 9](F))(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9])+(D[2, 6, 7](F))(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9])+(D[3, 4, 8](F))(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9])-(D[3, 5, 7](F))(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9])

 

 

Download mw2.mw

Answer...

Rouben Rostamian 7711
March 05 2018
1 0

I see that you still do unprotect(Pi) despite multiple warnings that you have received.

From now on you take the responsibility for that.

Regarding the differentiation problem, here is what to do.

> 

restart;
> 

F := x^2 + x*y + y^2;

x^2+x*y+y^2

> 

x := 10;
y := 15;

10

15

> 

diff(F, x);

Error, invalid input: diff received 10, which is not valid for its 2nd argument

That fails because x has been set to 10, and therefore diff(F, x)  is really diff(F, 10)

which makes no sense.

 

Here is how to fix it.

> 

restart;
> 

F := x^2 + x*y + y^2;

x^2+x*y+y^2

> 

params := { x=10, y=15 };

{x = 10, y = 15}

> 

diff(F, x);
eval(%, params);

2*x+y

35

Answer...

Rouben Rostamian 7711
March 04 2018
0 0

I don't know what you mean by "normalised forward sensitivity index formula".  

Let me tell you what I mean by sensitivity and you decide if that's what you want.

 

If we have a function of three variable, let's say F(x, y, z), then the differential of F is
"dF = (∂F)/(∂x)dx + (∂F)/(∂y)dy+(∂F)/(∂z)dz".
We rearrange this into
"dF/(F)=x/(F)*(∂F)/(∂x)*dx/(x)+y/(F)*(∂F)/(∂y)*(dy)/(y)+z/(F)*(∂F)/(∂z)*(dz)/(z)".

The fractions dx/x, dy/y, dz/z measure the relative changes in x, y, z.

The fraction dF/F measures the corresponding relative change in "F."

The coefficient "(x/(F)*(∂F)/(∂x))" of dx/xmeasures the sensitivity of F relative to "x."

If that's what you understand by sensitivity, then read on.  Otherwise explain

what you mean.

 

As to your code, you have been told already to never unprotect Pi. If you

want the Greek letter Pi, just write pi with lowercase p.

 

Additionally, you may have noticed that square brackets [ ]  have special meaning

in Maple.  The square brackets in the definition of your R__0 should be round

parentheses. Here is the fixed version of you R__0:

> 

restart;
> 

R[0] := k*tau*(rho*(Upsilon*(mu+alpha+eta)+chi)+(1-rho)*(Upsilon*(1-q)*eta+mu+beta+chi))*(pi*(-mu*p+mu+phi)/(mu*(vartheta+mu+phi))+epsilon*pi*(mu*p+vartheta)/(mu*(vartheta+mu+phi)))/((mu+beta+chi)*(mu+alpha+eta)-chi*eta*(1-q));

k*tau*(rho*(Upsilon*(mu+alpha+eta)+chi)+(1-rho)*(Upsilon*(1-q)*eta+mu+beta+chi))*(pi*(-mu*p+mu+phi)/(mu*(vartheta+mu+phi))+epsilon*pi*(mu*p+vartheta)/(mu*(vartheta+mu+phi)))/((mu+beta+chi)*(mu+alpha+eta)-chi*eta*(1-q))

From what we learned in this article's introduction, to determine the sensitivity

of R__0 relative to alpha, we do

> 

factor(alpha/R[0]*diff(R[0],alpha));

-(beta*rho+mu*rho-beta-chi-mu)*(Upsilon*eta*q-Upsilon*eta-beta-chi-mu)*alpha/((chi*eta*q+alpha*beta+alpha*chi+alpha*mu+beta*eta+beta*mu+chi*mu+eta*mu+mu^2)*(Upsilon*eta*q*rho+Upsilon*alpha*rho-Upsilon*eta*q+Upsilon*mu*rho+Upsilon*eta-beta*rho-mu*rho+beta+chi+mu))

> 

 

You may determine the sensitivity of R__0 relative to the other 14 variables in the same way.

Download mw2.mw

Solution...

Rouben Rostamian 7711
March 04 2018
0 4
> 

restart;

I will do the calculations with rmax = π.  Modify as needed.

 

I expand the complex basis functions c[k]*e^(I*k*x) into (a[k]+I*b[k])*(cos(k*x)+I*sin(k*x)) 

through Euler's formula, where a[k] and b[k] are real.  The index k goes from -N to "N."

> 

N := 2;

2

Here is the function f represented in the selected basis:

> 

f := add((a[k]+I*b[k])*(cos(k*x) + I*sin(k*x)), k=-N..N);

(a[-2]+I*b[-2])*(cos(2*x)-I*sin(2*x))+(a[-1]+I*b[-1])*(cos(x)-I*sin(x))+a[0]+I*b[0]+(a[1]+I*b[1])*(cos(x)+I*sin(x))+(a[2]+I*b[2])*(cos(2*x)+I*sin(2*x))

Evaluate abs(f(x))^2:

> 

conjugate(f)*f assuming real:
f2x := simplify(%):

Evaluate abs(f(y))^2:

> 

f2y := subs(x=y, f2x):

Evaluate  "|f '(x)|^(2)":

> 

diff(f,x):
conjugate(%)*% assuming real:
fp2 := simplify(%):

The constraint (change as needed):

> 

constr := add(a[k]^2 + b[k]^2, k=-N..N) = 1;

a[-2]^2+a[-1]^2+a[0]^2+a[1]^2+a[2]^2+b[-2]^2+b[-1]^2+b[0]^2+b[1]^2+b[2]^2 = 1

A simple choice for the function w

> 

w := s -> s^2;

proc (s) options operator, arrow; s^2 end proc

Here is the objective function.  Adjust as desired:

> 

E := int(fp2, x=0..Pi)
   + int(w(x-y)*f2x*f2y, x=0..Pi, y=0..Pi):
> 

infolevel[Optimization] := 2:

This will take a few minutes:

> 

ans := Optimization:-NLPSolve(E, {constr}, iterationlimit=100);

NLPSolve: calling NLP solver

NLPSolve: using method=sqp
NLPSolve: number of problem variables 10
NLPSolve: number of nonlinear inequality constraints 0
NLPSolve: number of nonlinear equality constraints 1
NLPSolve: number of general linear constraints 0
NLPSolve: trying evalhf mode
attemptsolution: number of major iterations taken 51

[.358565694712160155, [a[-2] = HFloat(0.033202595604861884), a[-1] = HFloat(0.42176063010072834), a[0] = HFloat(-0.12728560464292304), a[1] = HFloat(-0.4217606642506425), a[2] = HFloat(0.033202596782105666), b[-2] = HFloat(-0.1919507392055417), b[-1] = HFloat(0.07295391355332227), b[0] = HFloat(0.7358632696288413), b[1] = HFloat(-0.07295383471334241), b[2] = HFloat(-0.19195074805401965)]]

Here is the minimizing function f

> 

eval(f, ans[2]):
Re(%) + I*Im(%) assuming x::real;

HFloat(0.06640519238696754)*cos(2*x)+HFloat(8.848477950351707e-9)*sin(2*x)-HFloat(3.414991417427515e-8)*cos(x)+HFloat(0.1459077482666647)*sin(x)-HFloat(0.12728560464292304)+I*(HFloat(1.1772437816248704e-9)*sin(2*x)-HFloat(0.38390148725956136)*cos(2*x)-HFloat(0.8435212943513708)*sin(x)+HFloat(7.883997986402047e-8)*cos(x)+HFloat(0.7358632696288413))

The graphs of the real and imaginary parts of f:

> 

eval(f, ans[2]):
[Re(%), Im(%)] assuming x::real:
plot(%, x=0..Pi, color=[red,blue]);

Download mw.mw

Variation...

Rouben Rostamian 7711
March 02 2018
0 0

To add to what Kitonum has already said:

restart;
with(plots): with(plottools):
f := (x,y)->sin(x)*cos(y):
p1 := plot3d(f, 0..2*Pi, 0..2*Pi, style=patchcontour, shading=zhue):
p2 := contourplot(f, 0..2*Pi, 0..2*Pi, filledregions=true):
p3 := transform((x,y)->[x,y,-1.5])(p2):
display([p1,p3]);

 

Solution...

Rouben Rostamian 7711
February 28 2018
1 0

@amcmaple

This is likely too complex as a first exercise for plotting phase portraits, but here it is,

just to show you the sort of stuff that you need to learn.

> 

restart;
> 

with(DEtools):

What you have is Van der Pol's differential equation:

> 

de := diff(u(t),t,t) - epsilon*(1-64*u(t)^2)*diff(u(t),t) + 16*u(t) = 0;

diff(diff(u(t), t), t)-epsilon*(1-64*u(t)^2)*(diff(u(t), t))+16*u(t) = 0

For doing the phase portrait, it is convenient to convert this second order differential equation to an
equivalent system of two first order equations through introducing the auxiliary variable v = du/dt

> 

de1 := diff(u(t),t) = v(t);
de2 := diff(v(t),t) = epsilon*(1-64*u(t)^2)*v(t) - 16*u(t);

diff(u(t), t) = v(t)

diff(v(t), t) = epsilon*(1-64*u(t)^2)*v(t)-16*u(t)

We select a good choice of initial conditions:

> 

ic := seq([u(0)=0, v(0) = s], s=-2..2, 0.4);

[u(0) = 0, v(0) = -2], [u(0) = 0, v(0) = -1.6], [u(0) = 0, v(0) = -1.2], [u(0) = 0, v(0) = -.8], [u(0) = 0, v(0) = -.4], [u(0) = 0, v(0) = 0.], [u(0) = 0, v(0) = .4], [u(0) = 0, v(0) = .8], [u(0) = 0, v(0) = 1.2], [u(0) = 0, v(0) = 1.6], [u(0) = 0, v(0) = 2.0]

Here is my choice of epsilon:

> 

epsilon := 2.0;

2.0

and here the phase portrait:

> 

DEplot([de1,de2], [u(t),v(t)], t=0..3, [ic],
  linecolor=black, thickness=1, stepsize=0.02,
  u=-0.5..0.5, v=-2.5..2.5, labels=[u(t),v(t)]);

Download mw.mw

Linearization...

Rouben Rostamian 7711
February 25 2018
3 1

Your question is somewhat ambiguous.  You have quite a few variables in your expression.

You want to linearize with respect to which one?  I suspect that you have the variable psi1(t)

in mind, and that's what I will pick.  Additionally, "linearize" by itself is not well-defined"-"you

need say linearize about something!  Since you haven't specified that, I will

linearize about an arbitrary (but known) function `ψ__0`(t).  So here we go.

> 

restart;

The expression

> 

expr[0] := sin(phi2)*sin(psi1(t))*lk*m1*(diff(psi1(t), t, t))+sin(phi2)*(diff(psi1(t), t))^2*cos(psi1(t))*lk*m1;

sin(phi2)*sin(psi1(t))*lk*m1*(diff(diff(psi1(t), t), t))+sin(phi2)*(diff(psi1(t), t))^2*cos(psi1(t))*lk*m1

may be factorized into

> 

expr[1] := factor(expr[0]);

sin(phi2)*lk*m1*((diff(psi1(t), t))^2*cos(psi1(t))+sin(psi1(t))*(diff(diff(psi1(t), t), t)))

For linearization purposes the factor sin(phi2)*lk*m1is immaterial, therefore we drop it for now:

> 

expr[2] := op(-1, expr[1]);

(diff(psi1(t), t))^2*cos(psi1(t))+sin(psi1(t))*(diff(diff(psi1(t), t), t))

Through visual inspection we recognize this to be

> 

expr[3] := -Diff(cos(psi1(t)), t,t);

-(Diff(cos(psi1(t)), t, t))

Let's verify that claim:

> 

simplify(expr[3] - expr[2]);

0

OK then, your question comes down to: How do we linearize expr[3]?

 

Let's says we wish to linearize the expression about some known function `ψ__0`(t), that is,

we are considering psi1(t) near `ψ__0`(t).  We express this mathematically by letting

> 

ansatz := psi1(t) = psi__0(t) + epsilon*u(t);

psi1(t) = psi__0(t)+epsilon*u(t)

where `ε` is a small constant (independent of "t)" and u(t) is the linearization variable.

This is like the first two terms in a Taylor series.  Plug the ansatz into expr[3] and then

expand into a series in `ε`

> 

eval(expr[3], ansatz);
series(%, epsilon, 2);

-(Diff(cos(psi__0(t)+epsilon*u(t)), t, t))

series(-(Diff(Diff(cos(psi__0(t)), t), t))-(Diff(Diff(-sin(psi__0(t))*u(t), t), t))*epsilon+O(epsilon^2),epsilon,2)

The coefficient of `ε` in that expansion is the desired linearlization.

 

We conclude that the linearization of expr[0] about a given function `ψ__0`(t) is

> 

sin(phi2)*lk*m1*Diff(sin(psi__0(t))*u(t),t,t);

sin(phi2)*lk*m1*(Diff(sin(psi__0(t))*u(t), t, t))

where u(t) is the linearization variable.

 

 

Download mw.mw

Exact solution...

Rouben Rostamian 7711
February 24 2018
1 0

Applying Lagrange multipliers and numerical methods to this problem is overkill.  The problem is simple enough to solve exactly by Maple, or even by hand, if you want.

> 

restart;

We wish to find the extrema of the objective function:

> 

f := x*y + y*z;

x*y+y*z

subject to constraints:  x*y = 1,  y^2+z^2 = 1.

 

The constraint y^2+z^2 = 1 may be equivalently expressed as y = sin(t),  z = cos(t) in terms

of a parameter "t."  But then the constraint x*y = 1 takes the form x = 1/sin(t).  In short,

we express x, y, z in terms of t

> 

rels := {(x,y,z) =~ (1/sin(t), sin(t), cos(t))};

{x = 1/sin(t), y = sin(t), z = cos(t)}

Then the objective function reduces to a function of a single parameter t

> 

g := unapply(eval(f, rels), t);

proc (t) options operator, arrow; 1+sin(t)*cos(t) end proc

We set the derivative of g to zero to obtain its extrema:

> 

D[1](g)(t) = 0;
t_vals := solve(%, t);

cos(t)^2-sin(t)^2 = 0

(1/4)*Pi, (3/4)*Pi

To determine the types of the extrema, we plug the t_vals into the second derivative of g

> 

D[1,1](g)~([t_vals]);

[-2, 2]

We see that the second derivative is negative at t = (1/4)*Pi, therefore we have a local maximum there.

Similarly, t = 3*Pi*(1/4) corresponds to a local minimum.

Reverting to the originalx, y, z variables, we have:

> 

eval(rels, t=t_vals[1]);
'f(x,y,z)'=eval(f, %);  # maximum

{x = 2^(1/2), y = (1/2)*2^(1/2), z = (1/2)*2^(1/2)}

f(x, y, z) = 3/2

> 

eval(rels, t=t_vals[2]);
'f(x,y,z)'=eval(f, %);  # minimum

{x = 2^(1/2), y = (1/2)*2^(1/2), z = -(1/2)*2^(1/2)}

f(x, y, z) = 1/2

We may exhibit the max and min graphically by plotting g

> 

plot(g(t), t=0..Pi, view=0..2);

 

 

Download mw.mw

Answer...

Rouben Rostamian 7711
February 24 2018
0 0

Although you wrote F:=(x,y)->(x+y,2*x^2+y^2), I suspect that you would rather have F as a vector-valued function.  Additionally, since F is a function from R2 to R2, the proper name of the derivative is the Jacobian, not the gradient.  With these in mind, here is what we do:

> 

restart;
> 

sys := [x+y, 2*x+y^2];

[x+y, y^2+2*x]

> 

F := unapply(Vector(sys), [x,y]):
> 

J := unapply(VectorCalculus:-Jacobian(F(x,y), [x,y]), [x,y]):

Thus we have:

> 

F(x,y);

_rtable[18446884420085923838]

> 

J(x,y);

_rtable[18446884420085924678]

Idea...

Rouben Rostamian 7711
February 20 2018
0 0

Write a proc f which receives a guess for T and produces a new guess:

f := proc(T)
  local Tnew;
  blah blah ... compute the new T
  return Tnew;
end proc;

Then apply the proc in a loop, as in:

T1 := starting value of T;
T := f(T1);
while abs(T-T1)> 1e-3 do
  T1 := T;
  T := f(T1);
end do;

Upon exit from that loop, the value of T will be approximately what you expect.  The accuracy iis controlled by the value 1e-3 in the code above.  You may increase or decrease that as needed.

PS: I saw rlopez's note after I posted this.  And he is right; the way you have worded makes no sense.  In what I have shown I have assumed that we want two successive values of T to be the same, which agrees rlopez's guess.

PPS: In general there is no guarantee that such an iteration will converge.  You may want to add a counter which is incremented in each pass of the loop, and break out of the loop if the number of iterations exceeds a prescribed maximum.

First 32 33 34 35 36 37 38 Last Page 34 of 53