Advanced Search

Rouben Rostamian

8636 Reputation

23 Badges

19 years, 290 days
Contact Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

The cosine term...

Rouben Rostamian 8636
August 28 2019
3 3

The presence of the cosine term is more of an effect of the geometry than physics, although I dislike putting a hard barrier between the two subjects.

The Pythagorean Theorem says that the square of the hypotenuse is the sum of the squares of the two other sides in a right triangle.  The generalization to an arbitrary triangle is known as the law of cosines.  I trust that you already know that but in the unlikely case, look it up in Wikipedia.

The law of cosines enters your calculations because the two vectors that you have shown are not perpendicular. The pearl's velocity is the sum of those two vectors.  Calculating the sum relies on the law of cosines.

> 

restart;
> 

Typesetting:-Settings(typesetdot):
> 

 

> 

# the horizontal vector labled a*(diff(phi(t), t))in the diagram
< a*diff(phi(t),t), 0 >:

# the slanted vector labeled a*(diff(theta(t), t)) in the diagram
a*diff(theta(t),t) * < -cos(theta(t)), sin(theta(t)) >:

# the sum of the two vectors above, which is the pearl's absolute velocity
v := %% + %;

Vector[column](%id = 18446884233678926542)

The pearl's kinetic energy

> 

simplify(1/2*m*v^+ . v);

(1/2)*m*a^2*((diff(theta(t), t))^2-2*(diff(phi(t), t))*(diff(theta(t), t))*cos(theta(t))+(diff(phi(t), t))^2)

Download mw.mw

Solution...

Rouben Rostamian 8636
August 25 2019
3 0

This works in Maple 2019.  It may require minor adjustments in older versions.

> 

restart;
> 

A := < 1, 2; 3, 4 >;

Matrix(2, 2, {(1, 1) = 1, (1, 2) = 2, (2, 1) = 3, (2, 2) = 4})

> 

X := < x[1](t), x[2](t) >;

Vector(2, {(1) = x[1](t), (2) = x[2](t)})

> 

sys := diff(X,t) =~ A . X;

Vector(2, {(1) = diff(x[1](t), t) = x[1](t)+2*x[2](t), (2) = diff(x[2](t), t) = 3*x[1](t)+4*x[2](t)})

> 

dsolve(sys);

{x[1](t) = _C1*exp((1/2)*(5+33^(1/2))*t)+_C2*exp(-(1/2)*(-5+33^(1/2))*t), x[2](t) = -(1/4)*_C2*exp(-(1/2)*(-5+33^(1/2))*t)*33^(1/2)+(1/4)*_C1*exp((1/2)*(5+33^(1/2))*t)*33^(1/2)+(3/4)*_C2*exp(-(1/2)*(-5+33^(1/2))*t)+(3/4)*_C1*exp((1/2)*(5+33^(1/2))*t)}


 

Download mw.mw

Another way...

Rouben Rostamian 8636
August 22 2019
1 0

As tomleslie said, there are many ways of achieving the result.  Here is a straightforwad "bare hands" way which may be more accessible to a beginner.

> 

restart;

Angles of the two vectors with respect to the vertical, in radians:

> 

a1 := convert(30.0*degrees, radians);
a2 := convert(60.0*degrees, radians);

.5235987757

1.047197551

The two vectors:

> 

v1 := 1.0 * < sin(a1), cos(a1) >;
v2 := 1.5 * < sin(a2), cos(a2) >;

Vector(2, {(1) = .5000000001, (2) = .8660254037})

Vector[column](%id = 18446884770620611990)

Their sum:

> 

v := v1 + v2;

Vector(2, {(1) = 1.79903810565, (2) = 1.6160254040000002})

v's length

> 

sqrt(v[1]^2 + v[2]^2);

HFloat(2.418279597555689)

v's angle with respect to the vertical in radians and in degrees

> 

arctan(v[1],v[2]);
convert(%, degrees);

HFloat(0.838936788996444)

HFloat(48.06753728147721)*degrees

 

 

 

Download mw.mw

Loop construction...

Rouben Rostamian 8636
August 18 2019
0 0

