## 100 Reputation

2 years, 125 days

## How to solve system of equations explici...

Maple

Please I need your assistance. I want to solve for c__4, c__5, c__6, and c__8  from 4 systems of the equation: See my code below:

Since there 4 equations and 4 unknowns, is it possible to get the result explicitly without setting c__6=c__8 as maple did? The solution is at the end of the maple file.

 > ##
 > ###
 > ###
 > restart:
 > f__1 := Delta -(psi + mu)*S(t);
 (1)
 > f__2 := psi*S(t) -(delta + mu)*E(t);
 (2)
 > f__3 := Delta*E(t) -(gamma+gamma__1 + mu)*X(t);
 (3)
 > f__4 := gamma__1*X(t)-(eta + xi + mu)*H(t);
 (4)
 > f__5 := xi*H(t) - mu*R(t);
 (5)
 > f__6 := gamma*X(t)-eta*H(t) - d*D(t);
 (6)
 > f__7 := b*D(t) - b*B(t);
 (7)
 > f__8 := phi__p + sigma*X(t)+eta__1*H(t) +d__1*D(t)+ b__1*B(t) - alpha*P(t);
 (8)
 > S__1 := Delta/mu - a__1*E__1/mu: X__1:= delta*E__1/a__2: H__1 := delta*gamma__1*E__1/a__2*a__3: R__1 := delta*gamma__1*xi*E__1/a__2*a__3*mu: B__1 := 1/b*(gamma*delta/a__2 + delta*eta*gamma__1/a__2*a__3)*E__1: D__1 := 1/d*(gamma*delta/a__2 + delta*eta*gamma__1/a__2*a__3)*E__1: P__1 := (1/alpha)*(delta*sigma/a__2 + delta*eta*gamma__1/a__2*a__3 + (d__1/d - b__1/b)*(delta*gamma/a__2 + delta*eta*gamma__1/a__2*a__3))*E__1 + phi__p/alpha:
 > M__1 := c__1*f__1*(1-S__1/S(t));
 (9)
 > M__2 := c__2*f__2*(1-E__1/E(t));
 (10)
 > M__3 := c__3*f__3*(1-X__1/X(t));
 (11)
 > M__4 :=c__4*f__4*(1-H__1/H(t));
 (12)
 > M__5 := c__5*f__5*(1-R__1/R(t));
 (13)
 > M__6 := c__6*f__6*(1-D__1/D(t));
 (14)
 > M__7 := c__7*f__7*(1-B__1/B(t));
 (15)
 > M__8 := c__8*f__8*(1-P__1/P(t));
 (16)
 > restart:
 > u__1 := (psi+mu)*S__1;
 (17)
 > u__2 := psi*S__1;
 (18)
 > u__3 := delta*E__1;
 (19)
 > u__4 := gamma__1*X__1;
 (20)
 > u__5 := xi*H__1;
 (21)
 > u__6 := gamma*X__1 + eta*H__1;
 (22)
 > u__7 := b*D__1;
 (23)
 > u__8 := phi__p + sigma*X__1+eta__1*H__1 +d__1*D__1+ b__1*B__1;
 (24)
 > M__1 := c__1*(Delta-(psi+mu)*S(t))*(1-(Delta/mu-a__1*E__1/mu)/S(t));
 (25)
 > M__2 := c__2*(psi*S(t)-(delta+mu)*E(t))*(1-E__1/E(t));
 (26)
 > M__3 := c__3*(Delta*E(t)-(gamma+gamma__1+mu)*X(t))*(1-delta*E__1/(a__2*X(t)));
 (27)
 > M__4 := c__4*(gamma__1*X(t)-(eta+xi+mu)*H(t))*(1-delta*gamma__1*E__1*a__3/(a__2*H(t)));
 (28)
 > M__5 := c__5*(xi*H(t)-mu*R(t))*(1-delta*gamma__1*xi*E__1*a__3*mu/(a__2*R(t)));
 (29)
 > M__6 := c__6*(gamma*X(t)-eta*H(t)-d*D(t))*(1-(gamma*delta/a__2+delta*eta*gamma__1*a__3/a__2)*E__1/(d*D(t)));
 (30)
 > M__7 := c__7*(b*D(t)-b*B(t))*(1-(gamma*delta/a__2+delta*eta*gamma__1*a__3/a__2)*E__1/(b*B(t)));
 (31)
 > M__8 := c__8*(phi__p+sigma*X(t)+eta__1*H(t)+d__1*D(t)+b__1*B(t)-alpha*P(t))*(1-((delta*sigma/a__2+delta*eta*gamma__1*a__3/a__2+(d__1/d-b__1/b)*(gamma*delta/a__2+delta*eta*gamma__1*a__3/a__2))*E__1/alpha+phi__p/alpha)/P(t));
 (32)
 > J := M__1 + M__2 + M__3 + M__4 + M__5 + M__6 + M__7 + M__8;
 (33)
 > ##
 > L__1 := factor(expand(subs(Delta=u__1, M__1)));
 (34)
 > L__2 := expand(subs((delta+mu)*E(t)=u__2, M__2));
 (35)
 > L__3 := expand(subs((gamma+gamma__1+mu)*X(t)=u__3, M__3));
 (36)
 > L__4 := expand(subs((eta+xi+mu)*H(t)=u__4, M__4));
 (37)
 > L__5 := expand(subs(mu*R(t)=u__5, M__5));
 (38)
 > L__6 := expand(subs(d*D(t)=u__6, M__6));
 (39)
 > L__7 := expand(subs(b*B(t)=u__7, M__7));
 (40)
 > L__8 := expand(subs(alpha*P(t)=u__8, M__8));
 (41)
 > L := L__1 + L__2 + L__3 + L__4 + L__5 + L__6 + L__7 + L__8;
 (42)
 > ## Collecting the coefficients of X, H, D, B and E
 > k__1 := coeff(L,X(t));
 (43)
 > k__2 := coeff(L,H(t));
 (44)
 > k__3 := coeff(L, D(t));
 (45)
 > k__4 := coeff(L, B(t));
 (46)
 > k__5 := coeff(L, E(t));
 (47)
 > ##
 > ## L terms that not coeffiecient X, H, D, B and E
 > W__12 := coeff(L, X(t), 0):
 > W__1 := coeff(W__12, H(t), 0):
 > W__11 := coeff(W__1, D(t), 0):
 > W__12 := coeff(W__11, B(t), 0):
 > k__6 := coeff(W__12, E(t), 0);
 (48)
 > ## from k__5, c__3 = 0 and choose c__1 = c__2 and c__3 =c__7
 > c__3 := 0:  c__7:=0:
 > k__1;
 (49)
 > k__2;
 (50)
 > k__3;
 (51)
 > k__4;
 (52)
 > ## Also choose
 > sys := solve({k__1=0, k__2 =0, k__3=0, k__4=0}, [c__4,c__5, c__6, c__8]);
 (53)
 > #
 > #
 >

