Rouben Rostamian

Answer...

Answered: Rouben Rostamian 8596

June 27 2020

3 22

>

restart;

>	pde := diff(w(x,t), x$4) + diff(w(x,t), t$2)/c^2 = 0;

diff(diff(diff(diff(w(x, t), x), x), x), x)+(diff(diff(w(x, t), t), t))/c^2 = 0

We look for solutions of the form :

>	eval(pde, w(x,t)=u(x)cos(omegat)); de_tmp := expand(%/cos(omega*t));

(diff(diff(diff(diff(u(x), x), x), x), x))*cos(omega*t)-u(x)*omega^2*cos(omega*t)/c^2 = 0

diff(diff(diff(diff(u(x), x), x), x), x)-u(x)*omega^2/c^2 = 0

The coefficient is not convenient. It's better to work with dimensionless quantities.

Therefore we introduce a new unknown, , related to through

>	mu_def := omega = c/L^2*mu^2;

omega = c*mu^2/L^2

Then express the differential equation de_tmp in terms of

>	de := eval(de_tmp, mu_def);

diff(diff(diff(diff(u(x), x), x), x), x)-u(x)*mu^4/L^4 = 0

Here is the general solution of the de:

>	dsolve(de): U := unapply(rhs(%), x);

proc (x) options operator, arrow; _C1*exp(mu*x/L)+_C2*exp(-mu*x/L)+_C3*sin(mu*x/L)+_C4*cos(mu*x/L) end proc

We want to apply the boundary conditions , but of course
Maple will pick the trivial solution . The get around that, we temporarily suspend the
boundary condition and instead substitute a boundary condition , where
is an unspecified symbol. Applying the four boundary conditions, we obtain four equations in

the four unknowns which we solve:

>	U(0) = 0, D(U)(0)=0, (D@@2)(U)(0) = a, U(L) = 0: C_vals := solve({%}, {_C1, _C2, _C3, _C4});

${_C1 = -(1/2)*a*L^2*(-cos(mu)+sin(mu)+exp(-mu))/(mu^2*(exp(mu)-exp(-mu)-2*sin(mu))), _C2 = (1/2)*a*L^2*(exp(mu)-cos(mu)-sin(mu))/((exp(mu)-exp(-mu)-2*sin(mu))*mu^2), _C3 = (1/2)*a*L^2*(exp(mu)-2*cos(mu)+exp(-mu))/(mu^2*(exp(mu)-exp(-mu)-2*sin(mu))), _C4 = -(1/2)*a*L^2/mu^2}$

Plug these into the expression for , compute the derivative, and set the derivative to zero
in order to assert the boundary condition which we had left out.

>	eval(U(x), C_vals): diff(%, x): eval(%, x=L) = 0: simplify(%);

a*L*(exp(2*mu)*cos(mu)-2*exp(mu)+cos(mu))/(mu*(exp(2*mu)-2*sin(mu)*exp(mu)-1)) = 0

The denominator of that equation is a multiple of and therefore nonzero.
We conclude that

>	exp(2mu)cos(mu) - 2*exp(mu) + cos(mu) = 0; eq_tmp := isolate(%, cos(mu));

cos(mu) = 2*exp(mu)/(exp(2*mu)+1)

Let's verify that the right-hand side of this equation is

>	simplify(1/cosh(mu) - rhs(eq_tmp));

Okay, therefore eq_tmp is equivalent to:

>	eq := cos(mu) = 1/cosh(mu);

This equation has infinitely many roots, as can be seen in this graph:

>	plot([cos(mu),1/cosh(mu)], mu=0..5*Pi, color=[blue,red], thickness=3, size=[700,300]);

Here are the first few roots:

>	seq(fsolve(cos(mu)=1/cosh(mu), mu=(2n+1)Pi/2), n=1..5);

Call the roots .. For each root there is a mode function given by

>	subs(C_vals, mu=mu[j], U(x));

-(1/2)*a*L^2*(-cos(mu[j])+sin(mu[j])+exp(-mu[j]))*exp(mu[j]*x/L)/(mu[j]^2*(exp(mu[j])-exp(-mu[j])-2*sin(mu[j])))+(1/2)*a*L^2*(exp(mu[j])-cos(mu[j])-sin(mu[j]))*exp(-mu[j]*x/L)/((exp(mu[j])-exp(-mu[j])-2*sin(mu[j]))*mu[j]^2)+(1/2)*a*L^2*(exp(mu[j])-2*cos(mu[j])+exp(-mu[j]))*sin(mu[j]*x/L)/(mu[j]^2*(exp(mu[j])-exp(-mu[j])-2*sin(mu[j])))-(1/2)*a*L^2*cos(mu[j]*x/L)/mu[j]^2

