Advanced Search

Rouben Rostamian

8264 Reputation

23 Badges

18 years, 295 days
Contact Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

Try this simpler one...

Rouben Rostamian 8264
July 15 2023
6 4

If you don't know how so solve a problem, try a simpler one first.  See if you understand the following before you tackle your messy case.

> 

restart;
> 

with(LinearAlgebra):
> 

Digits := 30;

30

 (1)

Make a random matrix; it is very likely to have a nonzero determinant

and therefore be nonsingular:

> 

M__exact := RandomMatrix(4);
Determinant(M__exact);

Matrix(4, 4, {(1, 1) = 50, (1, 2) = -50, (1, 3) = -38, (1, 4) = -98, (2, 1) = 10, (2, 2) = -22, (2, 3) = -18, (2, 4) = -77, (3, 1) = -16, (3, 2) = 45, (3, 3) = 87, (3, 4) = 57, (4, 1) = -9, (4, 2) = -81, (4, 3) = 33, (4, 4) = 27})

  

-20185224

 (2)

Replace the (1,2)  element by x, and determine x to make M__exact into a singular matrix:

> 

M__exact[1,2] := x:
Determinant(M__exact):
solve(%, x):
M__exact[1,2] := %;

6503458/4499

 (3)

Now M__exact is singular:

> 

Determinant(M__exact);

0

 (4)

Reduce the accuracy by truncating the entries from 30 to 20 decimal points:

> 

M := evalf[20](M__exact);

Matrix(4, 4, {(1, 1) = 50., (1, 2) = 1445.5341186930428984, (1, 3) = -38., (1, 4) = -98., (2, 1) = 10., (2, 2) = -22., (2, 3) = -18., (2, 4) = -77., (3, 1) = -16., (3, 2) = 45., (3, 3) = 87., (3, 4) = 57., (4, 1) = -9., (4, 2) = -81., (4, 3) = 33., (4, 4) = 27.})

 (5)

Now M is approximately singular,

> 

Determinant(M);

-0.2952e-12

 (6)

but from the point of view of Digits = 30, M is not singular, and therefore

the null space is empty:

> 

NullSpace(M);

{}

 (7)

and that's the issue that you are having in your worksheet.
To calculate the null space of the approximately singular matrix M,

apply the singular value decomposition:

> 

S, Vt := SingularValues(M, output=['S','Vt']);

S, Vt := Vector(4, {(1) = 1453.21753260898162930059463838, (2) = 132.138542019696590676057858086, (3) = 47.0223025429904652916635508744, (4) = 3.26928153214428999570938012074*(1/100000000000000000000)}), Matrix(4, 4, {(1, 1) = -0.343263390130779463921470892225e-1, (1, 2) = -.996848014529597489318515478287, (1, 3) = 0.256899695701266591875628064244e-1, (1, 4) = 0.667515081615734136773577608453e-1, (2, 1) = -.160139286909303325787418878728, (2, 2) = 0.710671575010807034257774212978e-1, (2, 3) = .670551592300163005864255159870, (2, 4) = .720878235194136232433643942038, (3, 1) = 0.196400352119669472946663977236e-1, (3, 2) = 0.261339977153148193494660351144e-1, (3, 3) = -.730920435970437612722419386768, (3, 4) = .681679249692312424157563481466, (4, 1) = .986301870755100907032581022014, (4, 2) = -0.236750758506050268946221115509e-1, (4, 3) = .124321775276849024311433861988, (4, 4) = .105793226250403080836382592166})

 (8)

The returned vector S is the list of the singular values of M.  We wee that

the last singular value is of the order of 0.1e-19, which is almost

zero.  That says that the last row of the matrix Vt is the basis

of the corresponding approximate null space.  Let's check:

> 

null__vec := Vt[-1]^+;

Vector(4, {(1) = .986301870755100907032581022014, (2) = -0.236750758506050268946221115509e-1, (3) = .124321775276849024311433861988, (4) = .105793226250403080836382592166})

 (9)
> 

M . null__vec;

Vector(4, {(1) = 0.2064116200e-20, (2) = 0.2275617400e-21, (3) = -0.1073762554e-19, (4) = 0.3080927157e-19})

 (10)

We see that the product is essentially the zero vector, which confirms

that null__vec is in the null space.

 
 

Download singular_values.mw

 

Corrections...

Rouben Rostamian 8264
July 12 2023
5 6
 
> 

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

 

The proper demo...

Rouben Rostamian 8264
June 30 2023
3 2

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

 

There is no neat solution...

Rouben Rostamian 8264
June 23 2023
2 2
 
> 

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

A clean way...

Rouben Rostamian 8264
June 21 2023
1 1

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

 

 

 

Use "select"...

Rouben Rostamian 8264
June 20 2023
1 1 
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.

Start with this and then build upon it....

Rouben Rostamian 8264
June 13 2023
1 7

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

 

Try this...

Rouben Rostamian 8264
June 13 2023
0 0
 
> 

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

 

 

Use the AllSolutions option to int...

Rouben Rostamian 8264
June 12 2023
2 1
> 

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)

Without the goto...

Rouben Rostamian 8264
June 04 2023
2 2

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

 

Maybe this?...

Rouben Rostamian 8264
June 03 2023
3 3

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

 

Here is how...

Rouben Rostamian 8264
May 21 2023
1 7
> 

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

Do it in cylindrical...

Rouben Rostamian 8264
May 20 2023
2 0

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

Perhaps this?...

Rouben Rostamian 8264
May 19 2023
0 0

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

 

Maybe this?...

Rouben Rostamian 8264
May 19 2023
0 0

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

 

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