Advanced Search

Rouben Rostamian

8531 Reputation

23 Badges

19 years, 75 days
Contact Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

Not a phase portrait, but......

Rouben Rostamian 8531
April 13 2018
1 0

Asking for a phase portrait for your equation is not exactly meaningful.  The responses that you have received present graphs of solutions y(x) versus x, which is not the usual meaning of a phase portrait.  A phase portrait is a plot in the phase space.

The concept of the phase space, however, pertains to autonomous differential equations.  Your equation is non-autonomous, so the graphs that you have been presented with are remotely reminiscent of phase portraits but you should be aware that they are not phase portraits.

That said, here is what I would do if I were to plot a set of graphs of y(x) versus x (not a phase portrait).

> 

restart;
> 

with(DEtools):
> 

de := diff(y(x),x)=x^2-y(x)^2;

diff(y(x), x) = x^2-y(x)^2

> 

a := 2.0;  b := 3.0;

2.0

3.0

> 

ic := seq(y(s)=+a, s=-a..a, 0.3),
      seq(y(s)=-a, s=-a..a, 0.3),
      seq(y(-a)=s, s=-a..a, 0.3),
      seq(y(+a)=s, s=-a..a, 0.3):
> 

DEplot(de, y(x), x=-b..b, [ic], y=-b..b,
    linecolor=black, thickness=1, arrows=none);


 

Download mw.mw

Answer...

Rouben Rostamian 8531
April 10 2018
3 5
> 

restart;
> 

with(LinearAlgebra):
> 

A := <
1, 0, 0, 1;
0, -1, 1, 0;
0, 1, -1, 0;
1, 0, 0, 1 >;

_rtable[18446884057743726822]

Here are A's eigenvalues and eigenvectors:

> 

evals, U := Eigenvectors(A);

evals, U := Vector(4, {(1) = -2, (2) = 2, (3) = 0, (4) = 0}), Matrix(4, 4, {(1, 1) = 0, (1, 2) = 1, (1, 3) = -1, (1, 4) = 0, (2, 1) = -1, (2, 2) = 0, (2, 3) = 0, (2, 4) = 1, (3, 1) = 1, (3, 2) = 0, (3, 3) = 0, (3, 4) = 1, (4, 1) = 0, (4, 2) = 1, (4, 3) = 1, (4, 4) = 0})

Make a diagonal matrix of the eigenvalues:

> 

Lambda := DiagonalMatrix(evals);

_rtable[18446884057743703326]

We have A = U.Lambda.(1/U).  Let's verify that:

> 

A - U . Lambda . U^(-1);

_rtable[18446884057743692134]

Here are two equivalent ways of calculating e^A.

 

Method 1:

> 

U . MatrixExponential(Lambda) . U^(-1);

_rtable[18446884057743684062]

Method 2 (the preferred way) does not require explicitly calculating eigenvalues and eigenvectors:

> 

MatrixExponential(A);

_rtable[18446884057743676950]

Download mw.mw

Answer...

Rouben Rostamian 8531
April 10 2018
1 0

Define φ and ψψ as you already have.

Then:

r := 2^k * M;
phi_vec := Vector[column](r, i-> phi[i]);
`&psi;&psi;_vec` := Vector[row](r, j-> `&psi;&psi;`[j]);
P := phi_vec . `&psi;&psi;_vec`;

 

Aside: The symbol `&psi;&psi;` in your calculations is too cumbersome and unnecessary. Why not use a single-letter symbol such as an uppercase Ψ insteaad?

Assembling matrices...

Rouben Rostamian 8531
April 08 2018
2 2

Here is a worksheet that shows the general method of assembling a matrix from smaller matrices.

I must alert you, however, to your use of M^(−1) *~ K in your worksheet. That certainly IS NOT the same as the matrix product M^(−1) . K. Which of those constructions do you have in mind?

> 

restart;
> 

with(LinearAlgebra):
> 

r := 2;

2

Putting fiour rxr arbitrary matrices together:

> 

< Matrix(r, symbol=a), Matrix(r, symbol=b);
  Matrix(r, symbol=c), Matrix(r, symbol=d) >;

_rtable[18446884266110072342]

Putting fiour rxr matrices together, the first two being the zero matrix and the identity matrix:

> 

< ZeroMatrix(r), IdentityMatrix(r);
  Matrix(r, symbol=c), Matrix(r, symbol=d) >;

_rtable[18446884266110064990]

Download mw.mw

freeze/thaw...

Rouben Rostamian 8531
April 08 2018
3 1
> 

restart;
> 

z := sin(2*x+3*y);

sin(2*x+3*y)

> 

