Rouben Rostamian

Heap size...

Answered: Rouben Rostamian 8259

October 19 2023

1 4

The default java heap size is 512. You may change it through the commandline flag −j, as in

maple -j 65536 -x filename.mw

I don't have an answer to your second question.

Maybe this?...

Answered: Rouben Rostamian 8259

October 05 2023

3 6

It's hard to understand what you are asking. Maybe this is what you are looking for?

Download mw.mw

PS: In this particular example the interpolating surface is vertical. You may construct an interpolating surface even when the two curves are not situated one directly above the other. Then the interpolating surface will be:

S := (1-z)*C1 + z*C2;
plot3d(S, z=0..1, x=-1..1, lightmodel=light4);

Here is how...

Answered: Rouben Rostamian 8259

September 17 2023

0 1

You didn't say maximum value of what, so here I show how to calculate the maxima (and minima) of S(t). Adjust as needed.

Continuation of the code in deplot3d-animated_SIA.mw.

Finding the critical points of S(t)

Here we will find the maxima (and minima) of S(t) when the solution starts at ICs8[1].

Change as needed.

The critical points of S occur when the derivative is zero. So introduce a fourth

differential equation, dS/dt = U, and solve the system of four differential equations:

>	de4 := diff(S(t),t) = U(t);

>	IC := ICs8[1][];

>	dsol := dsolve({de1,de2,de3,de4, IC}, numeric, output=operator);

>	my_S := subs(dsol, S); my_U := subs(dsol, U);

>	plot([my_S(t), my_U(t)], t=0..6, color=["Red", "Green"]);

The times when S(t) has a local maximum or minimum:

>	t_vals := [fsolve(my_U(t), t=0..6, maxsols=10)];

>	S_vals := map(my_S, t_vals);

You may determine which is min and which is max visually.

It's possible to automate that by intorducing a fifth differential

equation dU/dt = V, and picking min and max based on the

sign of V at the cricial points, but perhaps that's more that what

you would care to do.

>

deplot3d-animated_SIA-extra.mw

The problem is ill-posed...

Answered: Rouben Rostamian 8259

September 14 2023

3 0

Your initial value problem is ill-posed. Applying Dirac's delta at time

and specifying at the same time is not quite meaningful.

To see why, apply Dirac's delta not at , but at some unspecified .

Then solve your ODE:

>

restart;

>	de := D(y)(t) + y(t)/tau = Dirac(t-a)/tau;

>	ic := y(0) = 0;

>	sol := dsolve({de, ic});

y(t) = exp((-t+a)/tau)*(Heaviside(t-a)+Heaviside(a)-1)/tau

You are interested in the case of . But the limit of the solution

as goes to zero does not exist:

>	limit(sol, a=0, right);

y(t) = exp(-t/tau)*Heaviside(t)/tau

>	limit(sol, a=0, left);

y(t) = exp(-t/tau)*(Heaviside(t)-1)/tau

It seems that you favor one of the limits over the other but that

bias is only in your head; the mathematics that you are applying

is not aware of it.

Download mw.mw

Symbolic solution...

Answered: Rouben Rostamian 8259

September 03 2023

4 0

Maple can find a symbolic solution to your boundary value problem in terms of special functions. There is no need to resort to a numerical approximation.

PS: Note added later:

It appears that you have imposed the odd-looking boundary condition y(1e-14)=1 solely to avoid division by zero in your numerical calculation. The symbolic solution shows that the limit of the solution as x approaches zero from the right is 1, so you might as well change your boundary condition to y(0)=1. To verify that, set eps=0 in the worksheet that I posted earlier.

>

restart;

>	Digits := 30;

>	eps := 1.0e-14;

>	de := sqrt(x)diff(y(x),x,x) - 3/2y(x) - 1/2 = 0;

x^(1/2)*(diff(diff(y(x), x), x))-(3/2)*y(x)-1/2 = 0

>	dsol := dsolve({de, y(eps)=1, y(10)=0});

>	indets(dsol, name);

>	evalf(subs(x=eps, dsol));

>	evalf(subs(x=10, dsol));

