Advanced Search

Rouben Rostamian

8506 Reputation

23 Badges

19 years, 37 days
Contact Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

To get you started......

Rouben Rostamian 8506
March 17 2018
1 8
  • It appears that you expect an analytical solution to that PDE.  That's highly unlikely.  The best you can do is to specify the parameters Omega, Lambda, etc., then have maple compute a numerical approximation to the solution.
  • What is the cos(2*Pi*x) doing there?  What is x?
  • The boundary conditions should be entered as;
    bc := D[1](y)(0,t)=0,  D[1](y)(1,t)=0;
  • You need to supply an initial condition as well.  Something like:
    ic := u(x,0) = sin(Pi*x);    # change as needed
  • Once everything is set up, you will do:
    sol := pdsolve(pde, {ic, bc}, numeric);

 

Try this one...

Rouben Rostamian 8506
March 17 2018
0 1

Calculate Cf.  Then:

eval(Cf,x=q+2*h):
expand(%):
simplify(%,size):
simplify(%):
y[n+2]=collect(%, [y[n],f[n],f[n+1],f[n+3/2],f[n+2]], factor);

 

Analytic solution...

Rouben Rostamian 8506
March 16 2018
3 2

This is in response to your statement: "For me it is important to show that Maple can't give an analytic solution. In a perfect situation I need to understand why".

Maple is pretty good in solving problems but it has its limitations.  The current version (Maple 2017) seems to be unable to obtain an analytic solution to your PDE on its own, but it can do it with some manual intervention.  See if this helps.

The Wave Equation through Fourier Series

We calculate the Fourier series solution of a boundary value problem for the wave

equation with a nonzero forcing term.  Perhaps one day Maple will be able to do

this on its own without manual intervention.  Currently (Maple 2017) does it
with some human help.

> 

restart;
pde := diff(u(x,t),t,t) = diff(u(x,t),x,x) + 5*sin(3*x);
bc := u(0,t) = 0, u(Pi,t) = 0;
ic := u(x,0) = 0, (D[2](u))(x,0) = 1;

diff(diff(u(x, t), t), t) = diff(diff(u(x, t), x), x)+5*sin(3*x)

  

u(0, t) = 0, u(Pi, t) = 0

  

u(x, 0) = 0, (D[2](u))(x, 0) = 1

 (1)

Step 1: Solve the homogeneous equation, that is, the PDE without the forcing

term "5*sin(3*x),"together with the boundary conditions:

> 

pdsolve({diff(u(x,t),t,t) = diff(u(x,t),x,x), bc}, HINT=`*`);

u(x, t) = Sum(sin(n*x)*(_C1(n)*sin(n*t)+_C5(n)*cos(n*t)), n = 1 .. infinity)

 (2)

Step 2: Find a particular solution of {PDE, bc}.  Any solution would do.

An easy choice would be to look for a solution that does not depend on time:

> 

dsolve({rhs(pde), bc});

u(x, t) = (5/9)*sin(3*x)

 (3)

Step 3: Form a candidate for the final solution:

> 

U := rhs((2)) + rhs((3));

Sum(sin(n*x)*(_C1(n)*sin(n*t)+_C5(n)*cos(n*t)), n = 1 .. infinity)+(5/9)*sin(3*x)

 (4)

Step 4: The solution candidate Usatisfies the PDE and the boundary

conditions.  What remains is to select the coefficients _C1(n) and _C5(n)

so that the initial conditions are satisfied.

 

Step 4a: Let's look at the condition u(x,0)=0:

> 

eval(U, t=0) = 0;

Sum(sin(n*x)*_C5(n), n = 1 .. infinity)+(5/9)*sin(3*x) = 0

 (5)

This says that all __C5(n) are zero except for _C5(3) which is -5/9.

Thus, U reduces to:

> 

U := sum(sin(n*x)*_C1(n)*sin(n*t), n=1..infinity) -sin(3*x)*5/9*cos(3*t) + 5*sin(3*x)/9;

