Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

 

restart;

You are solving this initial value problem numerically:

f := r -> 1 - 1/r:
eqq := {RR(2) = 2 + log(1), diff(RR(r), r) = 1/f(r)};

{RR(2) = 2, diff(RR(r), r) = 1/(1-1/r)}

and then you write a (very inefficient) procedure for finding the inverse function
of that solution numerically.

Why not solve for the inverse right from the beginning?  It can be done
symbolically and it's very fast!  To explain, let me change your RR notation

to a single-letter symbol q. Then the differential equation is
dq/dr = 1/(1-1/r)that is dr/dq = 1-1/r.   Your initial condition is r(2) = 2"." 

de := diff(r(q),q) = 1 - 1/r(q);  ic := r(2)=2;

diff(r(q), q) = 1-1/r(q)

r(2) = 2

dsol := dsolve({de, ic});

r(q) = LambertW(exp(q-1))+1

There you have it.  That's equivalent to your numerical procedure "rt."

I will call it RT.

RT := unapply(rhs(dsol), q);

proc (q) options operator, arrow; LambertW(exp(q-1))+1 end proc

Plot your rt(R)when you get the chance, and compare it with:

plot(RT(R), R=-4..12);

The rest of your worksheet works correctly as is.

 

Download mw.mw

 

Your attempt to animate the tan function misses an essential ingredient.

Consider the trigonometric circle (a unit circle centered at the origin in Cartesian coordinates) and a ray connecting the origin O = (0,0) to a point P = (cos(t), sin(t)).  Then, the horizontal projection of P is (cos(t),0) and the vertical projection is (0,sin(t)).  This enables us to associate the sine and cosine functions with the diagram's geometry.

Question: Where does the tan function fit in that diagram?

Answer:  Draw a tangent line L to the circle at (1,0).  Extend the ray OP to intersect the line L at some point H. You will see right away that H = (1, tan(t)).  So the tan function is associated with the point H in our diagram.  This is an absolutely essential information that should be conveyed to the students when teaching them the trigonometric circle.  [The function cot also has a similar geometric interpretation; I will leave it to you to figure it out.]

The animation below demonstrates the tan function properly in the context of the trigonometric circle.  Note the vertical tangent line to the unit circle.

AnimationCercleTrigoFonctiTrigo-rr.mw

 

 

restart;
q1:=1-1/(w**2-3*sigma*k**2)-deltab*mu/((w-k*u0b)**2-3*mu*deltab*sigmab*k**2)+A/k**2=0;

1-1/(-3*k^2*sigma+w^2)-deltab*mu/((-k*u0b+w)^2-3*mu*deltab*sigmab*k^2)+A/k^2 = 0

You want to solve for w/k, so let Q = w/k.  Then w = k*Q, and the equation changes to

q2 := subs(w=Q*k, q1);

1-1/(Q^2*k^2-3*k^2*sigma)-deltab*mu/((Q*k-k*u0b)^2-3*mu*deltab*sigmab*k^2)+A/k^2 = 0

You want to solve eq2 for Q.  Let's try a very special case where most of the parameters are taken to be 1:

q3 := subs(deltab=1, sigmab=1, u0b=1, sigma=1, A=1, q2);

1-1/(Q^2*k^2-3*k^2)-mu/((Q*k-k)^2-3*mu*k^2)+1/k^2 = 0

solve(q3, Q);

RootOf((k^2+1)*_Z^4+(-2*k^2-2)*_Z^3+(-3*k^2*mu-2*k^2-4*mu-3)*_Z^2+(6*k^2+8)*_Z+9*mu*k^2-3*k^2+15*mu-4)

We see that even in this very special case, the solution Q is a root of a fourth degree polynomial.

We conclude that there is no hope of obtaining a useful symbolic solution for the general case.

 

Download mw.mw

Acer has offered an answer by modifying your worksheet as you have requested.  But your worksheet is much more complicated than it needs be.  Here is a clean way of producing your diagram.  You should be able to merge this with the improvements suggested by acer.

restart;

with(plots):

with(plottools):

angles := { seq(i*Pi/4, i=0..7), seq(i*Pi/6, i=0..11) };

{0, Pi, (1/2)*Pi, (1/3)*Pi, (1/4)*Pi, (1/6)*Pi, (2/3)*Pi, (3/2)*Pi, (3/4)*Pi, (4/3)*Pi, (5/3)*Pi, (5/4)*Pi, (5/6)*Pi, (7/4)*Pi, (7/6)*Pi, (11/6)*Pi}

