ijuptilk

How to solve system of equations explici...

Asked: ijuptilk 100 Product: Maple

July 13 2022

1 1

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.

LSA.mw

>

##

>

###

>

###

>

restart:

>	f__1 := Delta -(psi + mu)*S(t);

(1)

>	f__2 := psiS(t) -(delta + mu)E(t);

(2)

>	f__3 := DeltaE(t) -(gamma+gamma__1 + mu)X(t);

(3)

>	f__4 := gamma__1X(t)-(eta + xi + mu)H(t);

(4)

>	f__5 := xiH(t) - muR(t);

(5)

>	f__6 := gammaX(t)-etaH(t) - d*D(t);

(6)

>	f__7 := bD(t) - bB(t);

(7)

>	f__8 := phi__p + sigmaX(t)+eta__1H(t) +d__1D(t)+ b__1B(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__1f__1(1-S__1/S(t));

c__1*(Delta-(psi+mu)*S(t))*(1-(Delta/mu-a__1*E__1/mu)/S(t))

(9)

>	M__2 := c__2f__2(1-E__1/E(t));

(10)

>	M__3 := c__3f__3(1-X__1/X(t));

(11)

>	M__4 :=c__4f__4(1-H__1/H(t));

c__4*(gamma__1*X(t)-(eta+xi+mu)*H(t))*(1-delta*gamma__1*E__1*a__3/(a__2*H(t)))

(12)

>	M__5 := c__5f__5(1-R__1/R(t));

c__5*(xi*H(t)-mu*R(t))*(1-delta*gamma__1*xi*E__1*a__3*mu/(a__2*R(t)))

(13)

>	M__6 := c__6f__6(1-D__1/D(t));

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

(14)

>	M__7 := c__7f__7(1-B__1/B(t));

(15)

>	M__8 := c__8f__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 := gammaX__1 + etaH__1;

(22)

>	u__7 := b*D__1;

(23)

>	u__8 := phi__p + sigmaX__1+eta__1H__1 +d__1D__1+ b__1B__1;

(24)

>	M__1 := c__1(Delta-(psi+mu)S(t))(1-(Delta/mu-a__1E__1/mu)/S(t));

c__1*(Delta-(psi+mu)*S(t))*(1-(Delta/mu-a__1*E__1/mu)/S(t))

(25)

>	M__2 := c__2(psiS(t)-(delta+mu)E(t))(1-E__1/E(t));

(26)

>	M__3 := c__3(DeltaE(t)-(gamma+gamma__1+mu)X(t))(1-deltaE__1/(a__2X(t)));

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__1X(t)-(eta+xi+mu)H(t))(1-deltagamma__1E__1a__3/(a__2H(t)));

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(xiH(t)-muR(t))(1-deltagamma__1xiE__1a__3mu/(a__2R(t)));

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(gammaX(t)-etaH(t)-dD(t))(1-(gammadelta/a__2+deltaetagamma__1a__3/a__2)E__1/(d*D(t)));

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(bD(t)-bB(t))(1-(gammadelta/a__2+deltaetagamma__1a__3/a__2)E__1/(bB(t)));

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+sigmaX(t)+eta__1H(t)+d__1D(t)+b__1B(t)-alphaP(t))(1-((deltasigma/a__2+deltaetagamma__1a__3/a__2+(d__1/d-b__1/b)(gammadelta/a__2+deltaetagamma__1a__3/a__2))*E__1/alpha+phi__p/alpha)/P(t));

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

-c__1*(psi+mu)*(S(t)-S__1)*(a__1*E__1+S(t)*mu-S__1*mu-psi*S__1)/(S(t)*mu)

(34)

>	L__2 := expand(subs((delta+mu)*E(t)=u__2, M__2));

c__2*psi*S(t)-c__2*psi*S(t)*E__1/E(t)-c__2*psi*S__1+c__2*psi*S__1*E__1/E(t)

(35)

>	L__3 := expand(subs((gamma+gamma__1+mu)*X(t)=u__3, M__3));

c__3*Delta*E(t)-c__3*Delta*E(t)*delta*E__1/(a__2*X(t))-c__3*delta*E__1+c__3*delta^2*E__1^2/(a__2*X(t))

(36)

>	L__4 := expand(subs((eta+xi+mu)*H(t)=u__4, M__4));

c__4*gamma__1*X(t)-c__4*gamma__1^2*X(t)*delta*E__1*a__3/(a__2*H(t))-c__4*gamma__1*X__1+c__4*gamma__1^2*X__1*delta*E__1*a__3/(a__2*H(t))

(37)

>	L__5 := expand(subs(mu*R(t)=u__5, M__5));

c__5*xi*H(t)-c__5*xi^2*H(t)*delta*gamma__1*E__1*a__3*mu/(a__2*R(t))-c__5*xi*H__1+c__5*xi^2*H__1*delta*gamma__1*E__1*a__3*mu/(a__2*R(t))

(38)

>	L__6 := expand(subs(d*D(t)=u__6, M__6));

(39)

>	L__7 := expand(subs(b*B(t)=u__7, M__7));

c__7*b*D(t)-c__7*D(t)*E__1*gamma*delta/(B(t)*a__2)-c__7*D(t)*E__1*delta*eta*gamma__1*a__3/(B(t)*a__2)-c__7*b*D__1+c__7*D__1*E__1*gamma*delta/(B(t)*a__2)+c__7*D__1*E__1*delta*eta*gamma__1*a__3/(B(t)*a__2)

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