Probably this is what you want:

for n in (10*2^i $i=0..4) do
    print(n);  # or whatever
end do;

or equivalently

for i from 0 to 4 do
    n := 10*2^i;
    print(n); # or whatever
end do;

or

n := 10; 
while n <= 160 do
    print(n); # or whatever
    n := 2*n;
end do;

Pick whichever you like.

 

Answer...

Rouben Rostamian 8636
August 16 2019
2 3

 

Here is a possible way.

 

Example 1

 
> 

restart;
ode:=3*diff(y(x),x)^2+diff(y(x),x)^3+sin(x)+y(x)=x*y(x)+x*diff(y(x),x);

3*(diff(y(x), x))^2+(diff(y(x), x))^3+sin(x)+y(x) = x*y(x)+x*(diff(y(x), x))

> 

indets['flat'](ode,{`^`('identical'(diff(y(x),x)),'algebraic'),'identical'(diff(y(x),x))})

{(diff(y(x), x))^2, (diff(y(x), x))^3, diff(y(x), x)}

> 

ode_wanted:= 3*diff(y(x),x)^2+diff(y(x),x)^3-x*diff(y(x),x)=-sin(x)-y(x)+x*y(x)

3*(diff(y(x), x))^2+(diff(y(x), x))^3-x*(diff(y(x), x)) = -sin(x)-y(x)+x*y(x)

> 

selectremove(has, (lhs-rhs)~(ode), diff(y(x),x));
ode_new := %[1] = - %[2];
ode_new - ode_wanted;

3*(diff(y(x), x))^2+(diff(y(x), x))^3-x*(diff(y(x), x)), sin(x)+y(x)-x*y(x)

3*(diff(y(x), x))^2+(diff(y(x), x))^3-x*(diff(y(x), x)) = -sin(x)-y(x)+x*y(x)

0 = 0

> 

 

> 

 

Example 2

 
> 

restart;
ode:=3*diff(y(x),x)^2+diff(y(x),x)=x*diff(y(x),x)+5;

3*(diff(y(x), x))^2+diff(y(x), x) = x*(diff(y(x), x))+5

> 

indets['flat'](ode,{`^`('identical'(diff(y(x),x)),'algebraic'),'identical'(diff(y(x),x))})

{(diff(y(x), x))^2, diff(y(x), x)}

> 

ode_wanted:= 3*diff(y(x),x)^2+diff(y(x),x)-x*diff(y(x),x)=5

3*(diff(y(x), x))^2+diff(y(x), x)-x*(diff(y(x), x)) = 5

> 

selectremove(has, (lhs-rhs)~(ode), diff(y(x),x));
ode_new := %[1] = - %[2];
ode_new - ode_wanted;

3*(diff(y(x), x))^2+diff(y(x), x)-x*(diff(y(x), x)), -5

3*(diff(y(x), x))^2+diff(y(x), x)-x*(diff(y(x), x)) = 5

0 = 0

> 

 


 

Download lhs_rhs-new.mw

 

Solution...

Rouben Rostamian 8636
August 14 2019
3 2

I have made several changes to your worksheet and I hope that I have not introduced any errors.  Check to be sure.

The optimal diffusion coefficient turns out to be d = 0.14 (approximately).

> 

restart;
> 

t_number:=<0, 10, 20, 30, 40>:
m_number:=<11.50000000, 4.641182547, 1.273311993, 0.361845238, 0.288711649>:
> 

 

> 

Q:=proc(d)
    local pds, i, S, PDE, IBC,
          L := 2, Mx0 := 0.05, cx0 := Mx0/(1-Mx0),
          Mdb_i := m_number[1], ct0 := Mdb_i;

    if not type(d, numeric) then return 'procname(_passed)' end if;

    PDE := diff(C(x,t),t)=d*diff(C(x,t),x,x);
    IBC := {C(x,0)=ct0, C(0,t)=cx0, D[1](C)(L,t)=0};

    pds := pdsolve(PDE,IBC,numeric);
    S := 0;
    for i from 1 to 5 do
        # solution at the desired time
        pds:-value(t=t_number[i], output=listprocedure);
        eval(C(x,t), %);              # extract the C(x,t) at that time
        int(%, 0..2, numeric) / L;    # compute the average of C(x,t) at that time
        S += (% - m_number[i])^2;     # accumulate the residuals
     end do;  
     return sqrt(S);