>	plot(rhs(dsol), x=0..10);

>

Download zz.mw

A numerical solution...

Answered: Rouben Rostamian 8259

August 26 2023

1 7

If you do:

fsolve({Eq1=p_o-p_i,Eq2=F_a}, {A=0.167629..0.167634, B=1.000640 .. 1.000650});

You will get:

{A = 0.1676309260049029044111806984482988970179873516977323,
 B = 1.000644958082968954903634017487303572859626170659325}

I obtained the rough range of the roots by trial-and-error. There may be other roots.

Solution...

Answered: Rouben Rostamian 8259

August 24 2023

2 2

Calculating the tangent plane to a sphere is quite trivial. Why bother with the geometry package for it?

Just write a proc to receive a point on the sphere, calculate the tangent plane, and format the result according to your specification. Then map the proc over your list of points.

>

restart;

>	S := (x-1)^2 + (y-3)^2 + (z-5)^2 = 13^2;

>

pts := [[-12, 3, 5], [-11, -2, 5], [-11, -1, 2], [-11, -1, 8], [-11, 0,
        1], [-11, 0, 9], [-11, 3, 0], [-11, 3, 10], [-11, 6, 1], [-11, 6,
        9], [-11, 7, 2], [-11, 7, 8], [-11, 8, 5], [-4, -9, 5], [-4,
        3, -7], [-4, 3, 17], [-4, 15, 5], [-3, -9, 2], [-3, -9, 8], [-3,
        0, -7], [-3, 0, 17], [-3, 6, -7], [-3, 6, 17], [-3, 15, 2], [-3, 15,
        8], [-2, -9, 1], [-2, -9, 9], [-2, -1, -7], [-2, -1, 17], [-2,
        7, -7], [-2, 7, 17], [-2, 15, 1], [-2, 15, 9], [1, -10, 5], [1, -9,
        0], [1, -9, 10], [1, -2, -7], [1, -2, 17], [1, 3, -8], [1, 3,
        18], [1, 8, -7], [1, 8, 17], [1, 15, 0], [1, 15, 10], [1, 16,
        5], [4, -9, 1], [4, -9, 9], [4, -1, -7], [4, -1, 17], [4,
        7, -7], [4, 7, 17], [4, 15, 1], [4, 15, 9], [5, -9, 2], [5, -9,
        8], [5, 0, -7], [5, 0, 17], [5, 6, -7], [5, 6, 17], [5, 15, 2], [5,
        15, 8], [6, -9, 5], [6, 3, -7], [6, 3, 17], [6, 15, 5], [13, -2,
        5], [13, -1, 2], [13, -1, 8], [13, 0, 1], [13, 0, 9], [13, 3,
        0], [13, 3, 10], [13, 6, 1], [13, 6, 9], [13, 7, 2], [13, 7,
        8], [13, 8, 5], [14, 3, 5]]:

>

tangent_plane := proc(P)
        local C, N, ans;
        C := [1,3,5];      # the center of the sphere
        N := P - C;        # the normal vector at P
        (x - P[1])*N[1] + (y - P[2])*N[2] + (z - P[3])*N[3]; # <x,y,z> . N
        # divide the coefficients by their greatest common divisor
        igcd(coeffs(%));
        # we want the sign of the leading coefficient to be positive
        ans := %%/%;
        if coeff(%, x) > 0 then
                return [P,ans = 0];
        elif coeff(%, x) < 0 then
                return [P,-ans = 0];
        elif coeff(%, y) > 0 then
                return [P,ans = 0];
        elif coeff(%, y) < 0 then
                return [P,-ans = 0];
        elif coeff(%, z) >= 0 then
                return [P,ans = 0];
        else
                return [P,-ans = 0];
        end if;
end proc:

>	map(tangent_plane, pts);

Download mw.mw

Animate pursuit...

Answered: Rouben Rostamian 8259

August 18 2023

2 1

That's a rather challenging problem. It took me some time to figure it out, and therefore I wrote down the details so that I won't forget.

But...