where the factors may be dropped since a mode function is determined modulo a

multiplicative constant. Also for each there corresponds a frequency given through

the equation mu_def, that is, "`omega__j`=(c*`mu__j`^(2))/(L^(2))." You should be all set now for eigenfunction expansion

of the solutions of your PDE.

>

Download mw.mw

Solution...

Answered: Rouben Rostamian 8596

June 27 2020

4 1

>

restart;

>	d := 6/7; # the delay

Plot the history function

>	p[0] := plot(cos(t), t=-d..0):

The solution on the interval picks up the history function (colored green)

>	diff(x(t),t,t) + 2/25diff(x(t),t) + 4x(t) = (1/25diff(x(t),t,t) + 2/3diff(x(t),t) + 5/2x(t))cos(t-d), x(0)=1, D(x)(0)=0: dsol := dsolve({%}, numeric, output=operator, range=0..d): myx[1] := eval(x, dsol): mydx[1] := eval(D(x), dsol): p[1] := plot(myx[1](t), t=0..d):

On each subsequent interval , the differential equation picks up the solution

calculated in the preceding interval (colored green)

>

for i from 2 to 14 do
sys := diff(x(t),t,t) + 2/25*diff(x(t),t) + 4*x(t)
= (1/25*diff(x(t),t,t) + 2/3*diff(x(t),t) + 5/2*x(t))*myx[i-1](t-d),
x((i-1)*d)=myx[i-1]((i-1)*d), D(x)((i-1)*d)=mydx[i-1]((i-1)*d):
dsol := dsolve({sys}, numeric, output=operator, range=(i-1)*d..i*d):
myx[i] := eval(x, dsol):
mydx[i] := eval(D(x), dsol):
p[i] := plot(myx[i](t), t=(i-1)*d..i*d);
end do:
plots:-display(seq(p[j], j=0..i-1));

Download mw.mw

Solution...

Answered: Rouben Rostamian 8596

June 25 2020

2 4

Here is the solution with no explanations. See if you can figure it out.

>

restart;

>	C1 := x^2 + y^2 = L^2; C2 := (x-R)^2 + y^2 = R^2;

>	solve(C1, y): c1 := %[1];

>	solve(C2, y): c2 := %[1];

>	X := solve(c1=c2, x);

>	A1 := int(c1, x=X..L) assuming L>0, 2*R>L;

-(1/8)*L^2*(4*arcsin((1/2)*L/R)*R^2-2*Pi*R^2+L*(-L^2+4*R^2)^(1/2))/R^2

>	A2 := int(c2, x=0..X) assuming L>0, 2*R>L;

(1/8)*(2*R^4*Pi+4*R^4*arcsin((1/2)*(L^2-2*R^2)/R^2)+L^3*(-L^2+4*R^2)^(1/2)-2*L*(-L^2+4*R^2)^(1/2)*R^2)/R^2

>	A := simplify(2*(A1+A2));

(1/2)*(-2*arcsin((1/2)*L/R)+Pi)*L^2-(1/2)*L*(-L^2+4*R^2)^(1/2)+(arcsin((1/2)*(L^2-2*R^2)/R^2)+(1/2)*Pi)*R^2

>	A = Pi*R^2/2: eval(%, R=1): fsolve(%, L=0..2);

Does not exist...

Answered: Rouben Rostamian 8596

June 23 2020

0 1

Such a rational function does not exist! Is that a trick question?

Corrections...

Answered: Rouben Rostamian 8596

June 22 2020

4 3

Since you are not in the habit of posting a worksheet, I won't post a worksheet either. I will give my comments as text.

Make the following corrections in your code:

The mathematical π is written Pi (with uppercase P) in Maple.
The gamma function is written GAMMA (all uppercase) in Maple.
Call your differential equations de1 and de2.
In the line that defines M1, change "=" to ":=".
Enter the initial conditions as
ic := U1(0)=M1, phi(0)=0, D(phi)(0)=0.001;
Set T:=0.1; beta:=0.6; k:=3.5;