end proc:
> 

# Try out the proc
> 

Q(0.1), Q(0.2), Q(0.5), Q(1.0);

HFloat(2.0125989798253805), HFloat(1.9268252656215228), HFloat(4.344399869795477), HFloat(4.74473396546933)

> 

plot(Q(d), d=0.1 .. 0.3, numpoints=5, view=0..2);

> 

# We conclude that the best choice of d is approximately 0.14
> 

 


 

Download zz.mw

 

Answer...

Rouben Rostamian 8636
August 11 2019
1 1

To plot the derivatives, add the output=operator option to your dsolve(), as in
p := dsolve(%, vars, numeric, output = operator);

Do odeplot() as before to plot the solutions.  To plot the derivative of y(t) (also those of x[1](t) and x[2](t)) do
eval(diff(y(t), t), Sys);
eval(%, p);
plot(%, t = 0 .. 1500);

As to the higher order derivatives, y'', y''', etc., these are not quite well-defined because the right-hand sides of your differential equations contain discontinuities, which means that y' is discontinuous.  You don't want to differentiate a discontinuous function, do you?

You can actually see the discontinuity in y' by changing the plot(%, t = 0 .. 1500) in what I have shown to plot(%, t = 500 .. 1500).

As to the plotting of the bar graph, I don't think it's difficult but I don't have the time to go into it right now.  Perhaps someone else will.

 

Applying odetest...

Rouben Rostamian 8636
August 07 2019
1 0

Of the two solutions returned by dsolve() one corresponds to positive y and the other to negative y.

restart;
de := diff(y(x),x) = abs(y(x)) + 1;
                         d                   
                  de := --- y(x) = |y(x)| + 1
                         dx                  
dsol := dsolve(de);
                       exp(-x)                           
      dsol := y(x) = - ------- + 1, y(x) = exp(x) _C1 - 1
                         _C1                             
odetest(dsol[1], de) assuming y(x)<0;
                               0
odetest(dsol[2], de) assuming y(x)>0;
                               0

Solution...

Rouben Rostamian 8636
August 07 2019
5 5

Your differential equation has y(x) = 0 as a solution, and that's what
Maple returns by default.  However, for special choices of omega^2 there are
nonzero solutions.  Those special choices of omega^2 are called the problem's
eigenvalues, and their corresponding solutions are called the eigenfunctions.

Here is how we go about finding the eigenvalues and eigenfunctions.

We begin with a couple of self-evident observations.
   

1. 

If a function y(x) is an eigenfunction, then for any nonzero constant c
the function c*y(x) is also an eigenfunction.
2. 

The derivative "y '(0)" of an eigenfunction cannot be zero because the
conditions y(0) = 0 and "y '(0)=0" will imply that `&equiv;`(y(x), 0) which
is not possible since an eigenfunction is a nonzero function by definition.


Putting those two observations together, we may take "y '(0)=1" for
all eigenfunctions without a loss of generality.

This adds an extra boundary condition over and above what you
already have.  That gives us the flexibility of introducing a new
unknown in the system.  We take the new unknown to be the
constantomega and we solve the system for the unknowns {omega, y(x)}.

> 

restart;

The differential equation

> 

de := T*diff(y(x), x, x) + rho*omega^2*y(x) = 0;

T*(diff(diff(y(x), x), x))+rho*omega^2*y(x) = 0

The boundary conditions (note the extra third condition!)

> 

bc := y(0) = 0, y(L) = 0, D(y)(0)=1;

y(0) = 0, y(L) = 0, (D(y))(0) = 1

> 

dsol_tmp := dsolve({de,bc}, {y(x),omega});

{omega = Pi*(2*_Z1+_B1)*T^(1/2)/(rho^(1/2)*L), y(x) = L*sin(Pi*(2*_Z1+_B1)*x/L)/(Pi*(2*_Z1+_B1))}

Here omega and y(x) are expressed in terms of the arbitrary
constants _Z1 and _B1.  Let's see what they are:

> 

about(_Z1);
about(_B1);

Originally _Z1, renamed _Z1~:

  is assumed to be: integer

Originally _B1, renamed _B1~:
  is assumed to be: OrProp(0,1)

OK then.  _Z1 is an integer and _B1 is either zero or one.  It follows
that the combination 2*_Z1+_B1 which enters the solution is an
arbitrary integer.  We call it n

> 

dsol := eval(dsol_tmp, {2*_Z1+_B1=n});

{omega = Pi*n*T^(1/2)/(rho^(1/2)*L), y(x) = L*sin(Pi*n*x/L)/(Pi*n)}

Let's verify that this satisfies the differential equation

> 

eval(de, dsol);

0 = 0


 

Download mw.mw

Solution...

Rouben Rostamian 8636
August 06 2019
1 4

I don't understand what you mean by "returns only one option".  Here is the full solution.

> 

restart;
> 

wW := unapply(piecewise(0 <= x and x < L/2, 0, L/2 <= x and x <= L, w[0]), x);

wW := proc (x) options operator, arrow; piecewise(0 <= x and x < (1/2)*L, 0, (1/2)*L <= x and x <= L, w[0]) end proc

> 

eq := k*diff(y(x), x$4) = wW(x);

eq := k*(diff(y(x), x, x, x, x)) = piecewise(0 <= x and x < (1/2)*L, 0, (1/2)*L <= x and x <= L, w[0])

> 

dsolve({eq, y(0) = 0, D(y)(0) = 0, (D@@2)(y)(L) = 0, (D@@3)(y)(L)=0}, y(x))
  assuming 0 < L;

y(x) = piecewise(x < (1/2)*L, -L*w[0]*x^3/(12*k)+3*L^2*w[0]*x^2/(16*k), x < L, -L*w[0]*x^3/(6*k)+L^2*w[0]*x^2/(4*k)-L^3*w[0]*x/(48*k)+L^4*w[0]/(384*k)+w[0]*x^4/(24*k), L <= x, 7*L^3*w[0]*x/(48*k)-5*L^4*w[0]/(128*k))


 

Download mw.mw

Need plain text...

Rouben Rostamian 8636
July 30 2019
1 2

In 2.mw we have
read "1.mw"; 
but Maple's read command is meant for reading plain text files, not Maple worksheets.

To fix, copy and paste the contents of 1.mw as a plain text into a file, let's say 1.maple.  Then change the read command to 
read "1.maple"; 

Here is 1.maple for your convenience.

 

Solution with numbers...

Rouben Rostamian 8636
July 30 2019
1 0

 

Approximate solution

 

On dividing through by r, your differential equation changes to
d/(r*dr)*(r*d*H(r)/dr)+(k^2-b^2*r/R^2)*H(r) = 0.``
The first term is the second derivative of H(r) in polar coordinates.

With your numerical parameters, which are
R = 370, k = 100, b = 35,
the coefficient of the second term is 10000-(49/5476)*r, which varies
between 9996 and 10000 on the range of interest 1/R < r and r < R, and
therefore is essentially a constant. Consequently, for practical purposes
we may take b = 0 and reduce the differential equation to
d/(r*dr)*(r*d*H(r)/dr)+k^2*H(r) = 0,
We expect the solution to be something like cos(k*r), that is, cos(100*r),
and therefore highly oscillatory.  We conclude that searching for an
approximation in the form of power series is hopeless.

The good news is that Maple can obtain a symbolic solution for the
reduced differential equation.  Here it is.

> 

restart;

> 

de := 1/r*Diff(r*Diff(H(r),r),r) + k^2*H(r);

(Diff(r*(Diff(H(r), r)), r))/r+k^2*H(r)

> 

bc := H(R)=0, D(H)(1/R)=R;

H(R) = 0, (D(H))(1/R) = R

> 

dsol := dsolve({de,bc}, H(r));

H(r) = -R*BesselY(0, k*R)*BesselJ(0, k*r)/(k*(BesselJ(1, k/R)*BesselY(0, k*R)-BesselY(1, k/R)*BesselJ(0, k*R)))+BesselJ(0, k*R)*R*BesselY(0, k*r)/(k*(BesselJ(1, k/R)*BesselY(0, k*R)-BesselY(1, k/R)*BesselJ(0, k*R)))

Apply the parameters and plot:

> 

R := 370:  k := 100:

> 

plot(rhs(dsol), r=1/R..R, view=-0.1 .. 0.1, numpoints=10000);

We cannot see the graph's details because it oscillates so much, but we do see
the details if we plot it on a narrow range:

> 

plot(rhs(dsol), r=100 .. 101);

which shows fast oscillations as expected.


 

Download de2.mw

Solution in sine series...

Rouben Rostamian 8636
July 29 2019
3 1
> 

restart;
> 

de_orig := r*diff(H(r),r,r) + diff(H(r),r) + (k^2*r - b^2/R^2*r^2)*H(r);

r*(diff(diff(H(r), r), r))+diff(H(r), r)+(k^2*r-b^2*r^2/R^2)*H(r)

> 

bc := H(R)=0, H(1/R)=R;

H(R) = 0, H(1/R) = R

We change the independent variable from r to rho so that domain R, 1/R is mapped to "(0,Pi)."
Additionally, we change the dependent variable from H to G so that G(0) = G((π) = 0.

> 

change_of_vars := { r = (R^2 - 1)*rho/(Pi*R) + 1/R,
                    H(r) = G(rho) + R/Pi*(Pi-rho) };

{r = (R^2-1)*rho/(Pi*R)+1/R, H(r) = G(rho)+R*(Pi-rho)/Pi}

> 

tmp1 := PDEtools:-dchange(change_of_vars, [de_orig,bc], [G(rho), rho]);

[((R^2-1)*rho/(Pi*R)+1/R)*(diff(diff(G(rho), rho), rho))*Pi^2*R^2/(R^2-1)^2+(diff(G(rho), rho)-R/Pi)*Pi*R/(R^2-1)+(k^2*((R^2-1)*rho/(Pi*R)+1/R)-b^2*((R^2-1)*rho/(Pi*R)+1/R)^2/R^2)*(G(rho)+R*(Pi-rho)/Pi), G(Pi) = 0, G(0)+R = R]

The new differential equation is:

> 

de := tmp1[1];

((R^2-1)*rho/(Pi*R)+1/R)*(diff(diff(G(rho), rho), rho))*Pi^2*R^2/(R^2-1)^2+(diff(G(rho), rho)-R/Pi)*Pi*R/(R^2-1)+(k^2*((R^2-1)*rho/(Pi*R)+1/R)-b^2*((R^2-1)*rho/(Pi*R)+1/R)^2/R^2)*(G(rho)+R*(Pi-rho)/Pi)

Let's verify that the boundary conditions are G(0) = G(Pi) and G(Pi) = 0, as expected:

> 

tmp1[2];

G(Pi) = 0

> 

tmp1[3];

G(0)+R = R

Now that the domain is 0, Pi and the boundary conditions are zero, we may
express the solution in an infinite sine series:

> 

S := sum(a[n]*sin(n*rho), n=1..infinity);

sum(a[n]*sin(n*rho), n = 1 .. infinity)

To determine the coefficients a__n, plug the series into the differential equation:

> 

eval(de, G(rho)=S):
tmp2 := combine(%);

(sum((-sin(n*rho)*Pi^4*R^7*n^2*rho*a[n]+sin(n*rho)*Pi^2*R^9*k^2*rho*a[n]-sin(n*rho)*Pi*R^8*b^2*rho^2*a[n]+cos(n*rho)*Pi^4*R^7*n*a[n]-sin(n*rho)*Pi^5*R^5*n^2*a[n]+sin(n*rho)*Pi^4*R^5*n^2*rho*a[n]+sin(n*rho)*Pi^3*R^7*k^2*a[n]-3*sin(n*rho)*Pi^2*R^7*k^2*rho*a[n]-2*sin(n*rho)*Pi^2*R^6*b^2*rho*a[n]+4*sin(n*rho)*Pi*R^6*b^2*rho^2*a[n]-Pi^4*R^5*a[n]*n*cos(n*rho)-2*sin(n*rho)*Pi^3*R^5*k^2*a[n]+3*Pi^2*R^5*k^2*rho*a[n]*sin(n*rho)-sin(n*rho)*Pi^3*R^4*b^2*a[n]+6*sin(n*rho)*Pi^2*R^4*b^2*rho*a[n]-6*rho^2*R^4*b^2*Pi*a[n]*sin(n*rho)+Pi^3*R^3*k^2*a[n]*sin(n*rho)-R^3*k^2*rho*Pi^2*a[n]*sin(n*rho)+2*sin(n*rho)*Pi^3*R^2*b^2*a[n]-6*R^2*b^2*rho*Pi^2*a[n]*sin(n*rho)+4*rho^2*R^2*b^2*Pi*a[n]*sin(n*rho)-b^2*Pi^3*a[n]*sin(n*rho)+2*b^2*Pi^2*a[n]*sin(n*rho)*rho-b^2*Pi*a[n]*sin(n*rho)*rho^2)/(R^2-1), n = 1 .. infinity)-Pi*R^8*k^2*rho^2-Pi*R^7*b^2*rho^2-3*Pi^2*R^6*k^2*rho+2*Pi*R^6*k^2*rho^2-2*Pi^2*R^5*b^2*rho+5*Pi*R^5*b^2*rho^2+2*Pi^2*R^4*k^2*rho-Pi*R^4*k^2*rho^2+5*Pi^2*R^3*b^2*rho-7*Pi*R^3*b^2*rho^2-3*Pi^2*R*b^2*rho+3*Pi*R*b^2*rho^2+Pi^2*R^8*k^2*rho-3*R^5*b^2*rho^3+3*R^3*b^2*rho^3-R*b^2*rho^3-Pi^3*R^4*k^2-Pi^3*R^3*b^2+R^7*b^2*rho^3+Pi^3*R^6*k^2-Pi^3*R^6+Pi^3*R*b^2)/(Pi^3*R^6-Pi^3*R^4)

We want to determine the coefficients a__n so that tmp2 is identically zero on 0 < rho and rho < Pi.
For that purpose, we multiple tmp2 by sin(m*rho) and integrate the result from 0 to Pi.
That results in equations of the form
`&alpha;__m`+`&beta;__m`*a__m+sum(`&gamma;__m,n`*a__n, n = 1 .. infinity) = 0,    m = 1, 2, () .. ()
where the constants `&alpha;__m`, `&beta;__m`, `&gamma;__m,n` can be computed explicitly.  This is a coupled linear system
of infinitely many equations in the infinitely many unknowns a__m.  We cannot afford solving
that system but we can do the next best thing, which is to replace infinity with a finite number N 
and solve the resulting linear system of N equations in N unknowns.  That way you may obtain
a solution to any desired degree of accuracy.

The calculation is not difficult but it is tedious, so I leave it to you to finish it.

 


Download de.mw

Solution through finite differences...

Rouben Rostamian 8636
July 28 2019
3 1

Maple 2019 does not have a numeric solver for elliptic PDEs.
That does not prevent us from writing our own finite-difference
code.  Here is an illustration of solving Laplace's equation in
polar coordinates.

> 

restart;
> 

with(LinearAlgebra):
> 

with(plots):

The domain is a < r and r < b, c < t and t < d in polar coordinates.
For illustration we take:

> 

a := 1:  b :=2:
c := Pi/6:  d := Pi/2:

We subdivide the r interval into n__r subintervals of equal lengths.
We subdivide the tinterval into `n__t `subintervals of equal lengths.

> 

n__r := 10;
n__t := 15;

10

15

Calculate the resulting mesh sizes in the r and t directions:

> 

delta__r := (b-a)/n__r;
delta__t := (d-c)/n__t;

1/10

(1/45)*Pi

Here is what the mesh looks like:

> 

seq(seq([(a+i*delta__r)*cos(c+j*delta__t), (a+i*delta__r)*sin(c+j*delta__t)],
    i=0..n__r), j=0..n__t):
pointplot([%], symbol=solidcircle, color="Green", scaling=constrained);

We wish to solve the PDE

> 

pde := diff(u(r,t),r,r) + diff(u(r,t),r)/r + diff(u(r,t),t,t)/r^2 =  0;

diff(diff(u(r, t), r), r)+(diff(u(r, t), r))/r+(diff(diff(u(r, t), t), t))/r^2 = 0

on that domain, subject to the Dirichlet boundary conditions
"u(r,c)=`alpha__1`(r), u(r,d)=`alpha__2`(r), a<r<b,"
u(a, t) = `&alpha;__3`(t), u(b, t) = `&alpha;__4`(t), c < t and t < d.

We discretize the PDE in the interior points through central differencing.
Here U[i, j] represents the solution at the point of index i, j

> 

discretized_pde :=
   (U[i-1,j] - 2*U[i,j] + U[i+1,j])/delta__r^2
 + 1/(a+i*delta__r)*(U[i+1,j] - U[i-1,j])/(2*delta__r)
 + 1/(a+i*delta__r)^2*(U[i,j-1] - 2*U[i,j] + U[i,j+1])/delta__t^2 = 0;

100*U[i-1, j]-200*U[i, j]+100*U[i+1, j]+5*(U[i+1, j]-U[i-1, j])/((1/10)*i+1)+2025*(U[i, j-1]-2*U[i, j]+U[i, j+1])/(((1/10)*i+1)^2*Pi^2) = 0

That discretization results in (n__r-1)*(n__t-1) equations in the unknowns U[i, j]:

> 

eq_interior := seq(seq(discretized_pde, i=1..n__r-1), j=1..n__t-1):

Additionally, we have 2*n__t+2*n__r equations that assign boundary values to U[i, j] 

> 

eq_boundary :=
  seq(U[i,0]=1,                    i=1..n__r-1),        # u(r,c)=1
  seq(U[i,n__t]=0,                 i=1..n__r-1),        # u(r,d)=0
  seq(U[0,j]=0,                    j=0..n__t),          # u(a,t)=0
  seq(U[n__r,j]=c+j*delta__t,      j=0..n__t):          # u(b,t)=t

Altogether, these provide a system of (n__r+1)*(n__t+1)linear equations
in the (n__r+1)*(n__t+1) unknowns U[i, j]

> 

sys := eq_interior, eq_boundary:
> 

vars := seq(seq(U[i,j], j=0..n__t), i=0..n__r):
> 

nops([sys]), nops([vars]);

176, 176

Extract the system's coefficient matrix:

> 

LinSys := GenerateMatrix(evalf([sys]), [vars], augmented,
    datatype=float[8], storage=sparse):

Solve the system:

> 

V := LinearSolve(LinSys):

The vector Vcalculated above is of length (n__r+1)*(n__t+1).  We reformat
its entries into an `&x`(n__r+1, n__t+1)table T:

> 

Matrix(n__r+1, n__t+1, convert(V, list), datatype=float[8]):
T := rtable_redim(%, 0..n__r, 0..n__t):

Finally, we plot the solution:

> 

dataplot(a..b, c..d, T, coords=cylindrical,
    scaling=constrained, orientation=[-72,50,0]);


 

Download finite-difference-elliptic-polar.mw

typesetprime...

Rouben Rostamian 8636
July 27 2019
4 0
> 

restart;
> 

kernelopts(version);

`Maple 2019.1, X86 64 LINUX, May 21 2019, Build ID 1399874`

> 

Typesetting:-Settings(typesetprime=true):
> 

diff(y(x),x);

diff(y(x), x)

> 

diff(y(x),x,x);

diff(diff(y(x), x), x)

 

First 26 27 28 29 30 31 32 Last Page 28 of 58