Rouben Rostamian

This maybe useful...

Answered: Rouben Rostamian 7711

April 21 2019

3 0

I can't tell exactly what problem you are attempting to solve but I can sort of guess by looking through your worksheet. It is not clear to me why you discretize time. It would be more natural to look at time as a continuous variable. With that in mind, the calculations in the following worksheet may be relevant to your objectives.

A homogenous disk of radius and mass rolls without slipping
on a horizontal line. One end of a massless rod of length is attached
to the disk's center. A point of mass is attached to the other end
of the rod. We derive the equations of motion and animate.

>

restart;

>	with(plots): with(plottools):

The position vector to the center of the disk (radius=):

>	r[1] := x -> < x, R >;

The position vector to the pendulum's bob.

Here is the length of the rod and is the angle

of the rod with respect to the vertical.

>	r[1](x) + a*<sin(phi), -cos(phi)>; r[2] := unapply(%, [x,phi]):

Vector(2, {(1) = x+a*sin(phi), (2) = R-a*cos(phi)})

The velocities

>	v[1] := diff(r[1](x(t)), t);

Vector(2, {(1) = diff(x(t), t), (2) = 0})

>	v[2] := diff(r[2](x(t),phi(t)), t);

Vector(2, {(1) = diff(x(t), t)+a*(diff(phi(t), t))*cos(phi(t)), (2) = a*(diff(phi(t), t))*sin(phi(t))})

Assuming that the disk does not slip against the axis leads to this constraint:

>	de1 := diff(x(t), t) = -(diff(phi(t), t))*R;

Now we calculate the system's kinetic energy.
The mass of the disk is . The mass of the bob is .
The moment of inertia of a uniform circular disk of
mass is .

>	1/2M(v[1]^+ . v[1]) + 1/2Jdiff(phi(t),t)^2 + 1/2m(v[2]^+ . v[2]): eval(%, {J=1/2MR^2, de1}): T := simplify(%);

(1/4)*(-4*cos(phi(t))*R*a*m+(3*M+2*m)*R^2+2*a^2*m)*(diff(phi(t), t))^2

and the potential energy

>	V := mgr[2](x(t),phi(t))[2];

>	L := T - V; # the Lagrangian

(1/4)*(-4*cos(phi(t))*R*a*m+(3*M+2*m)*R^2+2*a^2*m)*(diff(phi(t), t))^2-m*g*(R-a*cos(phi(t)))

Apply the Euler-Lagrange formula to obtain the differential equation of motion

>	VariationalCalculus:-EulerLagrange(L, t, phi(t)): remove(type, %, equation): de2 := solve(%, diff(phi(t),t,t))[];

diff(diff(phi(t), t), t) = -2*sin(phi(t))*a*m*((diff(phi(t), t))^2*R+g)/(2*a^2*m-4*cos(phi(t))*R*a*m+3*M*R^2+2*R^2*m)

>	params := M=2, m=1, g=1, R=1, a=3:

>	ic := x(0)=0, phi(0)=0, D(phi)(0)=1;

Solve the system of two differential equations de1, de2 for the unknowns , .

>	dsol := dsolve(eval({de1,de2,ic}, {params}), numeric, output=operator);

>	odeplot(dsol, [t,phi(t)], t=0..20);

The oscillation period:

>	tmax := fsolve(eval(phi(t),dsol), t=10..20);

This proc produces a snapshot of the system at time

We apply the neutral operators %disk, etc., in order to prescribe

the shape of the object before applying the numerical values in params.

>

frame := proc(t)
  local A, B, Disk, Hub, Rod, Bob;
  A := convert(r[1](x(t)), list);
  B := convert(r[2](x(t),phi(t)), list);
  Disk := %disk(A, R, color="Wheat");
  Hub := %pointplot(A, symbol=solidcircle, symbolsize=30, color=black),
          %pointplot(A, symbol=solidcircle, symbolsize=20, color=white);
  Rod := %plot([A,B], thickness=8, color=black);
  Bob := %pointplot(B, symbol=solidcircle, symbolsize=50, color=red);
  eval([Disk, Rod, Bob, Hub], {dsol[], params});
  value(%);
  display(%);