## Pairing 2 lists into another list and us...

Maple

Hi please, how can I pair two lists and form another list?

P := [.6286420119, -.6286420119, 0., 0., 0., 0., 0., 0., 0., 0., 0.]

Q = [2.106333379, 2.106333379, 4.654463885, 7.843624703, 10.99193295, 14.13546782, 17.27782732, 20.41978346, 23.56157073, 26.70327712, 29.84494078]

I want to pair them into something like this:

W := [ [.6286420119,2.106333379], [-.6286420119,  2.106333379] ... ]

Then use pointplot(W) and  display(seq(pointplot(W[j], j = 1 .. 20), insequence = true):

## How to use sequence in a procedure...

Maple 18

Please how can I define two sequences in a procedure with two arguments?

ans:= Array([seq( [j, doCalc(j, u)], j=-2..0, 0.006, u=1..334)]):

The procedure is doCalc(j,u)  and I received this error: invalid input: seq expects between 1 and 3 arguments, but received 4

Find attached my complete code.Seq_Proc.mw

## How to find the minimum value of the rea...

Maple 18

Please, how do I find the minimum of the real part of a complex function? I tried min ( ) function it didn't work.

Find attached the fileFinding_min_zero.mw

Import packages

 > restart: with(ArrayTools): with(Student:-Calculus1): with(LinearAlgebra): with(ListTools):with(RootFinding):with(ListTools):

Parameters

 > gamma1 := .1093: alpha3 := -0.1104e-2: k[1] := 6*10^(-12): d:= 0.2e-3: xi:= -0.01: theta0:= 0.1e-3: eta[1]:= 0.240e-1: alpha:= 1-alpha3^2/(gamma1*eta[1]): c:= alpha3*xi*alpha/(eta[1]*(4*k[1]*q^2/d^2-alpha3*xi/eta[1])): theta_init:= theta0*sin(Pi*z/d): n:= 10:

Assign g for q and plot g

 > g := q-(1-alpha)*tan(q)-c*tan(q): plot(g, q = 0 .. 3*Pi, view = [DEFAULT, -30.. 10]);

Set q as a complex

Assume  and subsitute the result into g and equate the real and complex part to zero, and solve for x and y.

 > f := subs(q = x+I*y, g): b1 := evalc(Re(f)) = 0: b2 := evalc(Im(f)) = 0:

Compute the Special Asymptotes

This asymptote is coming from the  from the definition of

 > qstar := (fsolve(1/c = 0, q = 0 .. infinity)):

Compute Odd asymptote

First, Since , then an asymptote occurs at  In general, we have

Next, we compute the entry of the Oddasymptotes that is close to qstar (special asymptote) as assign it to
ModifiedOaddAsym, and then find the minimum of the ModifiedOaddAsym. Searchall Function returns

the index of an entry in a list.

 > OddAsymptotes := Vector[row]([seq(evalf((1/2*(2*j+1))*Pi), j = 0 .. n)]); ModifiedOddAsym := abs(`~`[`-`](OddAsymptotes, qstar)); qstarTemporary := min(ModifiedOddAsym); indexOfqstar2 := SearchAll(qstarTemporary, ModifiedOddAsym); qstar2 := OddAsymptotes(indexOfqstar2);
 (4.2.1)

Compute x and y

Here, we solve for and  within the min. and max. of qstar2 and qstar, and substitute the results into .

 > AreThereComplexRoots := type(true, 'truefalse'); try    soln1:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = 0 .. infinity});    soln2:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = -infinity .. 0});    qcomplex1 := subs(soln1, x+I*y);    qcomplex2 := subs(soln2, x+I*y); catch:    AreThereComplexRoots := type(FAIL, 'truefalse'); end try;
 (4.3.1)