After these changes you should be able to call Maple's dsolve to solve your system numerically:

dsol := dsolve({de1,de2,ic}, numeric);

To plot U1, phi, and the derivative of phi over the time interval 0 to 30, do

plots:-odeplot(dsol, [x,U1(x)], x=0..30, view=0..1);
plots:-odeplot(dsol, [x,phi(x)], x=0..30);
plots:-odeplot(dsol, [x,diff(phi(x),x)], x=0..30);

You will see that Maple issues a warning, saying that it cannot go beyond 23.53. If you look closely for the reason for that, you will see that at t=23.53 the value of phi is zero and the value of its derivative is negative, therefore at the next timestep the value of phi is going to turn negative. But that is problematic since your differential equation contains the 3/2 power of phi. What is the 3/2 power of a negative number? Is there a reason why you have that term in your equation? Perhaps it shouldn't be there. Or perhaps it should be the absolute value of phi raised to the 3/2 power. Go back to the origin of the equation and investigate.

The longitudinal motion of a rod...

Answered: Rouben Rostamian 8596

June 16 2020

3 1

Here is my paraphrasing of the problem's statement.

The longitudinal vibration of a homogeneous one-dimensional rod of fixed cross-sectional area , modulus of elasticity , and density (mass per unit volume) of , satisfies the PDE where is the displacement of the point at at time , and the subscripts and indicate derivatives. It is customary to let , and write the PDE in the form of the standard wave equation . The strain at any cross-section is , and therefore the stress is , and the force acting on that cross-section is .

Consider such a rod of length , one end of which is fixed at , and the other end is free to move along the axis. The rod's mass is We attach a block of mass to the moving tip of the rod and set it into motion. We wish to model the rod's motion. The initial conditions are and
"((&PartialD; u)/(&PartialD; t)) ? ()|() ? (t=0)={[[0, 0<=x<L],[V,x=L]]"

This is a somewhat tricky problem since the initial velocity is discontinuous at the boundary. In the attached worksheet I show that the displacement of the rod is given by
u(x, t) = r*L*V*(sum(K__j^2*sin(`μ__j`)*sin(`μ__j`*x/L)*sin(c*`μ__j`*t/L)/`μ__j`, j = 1 .. infinity))/c ,
where K__j^2 = 1/(r*sin(`μ__j`)^2+1/2*(1-sin(2*mu[j])/(2*mu[j]))) ,
and the are the solutions of the equation .
Here is a sample of what the solution looks like (u(x,t) versus x, at various times t):

and here is an animation of the rod itself:

The short vertical lines are equally-spaced markers etched on the rod to help visualize the rod's motion.

Download worksheet: longitudinal-motion.mw

Proof: Use v(i) instead of v[i]...

Answered: Rouben Rostamian 8596

June 14 2020

3 0

>

restart;

>	T := v -> sum(v(i)*x^(i-1),i=1..n);

proc (v) options operator, arrow; sum(v(i)*x^(i-1), i = 1 .. n) end proc

We want to show that , and for all
-tuples and , and all scalars . The first one is easy:

>	combine(T(x)+T(y)) = T(x+y); is(%);

sum(x(i)*x^(i-1)+y(i)*x^(i-1), i = 1 .. n) = sum((y(i)+x(i))*x^(i-1), i = 1 .. n)

The second one takes a bit more of effort since we need to tell Maple that
is an -tuple while is a scalar. We have

>

T(c*x);

sum(c(i)*x(i)*x^(i-1), i = 1 .. n)

which is not what we want since should be just . To fix that, we define
to be the constant function that assigns the constant c to any , that is,

>	eval(T(c*x), c=(t->c));

sum(c*x(i)*x^(i-1), i = 1 .. n)

Now that's better. Let's verify that :

>	eval(T(cx), c=(t->c)) = cT(x); is(expand(%));

sum(c*x(i)*x^(i-1), i = 1 .. n) = c*(sum(x(i)*x^(i-1), i = 1 .. n))

QED

Here is a phase portrait...

Answered: Rouben Rostamian 8596

June 08 2020

3 8

>

restart;

Don't begin with the differential equations. Define their right-hand sides first:

>	f[0] := (x,y1,y2) -> alphax(1 - x/N) - beta[1]sqrt(x)y1/(1 + h[1]beta[1]sqrt(x)) - beta[2]sqrt(x)y2/(1 + h[2]beta[2]sqrt(x)) - dEx;

proc (x, y1, y2) options operator, arrow; alpha*x*(1-x/N)-beta[1]*sqrt(x)*y1/(1+h[1]*beta[1]*sqrt(x))-beta[2]*sqrt(x)*y2/(1+h[2]*beta[2]*sqrt(x))-d*E*x end proc

>	f[1] := (x,y1,y2) -> -omega[1]y1 + Mu[1]beta[1]sqrt(x)y1(1 - y1/(beta[1]sqrt(x)))/(1 + h[1]beta[1]sqrt(x));

proc (x, y1, y2) options operator, arrow; -omega[1]*y1+Mu[1]*beta[1]*sqrt(x)*y1*(1-y1/(beta[1]*sqrt(x)))/(1+h[1]*beta[1]*sqrt(x)) end proc

>	f[2] := (x,y1,y2) -> -omega[2]y2 + Mu[2]beta[2]sqrt(x)y2(1 - y2/(beta[2]sqrt(x)))/(1 + h[2]beta[2]sqrt(x));

proc (x, y1, y2) options operator, arrow; -omega[2]*y2+Mu[2]*beta[2]*sqrt(x)*y2*(1-y2/(beta[2]*sqrt(x)))/(1+h[2]*beta[2]*sqrt(x)) end proc

Now define the differential equations:

>	de[0] := diff(x(t),t) = f[0](x(t),y1(t),y2(t));

diff(x(t), t) = alpha*x(t)*(1-x(t)/N)-beta[1]*x(t)^(1/2)*y1(t)/(1+h[1]*beta[1]*x(t)^(1/2))-beta[2]*x(t)^(1/2)*y2(t)/(1+h[2]*beta[2]*x(t)^(1/2))-d*E*x(t)

>	de[1] := diff(y1(t),t) = f[1](x(t),y1(t),y2(t));

diff(y1(t), t) = -omega[1]*y1(t)+Mu[1]*beta[1]*x(t)^(1/2)*y1(t)*(1-y1(t)/(beta[1]*x(t)^(1/2)))/(1+h[1]*beta[1]*x(t)^(1/2))

>	de[2] := diff(y2(t),t) = f[2](x(t),y1(t),y2(t));

diff(y2(t), t) = -omega[2]*y2(t)+Mu[2]*beta[2]*x(t)^(1/2)*y2(t)*(1-y2(t)/(beta[2]*x(t)^(1/2)))/(1+h[2]*beta[2]*x(t)^(1/2))

Numbers come last:

>	params := alpha = 0.75, omega[1] = 0.15, omega[2] = 0.20, N = 0.8, beta[1] = 0.65, beta[2] = 0.70, Mu[1] = 0.75, Mu[2] = 0.74, h[1] = 0.005, h[2] = 0.004, d = 0.05, E = 0.35;

Here is the system:

>	sys := eval({de[0], de[1],de[2]}, {params});

${diff(x(t), t) = .75*x(t)*(1-1.250000000*x(t))-.65*x(t)^(1/2)*y1(t)/(1+0.325e-2*x(t)^(1/2))-.70*x(t)^(1/2)*y2(t)/(1+0.280e-2*x(t)^(1/2))-0.175e-1*x(t), diff(y1(t), t) = -.15*y1(t)+.4875*x(t)^(1/2)*y1(t)*(1-1.538461538*y1(t)/x(t)^(1/2))/(1+0.325e-2*x(t)^(1/2)), diff(y2(t), t) = -.20*y2(t)+.5180*x(t)^(1/2)*y2(t)*(1-1.428571429*y2(t)/x(t)^(1/2))/(1+0.280e-2*x(t)^(1/2))}$

Need some trial-and-error, and some insight, to find a good range for plotting.

>	xmin, xmax := 0.1, 0.6; y1min, y1max := 0.0, 0.6; y2min, y2max := 0.1, 0.2;

Select a good set of initial conditions:

>