z_tmp := subs(2*x=freeze(2*x), 3*y=freeze(3*y), z);

sin(`freeze/R0`+`freeze/R1`)

> 

thaw(expand(z_tmp));

sin(2*x)*cos(3*y)+cos(2*x)*sin(3*y)

Animation...

Rouben Rostamian 8531
April 04 2018
1 3

Worksheet: cone-map.mw

Suggestion...

Rouben Rostamian 8531
March 31 2018
1 0
> 

restart;

Note: I advise against assigning values to parameters because it messes up the

worksheet. Instead, define the parameter values without assigning them,

like this:

> 

params := beta=1, alpha[1]=.5;

beta = 1, alpha[1] = .5

> 

f := beta+alpha[1]+sin(beta*x);

beta+alpha[1]+sin(beta*x)

> 

plot(eval(f,[params]), x=0..Pi,
    labels=[x, convert("f'", symbol)(x)]):
plots:-textplot([Pi, 2.5, typeset(beta=1, ", ", alpha[1]=0.5, "\n\n")],
    align=[left], font=["Arial", 10, Bold]):
plots:-display([%,%%], axes=boxed, size = [300, 270]);

 

Download mm.mw

conic sections demo...

Rouben Rostamian 8531
March 28 2018
3 1

In this worksheet I call Maple's Explore() to visualize conic sections.  The graphics produced by Explore() cannot be shown on this web page.  Download the worksheet and load into Maple to play with it.

> 

restart;
> 

alpha := Pi/6;  # cone's vertex's half-angle

(1/6)*Pi

> 

doit := proc(beta)
  plots:-display([
    plot3d([ s*cos(t)*sin(alpha), s*sin(t)*sin(alpha), s*cos(alpha)],
        s=-5..5, t=0..2*Pi, color=gold),
    plot3d([2*sin(alpha) + t*cos(beta), s, -2*cos(alpha) + t*sin(beta)],
        s=-3..3, t=-4..6, color=red)
    ], style=surface, scaling=constrained, axes=none
     , orientation=[115, 80, 0]
  );
end proc:
> 

Explore(doit(beta),
  parameters = [beta=0..Pi/2],
  initvalues = [beta=Pi/4]
);
 

 

Download mw.mw

Addendum: The graphics may be improved significantly by adding a view option next to the orientation, as in:
    orientation=[115, 80, 0], view=[-3..6, -3..3, -4..4]

 

add versus sum...

Rouben Rostamian 8531
March 28 2018
1 0

Maple's add() and sum() have different evaluation rules. In general, you should use add(), not sum(). if the summation range is known ahead of the time.  And if you do that, you won't need the quotation marks around a__||k.  Here is the fixed code:
restart:
p := x -> add(a__||k * x^k, k=0..5):
p(x);
p(1);
p(0);

Perhaps this?...

Rouben Rostamian 8531
March 27 2018
2 7

I can't tell the purpose of the {y=-20}.  What is it supposed to do?  I am going to drop that and then fix the rest of the code.  We get:

restart;
with(PolyhedralSets):
p := PolyhedralSet([-x >= 0, -y >= 0, -z >= 0, -(1/2)*x >= 0, (-x-y)*(1/2) >= 0, (-x-z)*(1/2) >= 0]);
Plot(p);  # note the uppercase P

Is this what you want?...

Rouben Rostamian 8531
March 26 2018
1 3
> 

restart;

Let f be our Fourier series:

> 

f := Sum(A(x,t,n), n=1..infinity);

Sum(A(x, t, n), n = 1 .. infinity)

Define F

> 

F := (X,T,N) -> eval(f, [x=X, t=T, infinity=N]);

proc (X, T, N) options operator, arrow; eval(f, [x = X, t = T, infinity = N]) end proc

Then we will have:

> 

F(x,t,8);

Sum(A(x, t, n), n = 1 .. 8)

Download mw.mw

Answer...

Rouben Rostamian 8531
March 24 2018
2 14
> 

restart;
> 

G := (x,u) -> piecewise(x<u, 1/6*x^2*(1-u)^2*(3*u-2*u*x-x),
                             1/6*u^2*(1-x)^2*(3*x-2*x*u-u));

G := proc (x, u) options operator, arrow; piecewise(x < u, (1/6)*x^2*(1-u)^2*(-2*u*x+3*u-x), (1/6)*u^2*(1-x)^2*(-2*u*x-u+3*x)) end proc

> 

int(G(x,u)*f(u), u=0..1) assuming x>0, x<1;

int(-(1/6)*f(u)*u^2*(-1+x)^2*(2*u*x+u-3*x), u = 0 .. x)+int(-(1/6)*f(u)*x^2*(-1+u)^2*(2*u*x-3*u+x), u = x .. 1)