Compute the rest of the Roots

In this section we compute the roots between each asymptotes.

 > OddAsymptotes := Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]); AllAsymptotes := sort(Vector[row]([OddAsymptotes, qstar])); if AreThereComplexRoots then gg := [qcomplex1, qcomplex2, op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)]; elif not AreThereComplexRoots then gg := [op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)]; end if:
 (4.4.1)

Remove the repeated roots if any

 > qq := MakeUnique(gg):

Redefine n

 > m := numelems(qq):

Compute the time constants

 > for i to m do p[i] := gamma1*alpha/(4*k[1]*qq[i]^2/d^2-alpha3*xi/eta[1]); end do;
 (4.7.1)

Minimum of the Re()

 > for i to m do p[i] := min(Re(gamma1*alpha/(4*k[1]*qq[i]^2/d^2-alpha3*xi/eta[1]))); end do;
 (4.7.1.1)
 > ## I expected 0.20459 but return all the entries in the list.

## Francis QR algorithm...

Maple 18

I have tried a code in python for Francis QR algorithm but didn't desire the result. I don't know if it is possible to code in maple.

Given that A^0 = [[3.0, 1.0, 4.0], [1.0, 2.0, 2.0], [0., 13.0, 2]].

1. Write a little program that computes 1 step of Francis QR and test your program by starting from

A^0 = [[3.0, 1.0, 4.0], [1.0, 2.0, 2.0], [0., 13.0, 2]]  and compute A^1, A^2 ...A^6.

I expected to get:

A^0 = [[3.0, 1.0, 4.0], [1.0, 2.0, 2.0], [0., 13.0, 2]],

A^1 = [[3.5,  -4.264, 0.2688], [-9.206, 1.577, 9.197], [0., -1.41, 1.923]],

... A^6 = [[8.056,  1.596, 8.584], [0.3596, -2.01, -7.86], [0., 2.576*10^(-16), 0.9542]].

But didn't get the same results.

Here is my python code:

# Import packages
import numpy as np
from numpy.linalg import qr # QR from Linear Algebra Library
import scipy.linalg   # SciPy Linear Algebra Library

A = np.array([[3.0, 1.0, 4.0], [1.0, 2.0, 2.0], [0., 13.0, 2]])
p = [1, 2, 3, 4, 5, 6]
for i in range(30):
q, r = qr(A)
a = np.dot(r, q)
if i+1 in p:
print("Iteration {i+1}")
print(A)

I would really appreciate your help.

Thank you.

 1 2 3 4 5 6 Page 5 of 6
﻿