end proc:

Produce a sequence of frames

>	frames := seq(frame(t), t=0..tmax, 0.15):

Animate

>	display([frames], insequence, scaling=constrained, tickmarks=[0,0]);

Download compound-pendulum.mw

Two equations, one unknown...

Answered: Rouben Rostamian 7711

April 15 2019

0 3

You are attempting to solve two equations in the one unknown f. I don't think that pdsolve is the right tool for such a thing.

It just happens that the two equations are compatible and you can solve them by hand if you do things carefully. Here are the steps.

>

restart;

Writing x and xp for the variables is a good way of confusing

yourself. I'd rather use x and y.

>	de1 := diff(f(x,y),x)=-(1/2)(Ly+2x)kl; de2 := diff(f(x,y),y)=-(1/4)klL(Ly+2*x);

Integrate de1 with respect to x. The integration "constant"

can be any function of y:

>	tmp1 := int(lhs(de1), x) = int(rhs(de1),x) + g(y);

Differentiate the result, and then subtract from de2:

>	diff(tmp1, y); tmp2 := simplify(de2 - %);

This is an ODE in g(y). Solve it:

>	tmp3 := dsolve(tmp2);

Plug this back into tmp1:

>	eval(tmp1, tmp3);

This is the function you are looking for, where _C1 is an

arbitrary constant.

Download mw.mw

Solution by recursion...

Answered: Rouben Rostamian 7711

April 12 2019

1 0

In addition to Kitonum's solution, you may want to also look at the following recursive solution:

restart;
u := n -> sin(u(n-1));
u(0) := 1.5;

Then:

seq(u(n), n=0 .. 8);

Here it is...

Answered: Rouben Rostamian 7711

April 11 2019

2 0

Here is how it's done but check to be sure that what I have typed agrees with yours.

Download mw.mw

Need numeric values...

Answered: Rouben Rostamian 7711

April 09 2019

2 1

You have lambda, rho and p in your equations. To solve numerically, you need to assign numeric values to those.

Try this...

Answered: Rouben Rostamian 7711

April 08 2019

0 2

P := seq([a,a^2], a=1..10, 0.5);
plot([P], style=point);

Solving on 0 .. infinity...

Answered: Rouben Rostamian 7711

April 08 2019

3 1

As an alternative to tomleslie's suggestion, here I show how to solve the differential equation numerically on the unbounded interval (0,∞).

>

restart;

>	plots:-setoptions(gridlines=false); # for posting on MaplePrimes

The exact solution

We calculate and plot the exact solution.

>	de := diff(theta(eta),eta,eta) + 2etadiff(theta(eta),eta) = 0;

diff(diff(theta(eta), eta), eta)+2*eta*(diff(theta(eta), eta)) = 0

>	bc := theta(0)=1, theta(infinity)=0;

>	sol_exact := dsolve({de,bc});

>	p1 := plot(rhs(sol_exact), eta=0..infinity, color="Gold", thickness=11, labels=[eta, theta]);

>

The numerical solution

Now we calculate the same solution numerically.

Toward that end, we introduce the change of variable

which maps the infinite interval to the finite interval .

Let . The differential equation changes to:

>	DE_tmp := PDEtools:-dchange( {eta=tan(Pi/2*x), theta(eta)=y(x)}, de, [y(x), x]);

-4*((diff(y(x), x))*sin((1/2)*Pi*x)*Pi-(diff(diff(y(x), x), x))*cos((1/2)*Pi*x))*cos((1/2)*Pi*x)^3/Pi^2+4*tan((1/2)*Pi*x)*(diff(y(x), x))*cos((1/2)*Pi*x)^2/Pi = 0

We simplify the differential equation obtained above.