Distance of angle lables from the center.  Adjust to taste.

r := 1.12;

1.12

display(
        circle(color=blue, thickness=2),
        seq(textplot([r,a,a],coords=polar), a in angles),
        seq(line([0,0], [cos(a),sin(a)], color=red, linestyle=dash), a in angles),
        seq(line([0,sin(a)],[cos(a),sin(a)], color=gray), a in {Pi/6, Pi/4, Pi/3}),
        seq(line([cos(a),0],[cos(a),sin(a)], color=gray), a in {Pi/6, Pi/4, Pi/3}),
        textplot([1.2,0,typeset(", ", 2*Pi)], align=right),
        textplot([[1/2,0,1/2],[sqrt(2)/2,0,sqrt(2)/`2`],[sqrt(3)/2,0,sqrt(3)/`2`]], align=below),
        textplot([[0,1/2,1/2],[0,sqrt(2)/2,sqrt(2)/2],[0,sqrt(3)/2,sqrt(3)/2]], align=left),
axes=none, scaling=constrained, size=[800,800]);


 

Download mw.mw

 

 

 

restart;
A := {3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101};
select(x -> is(modp(x,3)=1), A);

The answer is

{7, 13, 19, 31, 37, 43, 61, 67, 73, 79, 97}

By the way, if you intend A to represent the list of primes up to 101, then A should include 2 as well.  Two is the only even prime.

As I noted earlier, it is a good idea to start with a simple case and build upon it after you gain confidence.  This toy example should be enough to get you started.

restart;

n := 2;

2

Z := Array(0..n-1,0..n-1, (i,j) -> z[i,j](t));

`Array(0..1, 0..1, {(1, 1) = z[0, 0](t)})`

I have defined the coefficients arbitrarily.  Change as needed.

a := (i,j,k,l) -> 1*i + 2*j + 3*k + 4*l;

proc (i, j, k, l) options operator, arrow; i+2*j+3*k+4*l end proc

The (i,j) entry of the matrix on the right-hand side of the system of differential equations

add(add(a(i,j,k,l)*z[k,l](t), k=0..n-1), l=0..n-1):
ij_entry := unapply(%, i, j);

proc (i, j) options operator, arrow; (2*j+i)*z[0, 0](t)+(i+2*j+3)*z[1, 0](t)+(i+2*j+4)*z[0, 1](t)+(i+2*j+7)*z[1, 1](t) end proc

Here is the right-hand side of the system of differential equations:

the_rhs := Array(0..n-1,0..n-1, ij_entry);

`Array(0..1, 0..1, {(1, 1) = \`+\`(\`*\`(3, \`*\`(z[1, 0](t))), \`*\`(4, \`*\`(z[0, 1](t))), \`*\`(7, \`*\`(z[1, 1](t))))})`

The system of differential equations:

DE := I*diff(Z, t) = the_rhs;

`Array(0..1, 0..1, {(1, 1) = \`*\`(I, \`*\`(diff(z[0, 0](t), t)))})` = `Array(0..1, 0..1, {(1, 1) = \`+\`(\`*\`(3, \`*\`(z[1, 0](t))), \`*\`(4, \`*\`(z[0, 1](t))), \`*\`(7, \`*\`(z[1, 1](t))))})`

Initial conditions

IC := eval(Z, t=0) =
        Array(0..n-1, 0..n-1, (i,j) -> Physics:-KroneckerDelta[i,j]);

`Array(0..1, 0..1, {(1, 1) = z[0, 0](0)})` = `Array(0..1, 0..1, {(1, 1) = 1})`

Solve the system for the unknowns z__ij(t):

dsol := dsolve({DE, IC});

{z[0, 0](t) = (1/3)*exp(-(22*I)*t)+(1/6)*exp((2*I)*t)+1/2, z[0, 1](t) = (1/2)*exp(-(22*I)*t)-1/2, z[1, 0](t) = (5/12)*exp(-(22*I)*t)+(1/12)*exp((2*I)*t)-1/2, z[1, 1](t) = (7/12)*exp(-(22*I)*t)-(1/12)*exp((2*I)*t)+1/2}

Here is the solution Z:

answer := simplify(subs(dsol, Z));

`Array(0..1, 0..1, {(1, 1) = \`+\`(\`*\`(\`/\`(1, 3), \`*\`(exp(\`+\`(\`-\`(\`*\`(\`+\`(\`*\`(22, \`*\`(I))), \`*\`(t))))))), \`*\`(\`/\`(1, 6), \`*\`(exp(\`*\`(\`*\`(2, \`*\`(I)), \`*\`(t))))), \`/\`(1, 2))})`