ic := seq([x(0)=xmax, y1(0)=y1max, y2(0)=(1-i/5)*y2min+i/5*y2max], i=0..5),
      seq([x(0)=xmax, y1(0)=0.5*y1max, y2(0)=(1-i/5)*y2min+i/5*y2max], i=0..5),
      seq([x(0)=0.5*xmax, y1(0)=y1max, y2(0)=(1-i/5)*y2min+i/5*y2max], i=0..5),
      seq([x(0)=xmax, y1(0)=0.1, y2(0)=(1-i/5)*y2min+i/5*y2max], i=0..5):

Now plot:

>	DEtools:-DEplot3d(sys, [x(t),y1(t),y2(t)], t=0..40, stepsize=0.1, [ic], x=xmin..xmax, y1=y1min..y1max, y2=y2min..y2max, arrows=cheap, linecolor=black, thickness=1, orientation=[30,72,0]);

We see that there is a stable equilibrium near x=0.3, y1=0.2, y2=0.15. Let's find it:

>	eval({f[0](x,y1,y2),f[1](x,y1,y2),f[2](x,y1,y2)}, {params}): fsolve(%, {x=0.3, y1=0.2, y2=0.15});

Download mw.mw

The Method of Lines...

Answered: Rouben Rostamian 8596

May 27 2020

4 2

I will show how to solve the PDE with the Method of Lines. For no particular reason I have changed all your coefficients to "1"s, and have added boundary conditions which are missing in your worksheet.

Additionally, I should point out that you may solve the PDE directly with Maple's pdsolve. In this worksheet I solve and plot the solutions both with pdsolve, and with Method of Lines, and of course the two solutions agree.

>

restart;

1. Solving the PDE with Maple's pdsolve

Here is how to solve the PDE with Maple's pdsolve.

>	pde := diff(u(x,t),t) = - diff(u(x,t),x$2) - diff(u(x,t),x$4) - u(x,t)*diff(u(x,t),x);

>	ic := u(x,0) = sin(Pi*x/L);

>	bc := u(0,t)=0, D[1,1](u)(0,t)=0, u(L,t)=0, D[1,1](u)(L,t)=0;

>	pdsol := pdsolve(pde, eval({ic,bc}, L=1), numeric, spacestep=0.01);

>	pdsol:-plot3d(t=0..0.020);

2. Solving the PDE with the Method of Lines

We divide the interval into equal subintervals and write

for the length of each subinterval. Let be the coordinates

of the endpoints (the nodes) of the intervals . Note that there are
nodes altogether, with nodes and on the boundary, and nodes
through in the interior. In my comments below, I write u_x
for the derivative of with respect to

The finite difference discretizations seen below are pretty standard

in numerical analysis.

>	# discretize u_x, i=1,..,n-1 d1 := add([-1,0,1] ~ [u[i+k](t) $k=-1..1])/(2h);

>	# discretize u_xx, i=1,..,n-1 d2 := add([1,-2,1] *~ [u[i+k](t) $k=-1..1])/h^2;

(u[i-1](t)-2*u[i](t)+u[i+1](t))/h^2

>	# discretize u_xx at i=0 d2_left := add([2,-5,4,-1] *~ [u[k](t) $k=0..3])/h^2;

(2*u[0](t)-5*u[1](t)+4*u[2](t)-u[3](t))/h^2

>	# discretize u_xx at i=n d2_right := add([-1,4,-5,2] *~ [u[n+k](t) $k=-3..0])/h^2;

(-u[n-3](t)+4*u[n-2](t)-5*u[n-1](t)+2*u[n](t))/h^2

>	# discretize u_xxxx, i=2,..,n-2 d4 := add([1,-4,6,-4,1] *~ [u[i+k](t) $k=-2..2])/h^4;

(u[i-2](t)-4*u[i-1](t)+6*u[i](t)-4*u[i+1](t)+u[i+2](t))/h^4

>	# the discretized PDE at any i=2..n-2 eq__i := diff(u[i](t),t) = - d2 - d4 - u[i](t)*d1;

diff(u[i](t), t) = -(u[i-1](t)-2*u[i](t)+u[i+1](t))/h^2-(u[i-2](t)-4*u[i-1](t)+6*u[i](t)-4*u[i+1](t)+u[i+2](t))/h^4-(1/2)*u[i](t)*(-u[i-1](t)+u[i+1](t))/h

Boundary conditions at the left end are u=0 and u_xx=0.
These determine u[0] and u[1] in terms of u[2] and u[3].

>	u[0](t)=0, d2_left = 0: bc_left := solve({%}, {u[0](t), u[1](t)});

