Advanced Search

Rouben Rostamian

7711 Reputation

22 Badges

19 years, 3 days
Contact Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

Animation...

Rouben Rostamian 7711
April 04 2018
1 3

Worksheet: cone-map.mw

Suggestion...

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

To get you started......

Rouben Rostamian 7711
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 7711
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 7711
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 7711
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);

 

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