This is not strictly necessary -- we do it to show that the new

differential equation is actually quite simple, but the simplification

is immaterial to the numeric solver.

>	isolate(DE_tmp, diff(y(x),x,x)): simplify(%): DE := convert(%, tan);

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

The boundary conditions change to . Here is the numeric solution:

>	dsol := dsolve({DE, y(0)=1, y(1)=0}, numeric, output=listprocedure);

We extract the solution :

>	tmp1 := eval(y(x), dsol);

Here is what the solution looks like:

>	plot(tmp1(x), x=0..1, labels=[x,y]);

We find the inverse transformation :

>	isolate(eta=tan(Pi/2*x), x): tmp2 := unapply(rhs(%), eta);

Convert the solution to

>	my_theta := tmp1@tmp2;

Here is what the numerical solution looks like:

>	p2 := plot(my_theta(eta), eta=0..infinity, color=red, labels=[eta,theta]);

>

Comparison of the exact and numerical solutions

We plot the exact and numeric solutions together.

>	plots:-display([p1,p2]);

Download mw.mw

Yes, there is. &n...

Answered: Rouben Rostamian 7711

March 30 2019

0 7

Yes, there is.

Default Zoom...

Answered: Rouben Rostamian 7711

March 28 2019

0 3

Maple's setting are saved in the file $HOME/.maple/2018/maplerc in Linux. Perhaps there are the equivalent files under other operating systems, but I wouldn't know that.

Anyway, in my maplerc I have:

Default Zoom=150

That sets a worksheet's default zoom. I suppose the 150 could be changed to anything else if I wanted to.

One line...

Answered: Rouben Rostamian 7711

March 28 2019

1 0

plot([seq([n*(n+1),n], n=2..16)], style=point);

The general solution...

Answered: Rouben Rostamian 7711

March 27 2019

1 2