Download mw.mw

 

 

restart;

sigma :=  t-> (47.965 + 3.08*0.2*t^(1 - 0.7)/GAMMA(3 - 0.7))*epsilon;

proc (t) options operator, arrow; (47.965+.616*t^.3/GAMMA(2.3))*epsilon end proc

Not much happens between t = 1 and t = 100:

plot([sigma(1), sigma(100)], epsilon=0..1);

To see a noticeable change in time, look for the time interval of zero to a million.
Here is what you get if you look at uniformly spaced time intervals:

tvals := seq(i, i=0..1e6, 1e5);

0, 0.1e6, 0.2e6, 0.3e6, 0.4e6, 0.5e6, 0.6e6, 0.7e6, 0.8e6, 0.9e6, 0.10e7

plot([seq(sigma(t), t in tvals)], epsilon=0..1);

And here is what you get with a variable time spacing:

tvals := seq(i^(1/(1-0.7)), i=0..100, 10);

0, 2154.434688, 21715.34091, 83895.27757, 218876.9209, 460503.9367, 845611.4093, 1413600.856, 2206141.119, 3266944.055, 4641588.826

plot([seq(sigma(t), t in tvals)], epsilon=0..1);

 

 

Download mw.mw

 

 

restart;

phi:= (mu,Q2)->sqrt(2/l)*sin(mu*Pi*(Q2+l/2)/l);  

proc (mu, Q2) options operator, arrow; sqrt(2/l)*sin(mu*Pi*(Q2+(1/2)*l)/l) end proc

A := diff(phi(mu,Q2),Q2);
B := diff(phi(nu,Q2),Q2);

2^(1/2)*(1/l)^(1/2)*mu*Pi*cos(mu*Pi*(Q2+(1/2)*l)/l)/l

2^(1/2)*(1/l)^(1/2)*nu*Pi*cos(nu*Pi*(Q2+(1/2)*l)/l)/l

-1/2/m__2*int(A*B, Q2=-l/2..l/2, AllSolutions)
    assuming mu::posint, nu::posint:

fh1 := combine(%);
 

fh1 := -piecewise(mu = nu, nu^2*Pi^2/(2*l^2), 0)/m__2

Download mw.mw

PS: Replacing combine with convert on the last line yields a better looking result:

-1/2/m__2*int(A*B, Q2=-l/2..l/2, AllSolutions)
                assuming mu::posint, nu::posint:
fh1 := convert(%, piecewise, mu);
 

fh1 := piecewise(mu = nu, -nu^2*Pi^2/(2*m__2*l^2), 0)

I don't know how the undocumented goto works.  Here is a  version without goto.  Have a look at the help page for break.

restart;

for i to 4 do
        for j to 4 do
                print([i,j]):
                if i = 2 and j = 3 then
                        break i;  # break out of the 'i' loop
                end if;
        end do;
end do;

[1, 1]

[1, 2]

[1, 3]

[1, 4]

[2, 1]

[2, 2]

[2, 3]

Download mw.mw

 

I see a lot of experiments in your worksheet but I didn't see a clear statement of the objective.  Perhaps this is what you have in mind?

restart;

de := diff(phi(s), s, s) + K*cos(phi(s)) = 0;

diff(diff(phi(s), s), s)+K*cos(phi(s)) = 0

bc := D(phi)(L) = phi__0;

(D(phi))(L) = phi__0

dsol := dsolve({de,bc}, phi(s), implicit);

Int(1/(-2*sin(_a)*K+2*sin(phi(L))*K+phi__0^2)^(1/2), _a = 0 .. phi(s))-s-c__2 = 0, Int(-1/(-2*sin(_a)*K+2*sin(phi(L))*K+phi__0^2)^(1/2), _a = 0 .. phi(s))-s-c__2 = 0

diff(dsol[1], s):  # either dsol[1] or dsol[2] will do
isolate(%, diff(phi(s), s))^2;

(diff(phi(s), s))^2 = -2*sin(phi(s))*K+2*sin(phi(L))*K+phi__0^2

 

Download mw.mw

 

restart;

f := t -> sin(t) - 1/2;

proc (t) options operator, arrow; sin(t)-1/2 end proc

plot(f(t), t=0..Pi);

plot([1, f(t) + 1], t=0..Pi, coords=polar, scaling=constrained, color=[blue,red], thickness=5);