sum(sin(n*x)*_C1(n)*sin(n*t), n = 1 .. infinity)-(5/9)*sin(3*x)*cos(3*t)+(5/9)*sin(3*x)

 (6)

Step 4b: Let's look at the condition D[2](u)(x,0)=1:

> 

diff(U, t):
eval(%, t=0) = 1;

sum(sin(n*x)*_C1(n)*n, n = 1 .. infinity) = 1

 (7)

To determine the coefficients _C1(n), we multiply this by sin(m*x) and

integrate over the interval 0, Pi

> 

sum(Int(sin(m*x)*sin(n*x)*_C1(n)*n, x=0..Pi), n=1..infinity) = Int(sin(m*x), x=0..Pi);

sum(Int(sin(m*x)*sin(n*x)*_C1(n)*n, x = 0 .. Pi), n = 1 .. infinity) = Int(sin(m*x), x = 0 .. Pi)

 (8)

It is easy to verify that

> 

Int(sin(m*x)*sin(n*x), x=0..Pi) = piecewise(m=n, Pi/2, 0);

Int(sin(m*x)*sin(n*x), x = 0 .. Pi) = piecewise(m = n, (1/2)*Pi, 0)

 (9)

Therefore the infinite sum in (8) reduces to a  (1/2)*Pi*_C1(m)*m. The right-hand

side integrates to:

> 

int(sin(m*x), x=0..Pi) assuming m::posint;

-(-1+(-1)^m)/m

 (10)

We conclude that

> 

Pi/2*_C1(m)*m = (10);

(1/2)*Pi*_C1(m)*m = -(-1+(-1)^m)/m

 (11)
> 

isolate((11), _C1(m));

_C1(m) = -2*(-1+(-1)^m)/(m^2*Pi)

 (12)

or equivalently:

> 

subs(m=n, (12));

_C1(n) = -2*(-1+(-1)^n)/(n^2*Pi)

 (13)

Step 5:  We are done.  Here is the final form of the solution:

> 

u(x,t) = eval(U, (13));

u(x, t) = sum(-2*sin(n*x)*(-1+(-1)^n)*sin(n*t)/(n^2*Pi), n = 1 .. infinity)-(5/9)*sin(3*x)*cos(3*t)+(5/9)*sin(3*x)

 (14)

Let's verify that we haven't goofed anywhere:

> 

pdetest((14), pde);

0

 (15)

That's OK!   Now let's plot the solution:

> 

subs(infinity=30, sol);
plot3d(%, x=0..Pi, t=0..Pi);

sum(-2*sin(n*x)*(-1+(-1)^n)*sin(n*t)/(n^2*Pi), n = 1 .. 30)-(5/9)*sin(3*x)*cos(3*t)+(5/9)*sin(3*x)

  

 
 

 

Worksheet: Download mw.mw

 

No need for diff_table...

Rouben Rostamian 8506
March 16 2018
1 5

Tomleslie has already shown the adjustments needed to make your code work.  The complications in that approach are due to the use of the diff_table construct.  Do you really need to use diff_table?  The problem may be solved in a more transparent manner without it:

restart;
pde := diff(u(x,t),t,t) = diff(u(x,t),x,x) + 5*sin(3*x);
bc := u(0,t)=0,  u(Pi,t)=0;
ic := u(x,0)=0,  D[2](u)(x,0)=1;
sol := pdsolve(pde, {bc, ic}, numeric);
sol:-plot3d(x=0..1, t=0..1);

 

Figures...

Rouben Rostamian 8506
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 8506
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 8506
March 13 2018
0 0

Worksheet: mw.mw

RootOf...

Rouben Rostamian 8506
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 8506
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 8506
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 8506
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 8506
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 8506
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 8506
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 8506
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

First 36 37 38 39 40 41 42 Last Page 38 of 58