Advanced Search

Rouben Rostamian

8506 Reputation

23 Badges

19 years, 37 days
Contact Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

Answer...

Rouben Rostamian 8506
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 8506
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 8506
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 8506
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 8506
April 04 2018
1 3

Worksheet: cone-map.mw

Suggestion...

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

Solution...

Rouben Rostamian 8506
March 18 2018
3 22

Solving a fractional differential equation

 

This worksheet implements that algorithm described in the article
115006_listening_sample_task_-_matching_example_1_.pdf

 cited earlier.

 

I have coded the algorithm literally as found in that article without verifying its

details.  I cannot be certain of the correctness of the algorithm since the article

seems to contain some careless errors.  In particular, the figures included in that

article clearly disagree with the given data, and therefore they cannot be

reproduced in this worksheet.

> 

restart;

We solve the initial value problem of the system of fractional

differential equation:

"(d^(alpha)x(t))/(dt^(alpha))=F(t,x(t),y(t)), ( d^(alpha)y(t))/(dt^(alpha))=G(t,x(t),y(t)), x(`t__0`)=`x__0`, y(`t__0`)=`y__0`, "
on the interval t__0 < t and t < `t__1 `by dividing the interval into N equal subintervals.

The solution is returned as an array of length N+1 of triplets "(`t__i`, `x__i`, `y__i`),"
where t__i = (1-i/N)*t__0+i*t__1/N, i = 0, () .. (), N are the dividing points of

the subintervals, and x__i, y__i is the computed solution at "t=`t__i`."

This proc implements the algorithm:

> 

fdsolve := proc({alpha:=NULL, t0:=NULL,
                 t1:=NULL, x0:=NULL, y0:=NULL,
                 N:=NULL}, params)
    local t, h, h1, h2, c, b, x, y, L, n, l, X, Y, f, g;
    eval(F(t,x,y), params);
    f := unapply(%, [t,x,y]);
    eval(G(t,x,y), params);
    g := unapply(%, [t,x,y]);
    L := floor(1/alpha);
    h := (t1 - t0)/N;
    h1 := h^alpha/GAMMA(alpha+1);
    h2 := h^alpha/GAMMA(alpha+2);
    c := (i,n) ->
        `if`(i=0,
            (n-1)^(alpha+1) - (n-1-alpha)*n^alpha,
            (n-i+1)^(alpha+1) + (n-i-1)^(alpha+1) - 2*(n-i)^(alpha+1));
    b := (i,n) -> (n-i)^alpha - (n-1-i)^alpha;
    t := Array(0..N, i-> (1-i/N)*t0 + i/N*t1, datatype=float[8]);
    x[0], y[0] := x0, y0;
    for n from 0 to N-1 do
        X[0], Y[0] :=
            x[0] + h1*add(b(i,n+1)*f(t[i],x[i],y[i]), i=0..n),
            y[0] + h1*add(b(i,n+1)*g(t[i],x[i],y[i]), i=0..n);
        for l from 1 to L do
            X[l], Y[l] :=
                x[0] + h2*add(c(i,n+1)*f(t[i],x[i],y[i]), i=0..n)
                     + h2*f(t[n+1], X[l-1], X[l-1]),
                y[0] + h2*add(c(i,n+1)*g(t[i],x[i],y[i]), i=0..n)
                     + h2*g(t[n+1], X[l-1], X[l-1]);
        end do;
        x[n+1], y[n+1] := X[L], Y[L];
    end do;
    return Array(0..N, i -> [t[i], x[i], y[i]]);
end proc:
> 

 

Example (taken from the cited article)

> 

F := (t,x,y) -> r*x*(1-x/k) - beta*x*y/(a+x^2);
G := (t,x,y) -> mu*beta*x*y/(a+x^2) - d*x - eta*x*y;

proc (t, x, y) options operator, arrow; r*x*(1-x/k)-beta*x*y/(a+x^2) end proc

  

proc (t, x, y) options operator, arrow; mu*beta*x*y/(a+x^2)-d*x-eta*x*y end proc

 (1)
> 

params := { r=0.05, a=0.8, mu=0.8, d=0.24,
            eta=0.01, beta=0.6, k=1.6 };

{a = .8, beta = .6, d = .24, eta = 0.1e-1, k = 1.6, mu = .8, r = 0.5e-1}

 (2)
> 

T := 10.0;  # solve over t=0 to t=T

10.0

 (3)

We produce several solutions starting from various initial data:

> 

# solution starting at (2.0, 0.5)
sol := fdsolve(alpha=0.98, t0=0.0, t1=T, x0=2.0, y0=0.5, N=50, params):
p1 := plot([seq([sol[i][2], sol[i][3]], i=0..50)]):
> 

# solution starting at (2.0, 1.0)
sol := fdsolve(alpha=0.98, t0=0.0, t1=T, x0=2.0, y0=1.0, N=50, params):
p2 := plot([seq([sol[i][2], sol[i][3]], i=0..50)]):
> 

# solution starting at (2.0, 2.0)
sol := fdsolve(alpha=0.98, t0=0.0, t1=T, x0=2.0, y0=2.0, N=50, params):
p3 := plot([seq([sol[i][2], sol[i][3]], i=0..50)]):
> 

# solution starting at (2.0, 3.0)
sol := fdsolve(alpha=0.98, t0=0.0, t1=T, x0=2.0, y0=3.0, N=50, params):
p4 := plot([seq([sol[i][2], sol[i][3]], i=0..50)]):

Display the phase portrait:

> 

plots:-display([p1,p2,p3,p4], color=[red,green,blue,gold], style=line);

 
 

 

Download frac-dsolve.mw

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