The solution obtained by Maple pdsolve (see vv's answer) corresponds to the initial condition u(x,t) = 2*A*cosh(c*x). Here I will show what the true general solution looks like. The underlying mathematical theory is explained in Wikipedia's page on Heat Equation.

>

restart;

We wish to obtain the general solution of

>	pde := diff(u(x,t),t) - K*diff(u(x,t),x,x);

on , with the boundary condition for all ,

and the initial condition for , where

is a given function.

The trick for solving the problem is to extend to the entire

real line as an even function, that is, for

solve the original PDE on the entire real line, and then throw

away the part of the solution corresponding to the negative .

The solution inherits the evenness of the initial condition, and

therefore it automatically satisfies the boundary condition at zero.

The fundamental solution of the heat equation (see Wikipedia) is:

>	G := (x,t) -> 1/sqrt(4PiKt)exp(-x^2/(4Kt));

proc (x, t) options operator, arrow; exp(-(1/4)*x^2/(K*t))/sqrt(4*Pi*K*t) end proc

The general solution of the PDE on the entire real line and
initial condition is obtained through a convolution:

>	u(x,t) = Int(G(x-y,t)*f(y), y=-infinity..infinity);

u(x, t) = Int((1/2)*exp(-(1/4)*(x-y)^2/(K*t))*f(y)/(Pi*K*t)^(1/2), y = -infinity .. infinity)

We split the integral into halves on the negative and positive axis,

apply the symmetry of to simplify the result, and arrive at

>	sol := u(x,t) = Int( (G(x-y,t) + G(x+y,t))*f(y), y=0..infinity);

u(x, t) = Int(((1/2)*exp(-(1/4)*(x-y)^2/(K*t))/(Pi*K*t)^(1/2)+(1/2)*exp(-(1/4)*(x+y)^2/(K*t))/(Pi*K*t)^(1/2))*f(y), y = 0 .. infinity)

Let us verify that the expression obtained above satisfies the PDE:

>	simplify(eval(pde, sol));

and satisfies the boundary condition:

>	diff(sol,x): eval(%, x=0): simplify(%);

eval(diff(u(x, t), x), {x = 0}) = 0

With some effort, it is also possible to show that the solution obtained above satisfies the initial condition in the sense that the limit of u(x,t) as t goes to zero from the right is f(x). I can do that by hand but I don't know how to do it in Maple. Perhaps someone else can.

Download general-solution.mw

Try this...

Answered: Rouben Rostamian 7711

March 25 2019

1 1

restart;
with(plots):
plt := inequal({sqrt((1/2)*x) < y and sqrt(Pi/(2*x)) < y and y < sqrt(x) and y < sqrt(Pi/x)}, x = 1 .. 3, y = 0 .. 2);
display([plt,
textplot([
[1.8,0.7, typeset('sqrt(x/2)')],
[2.5,0.6, typeset('sqrt(Pi/x^2)')]
])
]);

Square brackets are not parentheses!...

Answered: Rouben Rostamian 7711

March 23 2019

1 1

In several places in your code you have used square brackets [ ] for grouping. But square brackaets are not the same as parentheses ( ) in Maple or any other programming language. Here is the fixed code:

restart;
ddesys := {diff(pred(t), t) = pred(t)*(prey(t)/(prey(t)+10)-2/3), diff(prey(t), t) = prey(t)*(2-(1/20)*prey(t-tau))-prey(t)*(pred(t)/(prey(t)+10))-10, pred(0) = 1, prey(0) = 1};
dsn := dsolve(eval(ddesys, tau = 0), numeric);
plots[odeplot](dsn, [[t, prey(t), color = red], [t, pred(t), color = blue]], t=0 .. 0.3, legend = [prey, pred], labels = [t, ""]);

In the last line I have replaced your plotting range of 0..300 with 0..0.3. When you execute the code, you will see that the prey population drops down to zero and then becomes negative at time t=0.1. This indicates a flaw in your model. See how you can remedy that.

Solution via the Lagrangian formulation...

Answered: Rouben Rostamian 7711

March 22 2019

4 1

Your solution seems to be correct in general outline although I have not checked the details. Part of the complication in what you have presented is due to the need to handle the internal forces within the rod.

Lagrangian mechanics (analytical mechanics) introduces a method for deriving a system's equations of motion without the need to deal with the internal forces, and that significantly simplifies the analysis. Here I illustrate how to solve your problem through the Lagrangian formulation.

>

restart;

>	with(plots):

and : position vectors of the two balls.
: the angle of the rod measured counterclockwise

from the downward vertical direction
: length of the rod

>	< x(t), y(t) >; r[1] := unapply(%, t): r[1](t) + a*< sin(theta(t)), -cos(theta(t))>; r[2] := unapply(%, t):

Vector(2, {(1) = x(t), (2) = y(t)})

Vector[column](%id = 18446884489083941222)

Velocities of the balls

>	v[1] := diff~(r[1](t),t); v[2] := diff~(r[2](t),t);

Vector(2, {(1) = diff(x(t), t), (2) = diff(y(t), t)})

Vector[column](%id = 18446884489083937366)

Kinetic energy of the system

>	T := 1/2 * m[1] * (v[1]^+ . v[1]) + 1/2 * m[2] * (v[2]^+ . v[2]);

(1/2)*m[1]*((diff(x(t), t))^2+(diff(y(t), t))^2)+(1/2)*m[2]*((diff(x(t), t)+a*(diff(theta(t), t))*cos(theta(t)))^2+(diff(y(t), t)+a*(diff(theta(t), t))*sin(theta(t)))^2)

Potential energy of the system

>	V := m[1] * g * r[1](t)[2] + m[2] * g * r[2](t)[2];

The Lagrangian

>	L := T - V;

(1/2)*m[1]*((diff(x(t), t))^2+(diff(y(t), t))^2)+(1/2)*m[2]*((diff(x(t), t)+a*(diff(theta(t), t))*cos(theta(t)))^2+(diff(y(t), t)+a*(diff(theta(t), t))*sin(theta(t)))^2)-g*m[1]*y(t)-g*m[2]*(y(t)-a*cos(theta(t)))

The equations of motion and integrals of motion via Euler-Lagrange:

>	tmp := VariationalCalculus:-EulerLagrange(L, t, [x(t),y(t),theta(t)]):

The integrals of motion that appear above as something= and something=
are not of interest to us, so we get rid of them

>	EL := remove(type, tmp, 'equation'):

What remains are the equations of motion

>	de1 := simplify(EL[1]) = 0;

(-m[1]-m[2])*(diff(diff(x(t), t), t))+a*m[2]*((diff(theta(t), t))^2*sin(theta(t))-cos(theta(t))*(diff(diff(theta(t), t), t))) = 0

>	de2 := simplify(EL[2]) = 0;

-a*m[2]*((diff(diff(theta(t), t), t))*a+(diff(diff(x(t), t), t))*cos(theta(t))+sin(theta(t))*(diff(diff(y(t), t), t)+g)) = 0

>	de3 := simplify(EL[3]) = 0;

(-m[1]-m[2])*(diff(diff(y(t), t), t))-(diff(diff(theta(t), t), t))*sin(theta(t))*a*m[2]-cos(theta(t))*(diff(theta(t), t))^2*a*m[2]-g*(m[1]+m[2]) = 0

Pick parameters:

>	params := a=5, m[1]=1, m[2]=2, g=9.81:

Pick initial conditions:

>	ic := x(0)=0, y(0)=0, theta(0)=Pi/2, D(x)(0)=10, D(y)(0)=0, D(theta)(0)=2*Pi:

Here is the complete specification of the initial value problem:

>	sys := eval({de1,de2,de3,ic}, {params}):

Solve the initial value problem:

>	dsol := dsolve(sys, numeric, output=operator):

Find the time when ball #1 hits the ground:

>	eval(y(t),dsol): tmax := fsolve(%, t=1e-5..infinity);

Plot the paths of the two balls.

Note: In plotting the blue path we need "params" to get the

value of the bar's length.

>	paths := plot([ eval([r[1](t)[1], r[1](t)[2], t=0..tmax], dsol), eval([r[2](t)[1], r[2](t)[2], t=0..tmax], {dsol[],params}) ], color=["Pink","SkyBlue"]);

Produce a snapshot of the rod at time

>

frame := proc(t)
  pointplot(eval(r[1](t), dsol),
      symbol=solidcircle, symbolsize=20, color=red);
  pointplot(eval(r[2](t), {dsol[], params}),
      symbol=solidcircle, symbolsize=20, color=blue);
  pointplot(eval([r[1](t),r[2](t)], {dsol[], params}),
connect, color="Green");
  display([%%%,%%,%]);
end proc:

Produce animation:

>	seq(frame(t), t=0..tmax, 0.1): display(%, insequence): display([%, paths], scaling=constrained);

Download thrown-bola.mw

Is this good?...

Answered: Rouben Rostamian 7711

March 13 2019

0 0

There is no expression for a "general solution" to your equation. Normally you would want to solve a PDE along with some boundary conditions which you have not specified. If you are happy with just any particular solution, consider this one:

mysol := C0 + C1*r*cos(theta) + C2*r*sin(theta);

To verify that this is indeed a solution for any three constants C0, C1, and C2, do

simplify(subs(a(r,theta)=mysol, pde));

where "pde" stands for your PDE.

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

7711 Reputation

22 Badges

MaplePrimes Activity

These are answers submitted by Rouben Rostamian

This maybe useful...

Two equations, one unknown...

Solution by recursion...

Here it is...

Need numeric values...

Try this...

Solving on 0 .. infinity...

Yes, there is. &n...

Default Zoom...

One line...

The general solution...

Try this...

Square brackets are not parentheses!...

Solution via the Lagrangian formulation...

Is this good?...

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