Download mw.mw

Series expansion...

Rouben Rostamian 8531
March 23 2018
3 6

I am guessing that you are looking for a series expansion of the solution in terms of powers of δ.  If so, then this worksheet can help.

We obtain a series expansion in terms of a small parameter, of the solution of a second

order nonlinear differential equation.

> 

restart;
> 

de := diff(z(x),x,x) - z(x) - cos(2*x)/(1+delta*z(x)) = 0;

diff(diff(z(x), x), x)-z(x)-cos(2*x)/(1+delta*z(x)) = 0

> 

 bc := z(-Pi/4)=0, z(Pi/4)=0

z(-(1/4)*Pi) = 0, z((1/4)*Pi) = 0

We are going to expand the solution of this nonlinear boundary value problem in powers of delta,

and pick the leading N terms like this:

> 

N := 2;  # change as needed
S := z(x) = add(u[n](x)*delta^n, n=0..N-1);

2

z(x) = u[0](x)+u[1](x)*delta

Substitute the series expression S in the differential equation and simplify:

> 

eval(de, S):
series(lhs(%), delta, N):
collect(%, delata, simplify):
tmp := convert(%, polynom);

diff(diff(u[0](x), x), x)-u[0](x)-cos(2*x)+(diff(diff(u[1](x), x), x)-u[1](x)+cos(2*x)*u[0](x))*delta

This is supposed to be zero for all delta, therefore the coefficient of each power of delta is zero.
This gives us a system of N coupled differential equations in N unknowns:

> 

DEs := coeffs(tmp, delta);

diff(diff(u[0](x), x), x)-u[0](x)-cos(2*x), diff(diff(u[1](x), x), x)-u[1](x)+cos(2*x)*u[0](x)

Apply the given boundary conditions to generate boundary conditions for the new unknowns:

> 

BCs := seq(subs(z=u[n], [bc])[], n=0..N-1);

u[0](-(1/4)*Pi) = 0, u[0]((1/4)*Pi) = 0, u[1](-(1/4)*Pi) = 0, u[1]((1/4)*Pi) = 0

> 

sol := dsolve({DEs, BCs});

{u[0](x) = -(1/5)*cos(2*x), u[1](x) = -1/10-(1/170)*cos(4*x)+(4/85)*exp(-x)/cosh((1/4)*Pi)+(4/85)*exp(x)/cosh((1/4)*Pi)}

Now build the final solution:

> 

eval(S, sol):
sort(%, delta, ascending);

z(x) = -(1/5)*cos(2*x)+(-1/10-(1/170)*cos(4*x)+(4/85)*exp(-x)/cosh((1/4)*Pi)+(4/85)*exp(x)/cosh((1/4)*Pi))*delta

Remark: There is no impediment, in principle, to calculating higher order terms by increasing N.

Trying N = 3, however, reveals a Maple bug.  The solution fails at the dsolve() step.  We may

solve the system manually without much difficulty, because the coupling of the system's equations

is in triangular form and the solution may be achieved by solving the equations in sequence, a single

equation at a time.

Download mw.mw

PS: The graphs in tomleslie's post correspond to the expression obtained above for z(x), with various choices of δ.

Analytical solution...

Rouben Rostamian 8531
March 19 2018
3 3
> 

restart;
> 

pde:=diff(u(x,t),t$2)+diff(u(x,t),x$4)=0;

diff(diff(u(x, t), t), t)+diff(diff(diff(diff(u(x, t), x), x), x), x) = 0

> 

bc:=u(0,t)=-12*t^2,u(1,t)=1-12*t^2,D[1,1](u)(0,t)=0,D[1,1](u)(1,t)=12;

u(0, t) = -12*t^2, u(1, t) = -12*t^2+1, (D[1, 1](u))(0, t) = 0, (D[1, 1](u))(1, t) = 12

> 

ic:=u(x,0)=x^4,D[2](u)(x,0)=0;

u(x, 0) = x^4, (D[2](u))(x, 0) = 0

> 

sol:=pdsolve({pde,ic,bc},u(x,t),HINT=`+`);

u(x, t) = x^4-12*t^2

Download mw.mw

How about this?...

Rouben Rostamian 8531
March 18 2018
2 0

Since your purpose is to examine individual terms, perhaps you would prefer to do this:

> 

s := Sum(a*i, i=1..n);

Sum(a*i, i = 1 .. n)

> 

op(1,s) $i=1..3;

a, 2*a, 3*a
First 35 36 37 38 39 40 41 Last Page 37 of 58