When I display my worksheet here, it shows a lot of garbage. I don't know why. You may want to download the worksheet from the link below and view it in Maple. This anumation is extracted from the worksheet.

Download pursuit.mw

No solution?...

Answered: Rouben Rostamian 8259

August 14 2023

1 2

It appears that your system of equations has no solution. To see that, plot zz := rhs(eqE)-lhs(EqE) as a function of E1 and E2, with a variety of choices for nu and h1 within the specified ranges. You will see that the zz is always negative, meaning that eqE cannot be satisfied.

zz := unapply((rhs-lhs)(eqE), nu, h__1);

plot3d(zz(0.1, 500), E__1=0..5000, E__2=0..5000);

plot3d(zz(0.3, 500), E__1=0..5000, E__2=0..5000);

plot3d(zz(0.3, 500), E__1=0..5000, E__2=0..5000);

plot3d(zz(0.1, 1000), E__1=0..5000, E__2=0..5000);

plot3d(zz(0.3, 1000), E__1=0..5000, E__2=0..5000);

plot3d(zz(0.5, 1000), E__1=0..5000, E__2=0..5000);

See if this does what you want. ...

Answered: Rouben Rostamian 8259

August 09 2023

0 2

See if this does what you want.

>

restart;

>

w:=(1-delta[h]-delta[b])*(1-phi/(kc-3/2))^(-kc+3/2)*(-kc+3/2)/((kc-3/2)*(1-phi/(kc-3/2)))+delta[h]*(1-phi/(sigma[h]*(kh-3/2)))^(-kh+3/2)*(-kh+3/2)/((kh-3/2)*(1-phi/(sigma[h]*(kh-3/2))))-(1/18)*delta[b]*sqrt(3)*(3*sqrt((M-u+sqrt(3)*sqrt(mu*sigma[b]))^2+2*mu*phi)*mu-3*sqrt((M-u-sqrt(3)*sqrt(mu*sigma[b]))^2+2*mu*phi)*mu)/(mu*sqrt(mu*sigma[b]))-(1/18)*sqrt(3)*(-3*sqrt((M+sqrt(3)*sqrt(sigma[i]))^2-2*phi)+3*sqrt((M-sqrt(3)*sqrt(sigma[i]))^2-2*phi))/sqrt(sigma[i]);

>	e1 := eval([w,diff(w,phi)],phi=0);

>	e2 := simplify(e1) assuming positive;

>	e3 := simplify(e2) assuming positive, M - u - sqrt(3)sqrt(mu)sqrt(sigma[b]) > 0, M - u + sqrt(3)sqrt(mu)sqrt(sigma[b]) > 0, M - sqrt(3)*sqrt(sigma[i]) > 0;

Download New-eval-2.mw

Built-in...

Answered: Rouben Rostamian 8259

August 01 2023

1 3

Maple's OrthogonalExpansions library can calculate the Fourier series for you.

Download mw.mw

A solution...

Answered: Rouben Rostamian 8259

July 31 2023

0 1

The answer to your Question #1 is easy. Here it is. I haven't worked out the details of #2. Perhaps someone else will.

>

restart;

>	with(plots):

Cylinder of radius centered on the axis; parameters and

>	C1 := <p, acos(s), asin(s) >;

Vector(3, {(1) = p, (2) = a*cos(s), (3) = a*sin(s)})

Cylinder of radius centered on the axis; parameters and

>	C2 := <ccos(t), csin(t), q>;

Vector(3, {(1) = c*cos(t), (2) = c*sin(t), (3) = q})

Intersection curve

>	C1 =~ C2; solve(%, {s,p,q}); curve := eval(C2, %); # parametrized by t

Vector(3, {(1) = p = c*cos(t), (2) = a*cos(s) = c*sin(t), (3) = a*sin(s) = q})

${p = c*cos(t), q = a*sqrt(-(c^2*sin(t)^2-a^2)/a^2), s = arccos(c*sin(t)/a)}$

Vector[column](%id = 36893628579494067244)

>