Boundary conditions at the right end are u=0 and u_xx=0.
These determine u[n-1] and u[n] in terms of u[n-2] and u[n-3].

>	u[n](t)=0, d2_right = 0: bc_right := solve({%}, {u[n-1](t), u[n](t)});

We pick a reasonable and generate the system of ODEs

>	n := 15: L := 1: h := L/n: T := 0.02: # solution will be computed on the time interval 0..T sys := { seq(eq__i, i=2..n-2) } union bc_left union bc_right:

and here are the corresponding initial conditions

>	IC := {seq(u[i](0)=eval(rhs(ic), x=i*h), i=2..n-2)}:

At this point, the PDE has been converted into a system of ODEs

in unknowns. (That's the purpose of the Methods of Lines.)

We call Maple's dsolve to solve that system.

>	dsol := dsolve(sys union IC, numeric, output=operator, range=0..T, maxfun=0):

Here is the plot of

>	plots:-odeplot(dsol, [t, u[5](t)]);

This proc plots the solution produced by Method of Lines:

>	my_plot3d := proc() local i, j, V, m:=40, k:= T/m; V := Array(0..n, 0..m, datatype=float[8]): for i from 0 to m do V[0..n,i] := Array([eval('u[j](i*k)' $j=0..n, dsol)]); end do; plots:-display(GRID(0..L, 0..T, V), labels=['x','t','u'], _rest); end proc:

Here is what the solution looks like. Compare to that produced by Maple's pdsolve:

>	my_plot3d();

Download method-of-lines-illustration.mw

Another solution...

Answered: Rouben Rostamian 8596

May 24 2020

2 2

Here's an alternative to what Carl wrote. This one avoids implicitplot, and therefore produces crisper graphics.

>

restart;

>	with(plots):

>	with(plottools):

Let's see what the overall picture looks like

>	display([ sphere([0,0,0], 4, color=red), sphere([-4,0,2], 7, color=blue) ], scaling=constrained, transparency=0.6, labels=[x,y,z], orientation=[180,90,-90]);

Here is a 2D cross sectional diagram

>	display([ circle([ 0,0], 4, color=red), circle([-4,2], 7, color=blue), line([0,0], [-4,2], color="Green", thickness=3) ], scaling=constrained);

The green line segment in the diagram above connects the centers of the

circles and makes an angle of with respect to the vertical axis:

>	sigma := arctan(4/2);

If we rotate the above diagram by an angle so that the green line segment

becomes vertical, the diagram simplifies for subsequent calculations:

>	display([ circle([0,0], 4, color=red), circle([0,4], 7, color=blue), line([0,0], [0, sqrt(2^2+4^2)], color="Green", thickness=3), line([0,0], [7sqrt(15)/8, -17/8], color=red, thickness=3), line([0,4], [7sqrt(15)/8, -17/8], color=blue, thickness=3) ], scaling=constrained, labels=[x,r]);

Write down the equations of the two circles, and then find their intersection points:

>	x^2+r^2=16, x^2+(r-4)^2=49; rx := allvalues(solve({%}));

${r = -17/8, x = (7/8)*15^(1/2)}, {r = -17/8, x = -(7/8)*15^(1/2)}$

The angles of the red and blue line segments in the diagram above, relative

to the vertical axis:

>	alpha := Pi/2 + arctan( 17/8,7sqrt(15)/8); beta := Pi/2 + arctan(4+17/8,7sqrt(15)/8);

Now we plot the object in the spherical coordinates:

>	p1 := display([ plot3d([4sin(theta)cos(phi), 4sin(theta)sin(phi), 4cos(theta)], theta=0..alpha, phi=0..2Pi), plot3d([7sin(theta)cos(phi), 7sin(theta)sin(phi), 4+7cos(theta)], theta=beta..Pi, phi=0..2Pi) ], color=["Red","Green"]);

Restore the original orientation by rotating about the axis by

>	plottools:-rotate(p1, sigma, [[0,0,0],[0,1,0]]): display(%, labels=[x,y,z], scaling=constrained, style=surface, orientation=[-125,75,0]);

>

Download mw.mw

Testing for a solution...

Answered: Rouben Rostamian 8596

May 22 2020

2 0

Here is a slightly modified and extended version of your worksheet. The edited parts are in green.

>

restart:

>	with(PDEtools):

>	with(LinearAlgebra):

>	alias(f=f(x,t),g=g(x,t));

>	eq1:=diff(f,x)=-Ietaf +Iexp(-It)*g;

>	eq2:=diff(g,x)=-Ietag +Iexp(It)*f;

>	eq3:=diff(f,t)=(Ieta^2-I/2)f +Ietaexp(-It)g;

>	eq4:=diff(g,t)=(-Ieta^2+I/2)g +Ietaexp(It)f;

>	#### The solution of (2)-(5) is

>	#eq5:=f=I(c1exp(A)-c2exp(-A))exp(-it/2); eq5:=f=I(c1exp(A)-c2exp(-A))exp(-It/2); # changed i to I

>	#eq6:=g=(c2exp(A)-c1exp(-A))exp(it/2); eq6:=g=(c2exp(A)-c1exp(-A))exp(It/2); # changed i to I

>	#### where

>	#c1=sqrt(h-sqrt(h^2-1))/sqrt(h^2-1);c2=sqrt(h+sqrt(h^2-1))/sqrt(h^2-1);A=sqrt(h^2-1)(x+Ih*t);

>	# replaced the three semicolons with commas:

>	vals := c1=sqrt(h-sqrt(h^2-1))/sqrt(h^2-1), c2=sqrt(h+sqrt(h^2-1))/sqrt(h^2-1), A=sqrt(h^2-1)(x+Ih*t);

c1 = (h-(h^2-1)^(1/2))^(1/2)/(h^2-1)^(1/2), c2 = (h+(h^2-1)^(1/2))^(1/2)/(h^2-1)^(1/2), A = (h^2-1)^(1/2)*(x+I*h*t)

>	#### How to verify (6) and (7) is the solution of (2)-(5)?

>	z1 := eval({eq5, eq6}, {vals}): # plug in the values of c1, c2, A

>	# Test the proposed solution. If the test is successful, # we will get a list of zeros.

>	z2 := pdetest(z1, [eq1,eq2,eq3,eq4]);

>	# Those expressions are too complicated to see whether # or not they are zeros. We plug in numbers for x, t, eta, # and h, to see whether at least for these numbers we get # zeros.

>	evalf(eval(z2, {x=0, t=0, eta=0, h=5}));

>	# We see that the result is not a list of zeros. Therefore # the proposed "solution" is not a solution.

>

Download verification-2.mw

See if this helps...

Answered: Rouben Rostamian 8596

May 17 2020

1 0

The general solution...

Answered: Rouben Rostamian 8596

May 07 2020

2 0

Here is how to plot a general ellipsoid with semiaxis lengths a, b, c, centered at x₀, y₀, z₀. The variables t and s are the polar and azimuthal angles in spherical coordinates.

restart;
a, b, c, x0, y0, z0 := 2, 3, 1, 10, 11, 12:
plot3d([x0 + a*sin(t)*cos(s), y0 + b*sin(t)*sin(s), z0 + c*cos(t)],
    t=0..Pi, s=0..2*Pi, scaling=constrained, labels=[x,y,z]);

Tenth degree polynomial...

Answered: Rouben Rostamian 8596

May 06 2020

1 0

We can read off the y values of the function at 9 points x=-4, -3, ... , 4, 5 (except for x=2) on the graph that you have provided. To that collection I added estimated values of the function at -4.25 and 5.25, making for a total of eleven (x,y) pairs. Then we fit the data with a 10th degree polynomial (which has 11 coefficients) and solve for the coefficients.

Note that the graph does not quite look like the given graph. One reason for that is that in the given graph the coordinates are not spaced uniformly. For instance, the point with coordinate −5 on the y axis falls where −3 would be.

That said, I agree with Kitonum's recommendation of using splines in such problems. We are just lucky that a tenth degree polynomial provides a reasonable answer in this particular case. In general, fitting a single polynomial to arbitrary data will not produce good results.

Download mw.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

8596 Reputation

23 Badges

MaplePrimes Activity

These are answers submitted by Rouben Rostamian

Needs a refined initmesh...

Answer...

Solution...

Solution...

Does not exist...

Corrections...

The longitudinal motion of a rod...

Proof: Use v(i) instead of v[i]...

Here is a phase portrait...

The Method of Lines...

Another solution...

Testing for a solution...

See if this helps...

The general solution...

Tenth degree polynomial...

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