c__8*b__1^2*E__1*delta*eta*gamma__1*a__3/(P(t)*alpha*b*a__2)+c__8*b__1^2*E__1*gamma*delta/(P(t)*alpha*b*a__2)-c__8*b__1*phi__p/(P(t)*alpha)-c__8*b__1*E__1*delta*eta*gamma__1*a__3/(P(t)*alpha*a__2)-c__8*b__1*E__1*d__1*gamma*delta/(P(t)*alpha*d*a__2)+c__8*b__1-c__8*b__1*E__1*delta*sigma/(P(t)*alpha*a__2)-c__8*b__1*E__1*d__1*delta*eta*gamma__1*a__3/(P(t)*alpha*d*a__2)

(46)

>	k__5 := coeff(L, E(t));

-c__3*Delta*delta*E__1/(a__2*X(t))+c__3*Delta

(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;

-c__8*d__1^2*E__1*delta*eta*gamma__1*a__3/(P(t)*alpha*d*a__2)-c__8*d__1^2*E__1*gamma*delta/(P(t)*alpha*d*a__2)-c__8*d__1*phi__p/(P(t)*alpha)-c__8*d__1*E__1*delta*eta*gamma__1*a__3/(P(t)*alpha*a__2)+c__8*d__1*E__1*b__1*gamma*delta/(P(t)*alpha*b*a__2)+c__8*d__1-c__8*d__1*E__1*delta*sigma/(P(t)*alpha*a__2)+c__8*d__1*E__1*b__1*delta*eta*gamma__1*a__3/(P(t)*alpha*b*a__2)

(51)

>

k__4;

c__8*b__1^2*E__1*delta*eta*gamma__1*a__3/(P(t)*alpha*b*a__2)+c__8*b__1^2*E__1*gamma*delta/(P(t)*alpha*b*a__2)-c__8*b__1*phi__p/(P(t)*alpha)-c__8*b__1*E__1*delta*eta*gamma__1*a__3/(P(t)*alpha*a__2)-c__8*b__1*E__1*d__1*gamma*delta/(P(t)*alpha*d*a__2)+c__8*b__1-c__8*b__1*E__1*delta*sigma/(P(t)*alpha*a__2)-c__8*b__1*E__1*d__1*delta*eta*gamma__1*a__3/(P(t)*alpha*d*a__2)

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

[[c__4 = -H(t)*gamma*(E__1*a__3*delta*eta*gamma__1+E__1*gamma*delta-D(t)*a__2*d)*c__6/(D(t)*d*gamma__1*(E__1*a__3*delta*gamma__1-H(t)*a__2)), c__5 = R(t)*eta*(E__1*a__3*delta*eta*gamma__1+E__1*gamma*delta-D(t)*a__2*d)*c__6/(D(t)*d*xi*(E__1*a__3*delta*mu*xi*gamma__1-R(t)*a__2)), c__6 = c__6, c__8 = 0]]

(53)

>

#

>

#

>

Download LSA.mw

Pairing 2 lists into another list and us...

Asked: ijuptilk 100 Product: Maple

June 21 2022

3 3

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

Asked: ijuptilk 100 Product: Maple 18

May 29 2022

2 21

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

Asked: ijuptilk 100 Product: Maple 18

April 05 2022

1 32

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] := 610^(-12): d:= 0.2e-3: xi:= -0.01: theta0:= 0.1e-3: eta[1]:= 0.240e-1: alpha:= 1-alpha3^2/(gamma1eta[1]): c:= alpha3xialpha/(eta[1](4k[1]q^2/d^2-alpha3xi/eta[1])): theta_init:= theta0sin(Piz/d): n:= 10:

Assign g for q and plot g

>	g := q-(1-alpha)tan(q)-ctan(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
"q= ((2 k+1)Pi)/(2). "
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(2j+1))*Pi), j = 0 .. n)]); ModifiedOddAsym := abs(`~`[`-`](OddAsymptotes, qstar)); qstarTemporary := min(ModifiedOddAsym); indexOfqstar2 := SearchAll(qstarTemporary, ModifiedOddAsym); qstar2 := OddAsymptotes(indexOfqstar2);

$OddAsymptotes := Vector(4, {(1) = ` 1 .. 11 `*Vector[row], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})$

$ModifiedOddAsym := Vector(4, {(1) = ` 1 .. 11 `*Vector[row], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})$

(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:

$OddAsymptotes := Vector(4, {(1) = ` 1 .. 11 `*Vector[row], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})$

$AllAsymptotes := Vector(4, {(1) = ` 1 .. 12 `*Vector[row], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})$

(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] := gamma1alpha/(4k[1]qq[i]^2/d^2-alpha3xi/eta[1]); end do;

(4.7.1)

Minimum of the Re()

>	for i to m do p[i] := min(Re(gamma1alpha/(4k[1]qq[i]^2/d^2-alpha3xi/eta[1]))); end do;

(4.7.1.1)

>	## I expected 0.20459 but return all the entries in the list.

Download Finding_min_zero.mw

Francis QR algorithm...

Asked: ijuptilk 100 Product: Maple 18

March 19 2022

1 17

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.

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

100 Reputation

3 Badges

MaplePrimes Activity

These are questions asked by ijuptilk

How to solve system of equations explici...

Pairing 2 lists into another list and us...

How to use sequence in a procedure...

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

Francis QR algorithm...

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