%display(
        %tubeplot([seq](curve), t=-Pi..Pi, radius=0.03, color=yellow),
        %plot3d(C1, s=0..Pi, p=-a..a, color=blue),
        %plot3d(C2, t=-Pi..0, q=0..1.5, color=red, transparency=0.5),
        %plot3d(sqrt(a^2-y^2)+1e-3, x = -sqrt(c^2-y^2) .. sqrt(c^2-y^2), y=-c..c, color=green),
scaling=constrained, style=surface, axes=none, orientation=[20,65,0]):
value(subs(a=1.1, c=1, %));

>

Download intersection-of-cylinders.mw

Elliptic section...

Answered: Rouben Rostamian 8259

July 18 2023

3 0

This worksheet shows how to calculate the parametric equations of the cone, the plane, and the ellipse that lies in their intersection. The cases of parabola and hyperbola can be done in the same way.

Intersection of a plane and a cone

>

restart;

>	with(plots):

Parameters:

: half-angle of the cone's vertex

: angle of the plane's normal relative to the axis (set to get an ellipse)
: coordinate of the plane's intersection with the axis

>	params := beta=0.9*alpha, alpha=Pi/6, d=1/2;

Parametric equation of the cone ( and are the parameters)

>	C := s< tan(alpha)cos(t), tan(alpha)*sin(t), 1>;

Vector(3, {(1) = s*tan(alpha)*cos(t), (2) = s*tan(alpha)*sin(t), (3) = s})

>	p1 := plot3d(subs(params, C), s=-1..1, t=-Pi..Pi, scaling=constrained, color="Red", transparency=0.1);

The vectors form an orthonormal triad, with normal to the
plane and within the plane:

>	N := < sin(beta), 0, cos(beta) >; A := < 0, 1, 0 >; B := LinearAlgebra:-CrossProduct(N,A);

Vector(3, {(1) = sin(beta), (2) = 0, (3) = cos(beta)})

Vector(3, {(1) = 0, (2) = 1, (3) = 0})

Vector[column](%id = 36893628788357808124)

Parametric equation of the plane. ( and are the parameters)

>	P := Au + Bv + d*<0,0,1>;

Vector(3, {(1) = -v*cos(beta), (2) = u, (3) = v*sin(beta)+d})

>	p2 := plot3d(subs(params, P), u=-1..1, v=-1..1, color="Blue", transparency=0.7);

To find the parametric equation of the ellipse , we equate the parametric equations

of the cone and the plane. That gives us a system of three equations in the four unknowns

. We solve the system for in terms of : and then plug the result into the

equation of the cone to get the parametric equation of the ellipse, parametrized with :

>	C -~ P: seq(%): solve({%}, {u,v,s}): E := subs(%, C);

Vector(3, {(1) = d*cos(beta)*tan(alpha)*cos(t)/(sin(beta)*tan(alpha)*cos(t)+cos(beta)), (2) = d*tan(alpha)*sin(t)*cos(beta)/(sin(beta)*tan(alpha)*cos(t)+cos(beta)), (3) = d*cos(beta)/(sin(beta)*tan(alpha)*cos(t)+cos(beta))})

>	p3 := tubeplot([seq](subs(params, E)), t=-Pi..Pi, radius=0.02, color="Yellow");

>	display(p1,p2,p3,style=surface,scaling=constrained, axes=none);

Download elliptic-section.mw

Try this simpler one...

Answered: Rouben Rostamian 8259

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;

(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})

(2)

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

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

(3)

Now is singular:

>	Determinant(M__exact);

(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 is approximately singular,

>	Determinant(M);

(6)

but from the point of view of Digits = 30, 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 ,

apply the singular value decomposition:

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

(8)

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

the last singular value is of the order of , which is almost

zero. That says that the last row of the matrix 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 is in the null space.

Download singular_values.mw

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

8259 Reputation

23 Badges

MaplePrimes Activity

These are answers submitted by Rouben Rostamian

Heap size...

Use "Explore"...

Maybe this?...

Here is how...

The problem is ill-posed...

Symbolic solution...

A numerical solution...

Solution...

Animate pursuit...

No solution?...

See if this does what you want. ...

Built-in...

A solution...

Elliptic section...

Try this simpler one...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again