Download mw.mw

What you have is not an error, but a bad choice of coordinates.  The intersection of a cone and plane is best described in cylindrical coordinates.  You will get a perfect ellipse that way.

You have the cone z=x^2+y^2, and the plane z = 3 + y/3.  Switch to cylindrical coordinates with x=r*cos(t), y=r*sin(t).   Then the equations of the cone and the plane become z = r^2 and z=3+1/3*r*sin(t).  At the intersection we have r^2 = 2 + 1/3*r*sin(t).  Solve this quadratic for r in terms of t.  Pick either of the two solutions. Call it r(t). Then the intersection curve is the ellipse described as a space curve:

[  r(t)*cos(t), r(t)*sin(t), 3 + 1/3*r(t)*sin(t)] ,    -Pi < t < Pi

Double-check your equation.  Perhaps "y" should be "u"?

If so, your boundary value problem does have a solution.  Here it is.

restart;

de := diff(u(x),x,x) + 1/x*diff(u(x),x) + u(x) = 6 - 9*x + x^2 - x^3;

diff(diff(u(x), x), x)+(diff(u(x), x))/x+u(x) = -x^3+x^2-9*x+6

Solve with the boundary condition u(1) = 0.

We will deal with u(0) = 0 afterward.

dsol := dsolve({de, u(1)=0});

u(x) = BesselJ(0, x)*c__2+BesselY(0, x)*(-BesselJ(0, 1)*c__2-2)/BesselY(0, 1)-x^3+x^2+2

We pick c__2 by trial-and-error (there ought to

be a better way!) to enforce u(0) = 0:

subs(c__2=-2.5948, dsol);
plot(rhs(%), x=0..1, color=red, thickness=5);

u(x) = -2.5948*BesselJ(0, x)+BesselY(0, x)*(2.5948*BesselJ(0, 1)-2)/BesselY(0, 1)-x^3+x^2+2

 

Download mm.mw

 

I can't see a meaning in what you have written.  Perhaps this is what you really want?

restart;

n := [ 82, 79, 84, 85, 76, 76, 78, 82, 85, 77, 83, 84, 87, 86, 87, 88, 76 ];

[82, 79, 84, 85, 76, 76, 78, 82, 85, 77, 83, 84, 87, 86, 87, 88, 76]

nops(n);

17

for i from 3 to nops(n) do
        f[i] := n[i-2] - 3*n[i-1] + 3*n[i];
end do;

97

82

57

85

82

88

87

58

103

80

92

81

90

89

51

Download mw.mw

 

In your worksheet you have things like d*u/dt.  That says "d" times "u" divided by "dt". That's not how one expresses a derivative in Maple. Considering that, I suggest that you begin with learning about the Maple syntax, solve a few simple ordinary differential equations, and then a few simple partial differential equations, and only then attempt to solve a nonlinear system of PDEs.

For whatever it's worth, I have converted your sketches of equations to proper Maple syntax and solved them. I can't tell how useful will this be to you.

restart;

pde1 := diff(u(y,t),t) - s*diff(u(y,t),y) = diff(u(y,t),y,y) + delta*v(y,t);

diff(u(y, t), t)-s*(diff(u(y, t), y)) = diff(diff(u(y, t), y), y)+delta*v(y, t)

pde2 := diff(v(y,t),t) - s*diff(v(y,t),y) = Pr*diff(v(y,t),y,y) + lambda*exp(v(y,t)/(epsilon*v(y,t)+1))/Pr;

diff(v(y, t), t)-s*(diff(v(y, t), y)) = Pr*(diff(diff(v(y, t), y), y))+lambda*exp(v(y, t)/(epsilon*v(y, t)+1))/Pr

ic := u(y,0)=0, v(y,0)=0;

u(y, 0) = 0, v(y, 0) = 0

bc := u(0,t)=0, v(0,t)=0, u(1,t)=1, v(1,t)=0;

u(0, t) = 0, v(0, t) = 0, u(1, t) = 1, v(1, t) = 0

s := 0.5;
epsilon := 0.01;
Pr := 7;
lambda := 1;
delta := 1;

.5

0.1e-1

7

1

1

 

pdsol := pdsolve({pde1,pde2}, {ic, bc}, numeric, spacestep=0.005);

_m139674373634432

pdsol:-plot3d(u, t=0..0.3);

pdsol:-plot3d(v, t=0..0.3);

 

Download mw.mw

 

4 5 6 7 8 9 10 Last Page 6 of 57