Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

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

 

 

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.

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]);

(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]*3^(1/2)*(3*((M-u+3^(1/2)*(mu*sigma[b])^(1/2))^2+2*mu*phi)^(1/2)*mu-3*((M-u-3^(1/2)*(mu*sigma[b])^(1/2))^2+2*mu*phi)^(1/2)*mu)/(mu*(mu*sigma[b])^(1/2))-(1/18)*3^(1/2)*(-3*((M+3^(1/2)*sigma[i]^(1/2))^2-2*phi)^(1/2)+3*((M-3^(1/2)*sigma[i]^(1/2))^2-2*phi)^(1/2))/sigma[i]^(1/2)

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

[(1-delta[h]-delta[b])*(-kc+3/2)/(kc-3/2)+delta[h]*(-kh+3/2)/(kh-3/2)-(1/18)*delta[b]*3^(1/2)*(3*((M-u+3^(1/2)*(mu*sigma[b])^(1/2))^2)^(1/2)*mu-3*((M-u-3^(1/2)*(mu*sigma[b])^(1/2))^2)^(1/2)*mu)/(mu*(mu*sigma[b])^(1/2))-(1/18)*3^(1/2)*(-3*((M+3^(1/2)*sigma[i]^(1/2))^2)^(1/2)+3*((M-3^(1/2)*sigma[i]^(1/2))^2)^(1/2))/sigma[i]^(1/2), -(1-delta[h]-delta[b])*(-kc+3/2)^2/(kc-3/2)^2+(1-delta[h]-delta[b])*(-kc+3/2)/(kc-3/2)^2-delta[h]*(-kh+3/2)^2/(sigma[h]*(kh-3/2)^2)+delta[h]*(-kh+3/2)/((kh-3/2)^2*sigma[h])-(1/18)*delta[b]*3^(1/2)*(3*mu^2/((M-u+3^(1/2)*(mu*sigma[b])^(1/2))^2)^(1/2)-3*mu^2/((M-u-3^(1/2)*(mu*sigma[b])^(1/2))^2)^(1/2))/(mu*(mu*sigma[b])^(1/2))-(1/18)*3^(1/2)*(3/((M+3^(1/2)*sigma[i]^(1/2))^2)^(1/2)-3/((M-3^(1/2)*sigma[i]^(1/2))^2)^(1/2))/sigma[i]^(1/2)]

e2 := simplify(e1) assuming positive;

[-(1/6)*(-delta[b]*sigma[i]^(1/2)*(-3*mu^(1/2)*sigma[b]^(1/2)+3^(1/2)*(M-u))*signum(M-u-3^(1/2)*mu^(1/2)*sigma[b]^(1/2))+delta[b]*(3*mu^(1/2)*sigma[b]^(1/2)+3^(1/2)*(M-u))*sigma[i]^(1/2)*signum(M-u+3^(1/2)*mu^(1/2)*sigma[b]^(1/2))+((M*3^(1/2)-3*sigma[i]^(1/2))*signum(M-3^(1/2)*sigma[i]^(1/2))+(-6*delta[b]+3)*sigma[i]^(1/2)-M*3^(1/2))*mu^(1/2)*sigma[b]^(1/2))/(sigma[i]^(1/2)*sigma[b]^(1/2)*mu^(1/2)), -1+delta[h]+delta[b]+(-2+2*delta[h]+2*delta[b])/(2*kc-3)-delta[h]/sigma[h]-2*delta[h]/((2*kh-3)*sigma[h])-(1/6)*delta[b]*3^(1/2)*mu^(1/2)*(1/abs(M-u+3^(1/2)*mu^(1/2)*sigma[b]^(1/2))-1/abs(-M+u+3^(1/2)*mu^(1/2)*sigma[b]^(1/2)))/sigma[b]^(1/2)-(1/6)*3^(1/2)*(1/(M+3^(1/2)*sigma[i]^(1/2))-1/abs(-M+3^(1/2)*sigma[i]^(1/2)))/sigma[i]^(1/2)]

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;

[0, (((kc-1/2)*(-1+delta[h]+delta[b])*M^4-2*u*(kc-1/2)*(-1+delta[h]+delta[b])*M^3+(-3/2+(1+(-1+delta[h]+delta[b])*u^2+3*(1-delta[h]-delta[b])*sigma[i]-3*sigma[b]*mu*delta[h]+((-3*sigma[b]+1)*delta[b]+3*sigma[b])*mu)*kc+(1/2)*(1-delta[h]-delta[b])*u^2+(3/2)*(-1+delta[h]+delta[b])*sigma[i]+(3/2)*sigma[b]*mu*delta[h]+(3/2)*((sigma[b]-1)*delta[b]-sigma[b])*mu)*M^2+6*(1/2+(-1/3+(-1+delta[h]+delta[b])*sigma[i])*kc+(1/2)*(1-delta[h]-delta[b])*sigma[i])*u*M+((1+3*(1-delta[h]-delta[b])*sigma[i])*u^2+9*((sigma[b]*delta[h]+(sigma[b]-1/3)*delta[b]-sigma[b])*sigma[i]-(1/3)*sigma[b])*mu)*kc+(3/2)*(-1+(-1+delta[h]+delta[b])*sigma[i])*u^2-(9/2)*((sigma[b]*delta[h]+(sigma[b]-1)*delta[b]-sigma[b])*sigma[i]-sigma[b])*mu)*(kh-3/2)*sigma[h]-(kh-1/2)*(M^2-2*M*u-3*mu*sigma[b]+u^2)*(kc-3/2)*(M^2-3*sigma[i])*delta[h])/((M-u-3^(1/2)*mu^(1/2)*sigma[b]^(1/2))*(kc-3/2)*(M^2-3*sigma[i])*sigma[h]*(kh-3/2)*(M-u+3^(1/2)*mu^(1/2)*sigma[b]^(1/2)))]

 

 

Download New-eval-2.mw

 

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

restart;

with(OrthogonalExpansions):

f := x -> piecewise(x<0, 0, x);

f := proc (x) options operator, arrow; piecewise(x < 0, 0, x) end proc

F := FourierSeries(f(x), x=-Pi..Pi, n);

(1/4)*Pi+Sum(((-1)^i-1)*cos(i*x)/(i^2*Pi)-(-1)^i*sin(i*x)/i, i = 1 .. n)

plots:-animate(plot, [F, x=-Pi..3*Pi, color=blue], n=1..10, frames=10);

 

Download mw.mw

 

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 a centered on the x axis; parameters s and p

C1 := <p, a*cos(s), a*sin(s) >;

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

Cylinder of radius ccentered on the zaxis; parameters t and q

C2 := <c*cos(t), c*sin(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

 

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:

alpha: half-angle of the cone's vertex

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

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

beta = .9*alpha, alpha = (1/6)*Pi, d = 1/2

Parametric equation of the cone (s and t 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 N, A, B form an orthonormal triad, with N normal to the
plane and A, B 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. (u and v are the parameters)

P := A*u + B*v + 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

u, v, s, t.  We solve the system for u, v, s in terms of t: and then plug the result into the

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

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

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

 

 

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)

1 2 3 4 5 6 7 Last Page 3 of 54