MaplePrimes Questions

The so-called look-and-say sequence is something like:

# Starting digit: 0 
0→one 0→
10→one 1, then one 0→
1110→three 1's, then one 0→
3110→one 3, two 1's, then one 0→
132110→…

Rosetta Code gives an example of how to generate such sequences in Maple; however, the code in that example is lengthy. More to the point, compared to the subsequent Mathematica analogue, it is rather inefficient (though more readable, but that's not the point). 
I slightly modify the original code, and now the modified version is significantly faster than the original one: 
 

restart;

interface(version)NULL

`Standard Worksheet Interface, Maple 2023.1, Windows 10, July 7 2023 Build ID 1723669`

(1)
LookAndSay := proc(m::posint, n::nonnegint := 1, ` $`)::'Vector'(nonnegint);

LookAndSay(5, 0)NULL

Vector[column](%id = 36893490791509382076)

(2)

CodeTools['Usage'](LookAndSay(50), iterations = 4)

memory used=3.99GiB, alloc change=1.92MiB, cpu time=48.73s, real time=25.03s, gc time=39.05s

 

gc()time[real](LookAndSay(50))NULL

25.030

(3)


 

Download LookAndSay.mw

For the sake of convenience, I paste the original Maple code from the page above into here:

generate_seq := proc(s)
	local times, output, i;
	times := 1;
	output := "";
	for i from 2 to StringTools:-Length(s) do
		if (s[i] <> s[i-1]) then
			output := cat(output, times, s[i-1]);
			times := 1; # re-assign
		else 
			times ++;
		end if;
	end do;
	cat(output, times, s[i - 1]);
end proc:

Look_and_Say :=proc(n)
	local value, i;
	value := "1";
	print(value);
	for i from 2 to n do
		value := generate_seq(value);
		print(value);
	end do;
end proc:

#Test:
Look_and_Say(10);

Below is a modified version by me: 

LookAndSay := proc(m::posint, n::nonnegint := 1, $)::'Vector'(nonnegint);
    uses StringTools;
    local delete::'equationlist' := Repeat~(Explode(ExpandCharacterClass(":digit:")), 2), aux::procedure[nonnegint](posint) := proc(_::posint)::nonnegint;
        options cache;
        description "https://rosettacode.org/wiki/Look-and-say_sequence#Maple";
        if _ = 1 then
            n
        else
            local char::'nonemptystring' := String(thisproc(_ - 1)), temp::list(string)(*, str*);
            try
                temp := RegSplit(sprintf("(?=[%s])(?<=(?!%s)[%s])", Unique(char), Join(select[2](foldr, eval(apply), delete, rcurry(OrMap, Unique(char)), curry(curry, Has)), "|"), Unique(char)), char)
            catch "cannot compile regular expression: repetition-operator operand invalid":
                local rule::'equationlist' := (str -> str = Insert(str, 1, " "))~(subs(delete =~ NULL, Generate(2, Support(char))));
                temp := StringSplit(Subs(rule, char), " ") # to be replaced
            end;
            parse(String(seq('length(str), str[1]', str in temp)))
        fi;
    end;
    forget(aux);
    Vector[column](m, aux, datatype = nonnegint)
end:

The form LookAndSay(m, n); produces the first m terms in a "look-and-say" sequence starting with n (which is by default 1). Unfortunately, even if "the modified version is significantly faster than the original one", it is still much slower (25s versus 2.5s) than the uncompiled Mathematica implementation (on the same modern computer):

So, is there a workaround to evaluate LookAndSay(50, 1): in about three seconds (instead of half a minute) in the first call (without looking up a pre-calculated table) (on the same computer) as well? 

This is old question https://www.mapleprimes.com/questions/208909-Code-For-Integer-Points-On-Sphere. Now I see this question at here 
https://mathematica.stackexchange.com/questions/288956/how-can-i-get-all-squares-on-this-sphere-so-that-its-coordinates-are-integer-num
 My idea is select all diameters which diameters are perpendicularly from the sphere (x-2)^2 + (y-4)^2 + (z-6)^2 = 15^2. How can I tell Maple to do that?

how to plot graphs for both methods and comparison of different method values for Diff(f(eta),eta, eta) at eta =0

 

NULL

NULL

restart

F[0] := al

F[1] := a2

F[2] := a3

F[3] := a4

G[0] := a5

G[1] := a6

T[0] := a7

T[1] := a8

Q[0] := a9

Q[1] := a10

n[1] := 1

for k from 0 to n[1] do F[k+4] := solve((1+a)*(k+1)*(k+2)*(k+3)*(k+4)*F[k+4]-a*(k+1)*(k+2)*G[k+2]-R*(sum(F[k-m]*(m+1)*(m+2)*(m+3)*F[m+3], m = 0 .. k))+R*(sum((k-m+1)*F[k-m+1]*(m+1)*(m+2)*F[m+2], m = 0 .. k)), F[k+4]) end do

-(1/12)*(R*a2*a3-3*R*a4*al-a*G[2])/(1+a)

 

-(1/60)*(R^2*a2*a3*al-3*R^2*a4*al^2+2*R*a*a3^2-R*a*al*G[2]+2*R*a3^2-3*a^2*G[3]-3*a*G[3])/(1+a)^2

(1)

n[2] := 3

for k from 0 to n[2] do G[k+2] := solve(b*(k+1)*(k+2)*G[k+2]+a*(k+1)*(k+2)*F[k+2]-2*a*G[k]-c*R*(sum((m+1)*G[m+1]*F[k-m], m = 0 .. k))+c*R*(sum(G[k-m]*(m+1)*F[m+1], m = 0 .. k)), G[k+2]) end do

-(1/2)*(R*a2*a5*c-R*a6*al*c+2*a*a3-2*a*a5)/b

 

-(1/6)*(R^2*a2*a5*al*c^2-R^2*a6*al^2*c^2+2*R*a*a3*al*c-2*R*a*a5*al*c+2*R*a3*a5*b*c+6*a*a4*b-2*a*a6*b)/b^2

 

-(1/24)*(R^3*a*a2*a5*al^2*c^3-R^3*a*a6*al^3*c^3+R^3*a2*a5*al^2*c^3-R^3*a6*al^3*c^3+2*R^2*a^2*a3*al^2*c^2-2*R^2*a^2*a5*al^2*c^2+R^2*a*a2^2*a5*b*c^2-R^2*a*a2*a6*al*b*c^2+2*R^2*a*a3*a5*al*b*c^2+2*R^2*a*a3*al^2*c^2-2*R^2*a*a5*al^2*c^2+R^2*a2^2*a5*b*c^2-R^2*a2*a6*al*b*c^2+2*R^2*a3*a5*al*b*c^2+2*R*a^2*a2*a3*b*c-R*a^2*a2*a5*b*c+6*R*a^2*a4*al*b*c-3*R*a^2*a6*al*b*c+2*R*a*a3*a6*b^2*c+6*R*a*a4*a5*b^2*c-2*R*a*a2*a3*b^2+2*R*a*a2*a3*b*c+6*R*a*a4*al*b^2+6*R*a*a4*al*b*c-4*R*a*a6*al*b*c+2*R*a3*a6*b^2*c+6*R*a4*a5*b^2*c+2*a^3*a3*b-2*a^3*a5*b+4*a^2*a3*b-4*a^2*a5*b)/(b^3*(1+a))

 

-(1/120)*(R^4*a^2*a2*a5*al^3*c^4-R^4*a^2*a6*al^4*c^4+2*R^4*a*a2*a5*al^3*c^4-2*R^4*a*a6*al^4*c^4+R^4*a2*a5*al^3*c^4-R^4*a6*al^4*c^4+2*R^3*a^3*a3*al^3*c^3-2*R^3*a^3*a5*al^3*c^3+3*R^3*a^2*a2^2*a5*al*b*c^3-3*R^3*a^2*a2*a6*al^2*b*c^3+2*R^3*a^2*a3*a5*al^2*b*c^3+4*R^3*a^2*a3*al^3*c^3-4*R^3*a^2*a5*al^3*c^3+6*R^3*a*a2^2*a5*al*b*c^3-6*R^3*a*a2*a6*al^2*b*c^3+4*R^3*a*a3*a5*al^2*b*c^3+2*R^3*a*a3*al^3*c^3-2*R^3*a*a5*al^3*c^3+3*R^3*a2^2*a5*al*b*c^3-3*R^3*a2*a6*al^2*b*c^3+2*R^3*a3*a5*al^2*b*c^3+6*R^2*a^3*a2*a3*al*b*c^2-4*R^2*a^3*a2*a5*al*b*c^2+6*R^2*a^3*a4*al^2*b*c^2-4*R^2*a^3*a6*al^2*b*c^2+4*R^2*a^2*a2*a3*a5*b^2*c^2-R^2*a^2*a2*a5^2*b^2*c^2+2*R^2*a^2*a3*a6*al*b^2*c^2+6*R^2*a^2*a4*a5*al*b^2*c^2+R^2*a^2*a5*a6*al*b^2*c^2-2*R^2*a^2*a2*a3*al*b^2*c+12*R^2*a^2*a2*a3*al*b*c^2-R^2*a^2*a2*a5*al*b^2*c-6*R^2*a^2*a2*a5*al*b*c^2+6*R^2*a^2*a4*al^2*b^2*c+12*R^2*a^2*a4*al^2*b*c^2+R^2*a^2*a6*al^2*b^2*c-10*R^2*a^2*a6*al^2*b*c^2-2*R^2*a*a2*a3*a5*b^3*c+8*R^2*a*a2*a3*a5*b^2*c^2-R^2*a*a2*a5^2*b^2*c^2+4*R^2*a*a3*a6*al*b^2*c^2+6*R^2*a*a4*a5*al*b^3*c+12*R^2*a*a4*a5*al*b^2*c^2+R^2*a*a5*a6*al*b^2*c^2-2*R^2*a*a2*a3*al*b^3-2*R^2*a*a2*a3*al*b^2*c+6*R^2*a*a2*a3*al*b*c^2-2*R^2*a*a2*a5*al*b*c^2+6*R^2*a*a4*al^2*b^3+6*R^2*a*a4*al^2*b^2*c+6*R^2*a*a4*al^2*b*c^2-6*R^2*a*a6*al^2*b*c^2-2*R^2*a2*a3*a5*b^3*c+4*R^2*a2*a3*a5*b^2*c^2+2*R^2*a3*a6*al*b^2*c^2+6*R^2*a4*a5*al*b^3*c+6*R^2*a4*a5*al*b^2*c^2+4*R*a^4*a3*al*b*c-4*R*a^4*a5*al*b*c+12*R*a^3*a2*a4*b^2*c-4*R*a^3*a2*a6*b^2*c+2*R*a^3*a5^2*b^2*c+12*R*a^2*a4*a6*b^3*c-2*R*a^3*a3*al*b^2+12*R*a^3*a3*al*b*c+2*R*a^3*a5*al*b^2-12*R*a^3*a5*al*b*c+24*R*a^2*a2*a4*b^2*c-8*R*a^2*a2*a6*b^2*c-4*R*a^2*a3^2*b^3+4*R*a^2*a3*a5*b^2*c+2*R*a^2*a5^2*b^2*c+24*R*a*a4*a6*b^3*c+8*R*a^2*a3*al*b*c-8*R*a^2*a5*al*b*c+12*R*a*a2*a4*b^2*c-4*R*a*a2*a6*b^2*c-4*R*a*a3^2*b^3+4*R*a*a3*a5*b^2*c+12*R*a4*a6*b^3*c+6*a^4*a4*b^2-2*a^4*a6*b^2+18*a^3*a4*b^2-6*a^3*a6*b^2+12*a^2*a4*b^2-4*a^2*a6*b^2)/(b^4*(1+a)^2)

(2)

n[3] := 3

for k from 0 to n[3] do T[k+2] := solve((k+1)*(k+2)*T[k+2]+p3*(k+1)*(k+2)*Q[k+2]+p1*(sum((m+1)*F[m+1]*T[k-m], m = 0 .. k))-p1*(sum(F[k-m]*(m+1)*T[m+1], m = 0 .. k)), T[k+2]) end do

-(1/2)*p1*a2*a7+(1/2)*p1*al*a8-p3*Q[2]

 

-(1/6)*a2*a7*al*p1^2+(1/6)*a8*al^2*p1^2-(1/3)*al*p1*p3*Q[2]-(1/3)*a3*a7*p1-p3*Q[3]

 

-p3*Q[4]-(1/24)*p1^2*a2^2*a7+(1/24)*a2*p1^2*al*a8-(1/12)*p1*a2*p3*Q[2]-(1/12)*p1*a3*a8-(1/4)*p1*a4*a7-(1/24)*a2*a7*al^2*p1^3+(1/24)*a8*al^3*p1^3-(1/12)*al^2*p1^2*p3*Q[2]-(1/12)*al*a3*a7*p1^2-(1/4)*p1*al*p3*Q[3]

 

(1/120)*(-a*a2*a7*al^3*b*p1^4+a*a8*al^4*b*p1^4-2*a*al^3*b*p1^3*p3*Q[2]-a2*a7*al^3*b*p1^4+a8*al^4*b*p1^4-3*a*a2^2*a7*al*b*p1^3+3*a*a2*a8*al^2*b*p1^3-2*a*a3*a7*al^2*b*p1^3-2*al^3*b*p1^3*p3*Q[2]-6*a*a2*al*b*p1^2*p3*Q[2]-6*a*al^2*b*p1^2*p3*Q[3]-3*a2^2*a7*al*b*p1^3+3*a2*a8*al^2*b*p1^3-2*a3*a7*al^2*b*p1^3+R*a*a2*a5*a7*c*p1-R*a*a6*a7*al*c*p1-4*a*a2*a3*a7*b*p1^2-2*a*a3*a8*al*b*p1^2-6*a*a4*a7*al*b*p1^2-6*a2*al*b*p1^2*p3*Q[2]-6*al^2*b*p1^2*p3*Q[3]+2*R*a2*a3*a7*b*p1-6*R*a4*a7*al*b*p1-12*a*a2*b*p1*p3*Q[3]-24*a*al*b*p1*p3*Q[4]-4*a2*a3*a7*b*p1^2-2*a3*a8*al*b*p1^2-6*a4*a7*al*b*p1^2+2*a^2*a3*a7*p1-2*a^2*a5*a7*p1-12*a*a4*a8*b*p1-12*a2*b*p1*p3*Q[3]-24*al*b*p1*p3*Q[4]-120*a*b*p3*Q[5]-12*a4*a8*b*p1-120*b*p3*Q[5])/(b*(1+a))

(3)

n[4] := 3

for k from 0 to n[4] do Q[k+2] := solve((k+1)*(k+2)*Q[k+2]+p4*(k+1)*(k+2)*Q[k+2]+p2*(sum((m+1)*F[m+1]*Q[k-m], m = 0 .. k))-p2*(sum(F[k-m]*(m+1)*Q[m+1], m = 0 .. k)), Q[k+2]) end do

(1/2)*p2*(a10*al-a2*a9)/(p4+1)

 

(1/6)*p2*(a10*al^2*p2-a2*a9*al*p2-2*a3*a9*p4-2*a3*a9)/(p4+1)^2

 

(1/24)*p2*(a10*al^3*p2^2-a2*a9*al^2*p2^2+a10*a2*al*p2*p4-a2^2*a9*p2*p4-2*a3*a9*al*p2*p4+a10*a2*al*p2-2*a10*a3*p4^2-a2^2*a9*p2-2*a3*a9*al*p2-6*a4*a9*p4^2-4*a10*a3*p4-12*a4*a9*p4-2*a10*a3-6*a4*a9)/(p4+1)^3

 

(1/120)*p2*(a*a10*al^4*b*p2^3-a*a2*a9*al^3*b*p2^3+R*a*a2*a5*a9*c*p4^3-R*a*a6*a9*al*c*p4^3+3*a*a10*a2*al^2*b*p2^2*p4-3*a*a2^2*a9*al*b*p2^2*p4-2*a*a3*a9*al^2*b*p2^2*p4+a10*al^4*b*p2^3-a2*a9*al^3*b*p2^3+3*R*a*a2*a5*a9*c*p4^2-3*R*a*a6*a9*al*c*p4^2+2*R*a2*a3*a9*b*p4^3-6*R*a4*a9*al*b*p4^3+3*a*a10*a2*al^2*b*p2^2-2*a*a10*a3*al*b*p2*p4^2-3*a*a2^2*a9*al*b*p2^2-4*a*a2*a3*a9*b*p2*p4^2-2*a*a3*a9*al^2*b*p2^2-6*a*a4*a9*al*b*p2*p4^2+3*a10*a2*al^2*b*p2^2*p4-3*a2^2*a9*al*b*p2^2*p4-2*a3*a9*al^2*b*p2^2*p4+3*R*a*a2*a5*a9*c*p4-3*R*a*a6*a9*al*c*p4+6*R*a2*a3*a9*b*p4^2-18*R*a4*a9*al*b*p4^2+2*a^2*a3*a9*p4^3-2*a^2*a5*a9*p4^3-4*a*a10*a3*al*b*p2*p4-12*a*a10*a4*b*p4^3-8*a*a2*a3*a9*b*p2*p4-12*a*a4*a9*al*b*p2*p4+3*a10*a2*al^2*b*p2^2-2*a10*a3*al*b*p2*p4^2-3*a2^2*a9*al*b*p2^2-4*a2*a3*a9*b*p2*p4^2-2*a3*a9*al^2*b*p2^2-6*a4*a9*al*b*p2*p4^2+R*a*a2*a5*a9*c-R*a*a6*a9*al*c+6*R*a2*a3*a9*b*p4-18*R*a4*a9*al*b*p4+6*a^2*a3*a9*p4^2-6*a^2*a5*a9*p4^2-2*a*a10*a3*al*b*p2-36*a*a10*a4*b*p4^2-4*a*a2*a3*a9*b*p2-6*a*a4*a9*al*b*p2-4*a10*a3*al*b*p2*p4-12*a10*a4*b*p4^3-8*a2*a3*a9*b*p2*p4-12*a4*a9*al*b*p2*p4+2*R*a2*a3*a9*b-6*R*a4*a9*al*b+6*a^2*a3*a9*p4-6*a^2*a5*a9*p4-36*a*a10*a4*b*p4-2*a10*a3*al*b*p2-36*a10*a4*b*p4^2-4*a2*a3*a9*b*p2-6*a4*a9*al*b*p2+2*a^2*a3*a9-2*a^2*a5*a9-12*a*a10*a4*b-36*a10*a4*b*p4-12*a10*a4*b)/((p4+1)^4*b*(1+a))

(4)

U[1] := sum(F[r]*t^r, r = 0 .. n[1]+4)

p[1] := subs(R = 1, a = 1, b = 1, c = 1, p1 = 1, p2 = .8, p3 = .1, p4 = .1, U[1])

U[2] := sum(G[r]*t^r, r = 0 .. n[2]+2)

p[2] := subs(R = 1, a = 1, b = 1, c = 1, p1 = 1, p2 = .8, p3 = .1, p4 = .1, U[2])

U[3] := sum(T[r]*t^r, r = 0 .. n[2]+2)

p[3] := subs(R = 1, a = 1, b = 1, c = 1, p1 = 1, p2 = .8, p3 = .1, p4 = .1, U[3])

U[4] := sum(Q[r]*t^r, r = 0 .. n[2]+2)

p[4] := subs(R = 1, a = 1, b = 1, c = 1, p1 = 1, p2 = .8, p3 = .1, p4 = .1, U[4])

e1 := subs(t = -1, p[1]) = 0

e2 := subs(t = -1, diff(p[1], t)) = 0

e3 := subs(t = 1, diff(p[1], t)) = -1

e4 := subs(t = 1, p[1]) = 0

e5 := subs(t = -1, p[2]) = 0

e6 := subs(t = 1, p[2]) = 1

e7 := subs(t = -1, p[3]) = 1

e8 := subs(t = 1, p[3]) = 0

e9 := subs(t = -1, p[4]) = 1

e10 := subs(t = 1, p[4]) = 0

j := {e1, e10, e2, e3, e4, e5, e6, e7, e8, e9}

j := solve(j)

sj := evalf(j)

{a10 = -3.476623407, a2 = -5.754056209, a3 = .1776219452, a4 = 11.75811242, a5 = 1.324264301, a6 = -684.5523526, a7 = -.2700369914, a8 = 1.152227714, a9 = 2.191204245, al = 0.3618902741e-1}, {a10 = -.5218741555, a2 = .2575353882, a3 = -.2672619833, a4 = -.2650707765, a5 = 0.7065354871e-1, a6 = .1172581545, a7 = .6100817436, a8 = -.5277387253, a9 = .5842364534, al = .2586309916}, {a10 = -4.849411034, a2 = 11.61910224, a3 = -20.01600142, a4 = -22.98820448, a5 = -303.7401922, a6 = -153.4446663, a7 = -7.896832028, a8 = -4.917031955, a9 = -9.645684059, al = 10.13300071}, {a10 = -12.41434918+6.055636678*I, a2 = -6.912869603-3.362489448*I, a3 = -9.364948739-.7062944755*I, a4 = 14.07573921+6.724978896*I, a5 = -106.6284397-3.087774395*I, a6 = 184.4202683+38.56644530*I, a7 = 2.689687372-4.048821750*I, a8 = -4.715343127+5.167588829*I, a9 = 8.474095612-5.785653488*I, al = 4.807474369+.3531472377*I}, {a10 = -8.462156658-37.78952093*I, a2 = -22.10322629+.7748996783*I, a3 = -2.926063539-87.71943544*I, a4 = 44.45645258-1.549799357*I, a5 = 126.1645842+1357.517358*I, a6 = -880.5344239+73.01362458*I, a7 = -96.56841781+19.40514883*I, a8 = -11.30265439-58.49348719*I, a9 = -59.25678527+13.86225901*I, al = 1.588031769+43.85971772*I}, {a10 = 21.28781597+0.9115942334e-2*I, a2 = -2.190767380-.1297694199*I, a3 = 0.4834062985e-1-8.617807139*I, a4 = 4.631534761+.2595388398*I, a5 = -1.070222696-4.103740084*I, a6 = 28.93315819+1.060309794*I, a7 = -.6440073083+2.959900705*I, a8 = 3.178056838-1.712994921*I, a9 = -1.124006374+8.865509135*I, al = .1008296851+4.308903570*I}, {a10 = -2.226772562-4.893664011*I, a2 = -5.213384606-.4953312060*I, a3 = 1.881656676-24.64377975*I, a4 = 10.67676921+.9906624121*I, a5 = -5.922885277-14.38776520*I, a6 = 9.281006594-6.268746147*I, a7 = -8.563253672+2.519226454*I, a8 = -2.293245547-7.112743663*I, a9 = -4.948019289+2.035858706*I, al = -.8158283379+12.32188987*I}, {a10 = -3.311080211+1.380948844*I, a2 = -6.825505968+3.517539795*I, a3 = 10.11566715-.6387142267*I, a4 = 13.90101194-7.035079589*I, a5 = 106.6696011-4.144959139*I, a6 = 183.4179274-43.03852019*I, a7 = -1.117431335-0.4722817327e-1*I, a8 = -1.705921790+.2164542338*I, a9 = -2.431505210+.6185873236*I, al = -4.932833576+.3193571133*I}, {a10 = 1.720689325, a2 = 11.30494181, a3 = 20.89441402, a4 = -22.35988362, a5 = 304.5741226, a6 = -141.0519632, a7 = -3.607319024, a8 = 2.107261122, a9 = -3.764007990, al = -10.32220701}, {a10 = -3.311080211-1.380948844*I, a2 = -6.825505968-3.517539795*I, a3 = 10.11566715+.6387142267*I, a4 = 13.90101194+7.035079589*I, a5 = 106.6696011+4.144959139*I, a6 = 183.4179274+43.03852019*I, a7 = -1.117431335+0.4722817327e-1*I, a8 = -1.705921790-.2164542338*I, a9 = -2.431505210-.6185873236*I, al = -4.932833576-.3193571133*I}, {a10 = -2.226772562+4.893664011*I, a2 = -5.213384606+.4953312060*I, a3 = 1.881656676+24.64377975*I, a4 = 10.67676921-.9906624121*I, a5 = -5.922885277+14.38776520*I, a6 = 9.281006594+6.268746147*I, a7 = -8.563253672-2.519226454*I, a8 = -2.293245547+7.112743663*I, a9 = -4.948019289-2.035858706*I, al = -.8158283379-12.32188987*I}, {a10 = 21.28781597-0.9115942334e-2*I, a2 = -2.190767380+.1297694199*I, a3 = 0.4834062985e-1+8.617807139*I, a4 = 4.631534761-.2595388398*I, a5 = -1.070222696+4.103740084*I, a6 = 28.93315819-1.060309794*I, a7 = -.6440073083-2.959900705*I, a8 = 3.178056838+1.712994921*I, a9 = -1.124006374-8.865509135*I, al = .1008296851-4.308903570*I}, {a10 = -8.462156658+37.78952093*I, a2 = -22.10322629-.7748996783*I, a3 = -2.926063539+87.71943544*I, a4 = 44.45645258+1.549799357*I, a5 = 126.1645842-1357.517358*I, a6 = -880.5344239-73.01362458*I, a7 = -96.56841781-19.40514883*I, a8 = -11.30265439+58.49348719*I, a9 = -59.25678527-13.86225901*I, al = 1.588031769-43.85971772*I}, {a10 = -12.41434918-6.055636678*I, a2 = -6.912869603+3.362489448*I, a3 = -9.364948739+.7062944755*I, a4 = 14.07573921-6.724978896*I, a5 = -106.6284397+3.087774395*I, a6 = 184.4202683-38.56644530*I, a7 = 2.689687372+4.048821750*I, a8 = -4.715343127-5.167588829*I, a9 = 8.474095612+5.785653488*I, al = 4.807474369-.3531472377*I}

(5)

p[1] := subs(a10 = -.5218741555, a2 = .2575353882, a3 = -.2672619833, a4 = -.2650707765, a5 = 0.7065354871e-1, a6 = .1172581545, a7 = .6100817436, a8 = -.5277387253, a9 = .5842364534, al = .2586309916, p[1])

.2586309916+.2575353882*t-.2672619833*t^2-.2650707765*t^3+0.8630991633e-2*t^4+0.7535388242e-2*t^5

(6)

p[2] := subs(a10 = -.5218741555, a2 = .2575353882, a3 = -.2672619833, a4 = -.2650707765, a5 = 0.7065354871e-1, a6 = .1172581545, a7 = .6100817436, a8 = -.5277387253, a9 = .5842364534, al = .2586309916, p[2])

0.7065354871e-1+.1172581545*t+.3439809338*t^2+.3401058738*t^3+0.8536551748e-1*t^4+0.4263597162e-1*t^5

(7)

p[3] := subs(a10 = -.5218741555, a2 = .2575353882, a3 = -.2672619833, a4 = -.2650707765, a5 = 0.7065354871e-1, a6 = .1172581545, a7 = .6100817436, a8 = -.5277387253, a9 = .5842364534, al = .2586309916, p[3])

.6100817436-.5277387253*t-.1364241818*t^2+0.3945483872e-1*t^3+0.2634243820e-1*t^4-0.1171611337e-1*t^5

(8)

p[4] := subs(a10 = -.5218741555, a2 = .2575353882, a3 = -.2672619833, a4 = -.2650707765, a5 = 0.7065354871e-1, a6 = .1172581545, a7 = .6100817436, a8 = -.5277387253, a9 = .5842364534, al = .2586309916, p[4])

.5842364534-.5218741555*t-.1037943244*t^2+0.3134539737e-1*t^3+0.1955787096e-1*t^4-0.9471241840e-2*t^5

(9)

NULL

value*of*D@@2*F(0)*For*R = 1, 1.5, `and`(2*Using*Both*DTM*scheme, dsolve*method)

 

Download DTM_practice.mw

how I can extract data from plot3d such as x, y, phi2?

thanks

tez-1.mw

Download tez-1.mw
 

restart

``

beta := 2.5; lambda := 0.1e-1; b := Pi; a := Pi; alpha := 1; y[1] := 1.5; y[2] := 1.5; x[1] := -1; x[2] := 1; Q[1] := 40; Q[2] := 35

2.5

 

0.1e-1

 

Pi

 

Pi

 

1

 

1.5

 

1.5

 

-1

 

1

 

40

 

35

(1)

NULL

NULL

v := (2*n-1)*Pi/(2*b)

n-1/2

(2)

Delta := exp(2*v*a)*(alpha*v+beta)*(1+lambda)-(1-lambda)*(alpha*v-beta)

1.01*exp(2*(n-1/2)*Pi)*(n+2.000000000)-.99*n+2.970000000

(3)

g[22] := ((alpha*v+beta)*((1+lambda)*exp(-v*abs(x-xi))+(-1+lambda)*exp(-v*(x+xi)))*exp(2*v*a)+(alpha*v-beta)*((1+lambda)*exp(-v*(x+xi))+(-1+lambda)*exp(-v*abs(x-xi))))/(2*v*Delta)

g[21] := ((alpha*v+beta)*exp(v*(2*a+xi))+(alpha*v-beta)*exp(-v*xi))*exp(-v*x)/(v*Delta)

NULL

u[2] := int(2*g[21]*Q[1]*Dirac(xi-x[1])*sin(n*Pi*y[1]/b)/b, xi = -a .. 0)+int(2*g[22]*Q[2]*Dirac(xi-x[2])*sin(n*Pi*y[2]/b)/b, xi = 0 .. infinity)

NULL

phi[2] := sum(u[2](x)*sin(v*y), n = 1 .. 30)

NULL

``

plot3d(phi[2], x = 0 .. 5, y = 0 .. b)

 

NULL


 

Download tez-1.mw

 

 

Hello 
I have two questions if you please.
Q1)
Is there any kind of a ruler in the maple worksheet that helps me to Adjust the wideness of the line text? 0r even to  spilt the outputs in two columns 
you know as the ruler in the "Word document ".
Q2) Is there any way"method" to change the output type 
since the default adjustment of the outputs is a type of stretchable image? But I need the outputs of the text type. That may  help me to insert them in a given templet

Thank you

Hello everyone, 

I am trying to solve 6 equations with 6 unknowns in maple. I tried to simplify them and they are shown below. Unfortunately, after using solve I get:

Warning, solutions may have been lost
                            sol := ()

Can you help me with it, please? I tried to solve it with Matlab and got a warning that no explicit solution could be found and had no luck with mathematica either.  Thank you very much!!

restart;

omega_1 := 1;
beta_1 = 0.2;
omega_2 := 0.15;
beta_2 := 0.15;
mu := 0.4;
g := 9.81;
K := 5;
vb := 0.5;
tc := arccos(vb/(A*omega));

Ns := mu*(omega_2^2 - K)*B*(2*omega*tc - sin(2*omega*tc) - pi)/pi;

Nc := 4*mu*(omega_2^2 - K)*x2*sin(omega*tc)/pi + mu*(omega_2^2 - K)*C*(2*omega*tc + 2*sin(2*omega*tc) - pi)/pi;

a0 := -mu*(omega_2^2 - K)*x2*(1 - 2*omega*tc/pi) + 2*mu*(omega_2^2 - K)*C*sin(omega*tc)/pi;

eq1 := -A*omega^2 + A*omega_1^2 - B*K - Ns;

eq2 := 2*A*beta_1*omega*omega_1 - C*K - Nc;

eq3 := -2*C*beta_2*omega*omega_2 - B*omega^2 + B*omega_2^2 - A*K;

eq4 := 2*B*beta_2*omega*omega_2 - C*omega^2 + C*omega_2^2;

eq5 := omega_1^2*x1 - K*x2 - a0;

eq6 := omega_2^2*x2 - K*x1 - g;

sol := solve({eq1, eq2, eq3, eq4, eq5, eq6}, {A, B, C, omega, x1, x2});

Hello!

I have a Maple program for plotting from z plane to w plane under transformation w=1/z

I have lines in the z plane and should have had nice circles in the w plane, but the circles are not so good. It is a hole around zero.I can do it manually, but I like to do it with the Maple program.

What must be done in the program to get nice circles???

Thanks for any answers!!

Regards

Kjell

restart;

eq32r:=diff(B(r),`$`(r,3))+diff(B(r),r);

diff(diff(diff(B(r), r), r), r)+diff(B(r), r)

(1)

with(PDEtools, casesplit, declare); with(DEtools, gensys):declare(b(r), prime = r);

[casesplit, declare]

 

b(r)*`will now be displayed as`*b

 

`derivatives with respect to`*r*`of functions of one variable will now be displayed with '`

(2)

eq32r;

diff(diff(diff(B(r), r), r), r)+diff(B(r), r)

(3)

latex(eq32r);

{\frac {{\rm d}^{3}}{{\rm d}{r}^{3}}}B \left( r \right) +{\frac

{\rm d}{{\rm d}r}}B \left( r \right)

 

#instead of {\frac {{\rm d}^{3}}{{\rm d}{r}^{3}}}B I want to have B''' in my latex file.

``

``

Download convert-latex.mw

Hi everybody

I have a differential equation where derivatives (for example d/dr) have been displayed by prime notation.

I want to convert this equation latex format by keeping the prime notation not d/dr.

I would appreciate it if anyone could help me.

I have attached a sample code.

Thank you

Hadi

I'm having problems trying to plot a discrete plot for a Poisson Distribution CDF.  I followed a model found elsewhere that someone had used to do a discrete plot of a sum.  The OP, in that case, was looking to plot a partial sum of:

sum(1/(n^3*sin^2(n)),n=1..400);

The solution that was suggested to the OP was the following:
 

with(DynamicSystems):
f[N1_]=sum[1/(n^3*sin^2[n]), {n,1,N1}];
DiscretePlot[f[x], {x,0,400},PlotRange->All];



Using that solution as a guide, I tried this:
 

with(DynamicSystems):
lambda := 15.4;
f[N1_] = sum[exp(-lambda)*lambda^n/n!, {0, N1, n}];
DiscretePlot[f[x], {0, 50, x}, PlotRange -> All];

In the attached image, you can see the output from my attempt.

I am completely flummoxed. Any suggestions would be greatly appreciated.

Dear Maple Expeets, In the attached Maple file, I have three procedures that are run smoothly. I am trying to compare outputs from each procedure to see which one if larger than the others. For changes on two parameters (alpha and delta), I used implicitplot and inequal commands to identify each region by a specific color. But the outcome is not reasonable. Would you please help with this? Thanks in advance!

NULLNULL

restart

with(plots)

c := 1; cr := 0.3e-1*c; u := 1; sExp := 0.6e-1*c; s := .65*c; v := 3*c

NULL

FirmModelPP := proc (alpha, delta) local p0, xi0, q0, Firmpf0, G0, Recpf0, Unsold0, Environ0, SoldPre0; option remember; xi0 := 1; p0 := min(s+sqrt((v-s)*(c-s)), delta*v+sExp); q0 := u*(v-p0)/(v-s); f(N) := 1/u; F(N) := N/u; G0 := int(F(N), N = 0 .. q0); Firmpf0 := (p0-c)*q0-(p0-s)*G0; Recpf0 := (sExp-cr)*xi0*q0; Environ0 := q0+G0; Unsold0 := G0; SoldPre0 := 0; return p0, q0, Firmpf0, Recpf0, Environ0, Unsold0, SoldPre0 end proc

NULL

FirmModelFC := proc (alpha, beta, delta) local p00, xi00, q00, Firmpf00, G00, Recpf00, Unsold00, Environ00, pr00, SoldPre00; option remember; xi00 := 1; p00 := s+sqrt((v-s)*(c-s)); if p00 < delta*v+sExp then q00 := u*(v-p00)/(v-s); f(N) := 1/u; F(N) := N/u; G00 := int(F(N), N = 0 .. q00); Firmpf00 := (p00-c)*q00-(p00-s)*G00; Recpf00 := `&xi;00*q00*`(sExp-cr); Unsold00 := G00; Environ00 := q00+Unsold00 else q00 := alpha*u*(v-p00)/(v-s); f(N) := 1/u; F(N) := N/u; G00 := int(F(N), N = 0 .. q00/alpha); pr00 := p00-delta*v; Firmpf00 := (p00-c)*q00-alpha*(p00-s)*G00; Recpf00 := (sExp-cr)*xi00*q00+(beta*xi00*q00-(1/2)*beta^2*xi00^2*q00^2/(u*(1-alpha)))*(pr00-sExp); Unsold00 := G00; Environ00 := q00+Unsold00; SoldPre00 := beta*xi00*q00-(1/2)*beta^2*xi00^2*q00^2/(u*(1-alpha)) end if; return p00, q00, Firmpf00, Recpf00, Environ00, Unsold00, SoldPre00 end proc

NULLNULL

NULL

FirmModelHmax := proc (alpha, beta, delta) local q, p, pr, FirmpfSiS, F1, G1, G2, G3, RecpfSiS, sol, UnsoldSiS, EnvironSiS, p0, OldSoldPrim, xi, h; option remember; xi := 1; if alpha <= 1/(1+beta*xi) then p := max(`assuming`([solve(u*(psol-c+(psol-delta*v-sExp)*beta*xi)/(beta^2*xi^2*(psol-delta*v-sExp)/(1-alpha)-(beta^2*xi^2/(1-alpha)-(1+beta*xi)^2)*(psol-s)) = alpha*u*(v-psol)/(v-s), psol, useassumptions)], [0 < psol])); q := alpha*u*(v-p)/(v-s); G2 := (1/2)*beta^2*xi^2*q^2/(u*(1-alpha)^2); G3 := (1/2)*q^2*(1+beta*xi)^2/u; h := (p-delta*v-sExp)/(p-delta*v); FirmpfSiS := (p-c)*q+(p-s)*((1-alpha)*G2-G3)+h*(p-delta*v)*(beta*xi*q-(1-alpha)*G2); RecpfSiS := ((1-h)*(p-delta*v)-sExp)*(beta*xi*q-G2)+sExp*xi*q-cr*xi*q; UnsoldSiS := G3-(1-alpha)*G2; EnvironSiS := q+UnsoldSiS; OldSoldPrim := beta*xi*q-(1-alpha)*G2 else p := max(`assuming`([solve(u*(psol-c+(psol-delta*v-sExp)*beta*xi)/((psol-s)/u+beta^2*xi^2*(psol-delta*v-sExp)/(1-alpha)) = alpha*u*(v-psol)/(v-s), psol, useassumptions)], [0 < psol])); q := alpha*u*(v-p)/(v-s); F1 := beta*xi*q/(u*(1-alpha)); G1 := (1/2)*q^2/(u*alpha^2); G2 := (1/2)*beta^2*xi^2*q^2/(u*(1-alpha)^2); G3 := (1/2)*q^2*(1+beta*xi)^2/u; h := (p-delta*v-sExp)/(p-delta*v); FirmpfSiS := (p-c)*q-alpha*(p-s)*G1+h*(p-delta*v)*(beta*xi*q-(1-alpha)*G2); RecpfSiS := ((1-h)*(p-delta*v)-sExp)*(beta*xi*q-G2)+sExp*xi*q-cr*xi*q; UnsoldSiS := alpha*G1; EnvironSiS := q+UnsoldSiS; OldSoldPrim := beta*xi*q-(1-alpha)*G2 end if; return p, q, FirmpfSiS, RecpfSiS, EnvironSiS, h, UnsoldSiS, OldSoldPrim, xi end proc

NULL

NULLNULL

FirmModelH := proc (alpha, beta, delta, h) local q, p, pr, FirmpfSiS, F1, G1, G2, G3, RecpfSiS, sol, UnsoldSiS, EnvironSiS, p0, OldSoldPrim, xi; option remember; xi := 1; if alpha <= 1/(1+beta*xi) then p := max(`assuming`([solve(u*(psol-c+h*(psol-delta*v)*beta*xi)/(beta^2*xi^2*h*(psol-delta*v)/(1-alpha)-(beta^2*xi^2/(1-alpha)-(1+beta*xi)^2)*(psol-s)) = alpha*u*(v-psol)/(v-s), psol, useassumptions)], [0 < psol])); q := alpha*u*(v-p)/(v-s); G2 := (1/2)*beta^2*xi^2*q^2/(u*(1-alpha)^2); G3 := (1/2)*q^2*(1+beta*xi)^2/u; FirmpfSiS := (p-c)*q+(p-s)*((1-alpha)*G2-G3)+h*(p-delta*v)*(beta*xi*q-(1-alpha)*G2); RecpfSiS := ((1-h)*(p-delta*v)-sExp)*(beta*xi*q-G2)+sExp*xi*q-cr*xi*q; UnsoldSiS := G3-(1-alpha)*G2; EnvironSiS := q+UnsoldSiS; OldSoldPrim := beta*xi*q-(1-alpha)*G2 else p := max(`assuming`([solve(u*(psol-c+h*(psol-delta*v)*beta*xi)/((psol-s)/u+beta^2*xi^2*h*(psol-delta*v)/(1-alpha)) = alpha*u*(v-psol)/(v-s), psol, useassumptions)], [0 < psol])); q := alpha*u*(v-p)/(v-s); F1 := beta*xi*q/(u*(1-alpha)); G1 := (1/2)*q^2/(u*alpha^2); G2 := (1/2)*beta^2*xi^2*q^2/(u*(1-alpha)^2); G3 := (1/2)*q^2*(1+beta*xi)^2/u; FirmpfSiS := (p-c)*q-alpha*(p-s)*G1+h*(p-delta*v)*(beta*xi*q-(1-alpha)*G2); RecpfSiS := ((1-h)*(p-delta*v)-sExp)*(beta*xi*q-G2)+sExp*xi*q-cr*xi*q; UnsoldSiS := alpha*G1; EnvironSiS := q+UnsoldSiS; OldSoldPrim := beta*xi*q-(1-alpha)*G2 end if; return p, q, FirmpfSiS, RecpfSiS, EnvironSiS, h, UnsoldSiS, OldSoldPrim, xi end proc

NULL

NULL

NULL

NULL

NULL

NULL

NULL

diffr1 := proc (alpha, delta) if not [alpha, delta]::(list(numeric)) then return ('procname')(args) end if; FirmModelPP(alpha, delta)[3]-FirmModelHmax(alpha, .2, delta)[3] end proc

diffr2 := proc (alpha, delta) if not [alpha, delta]::(list(numeric)) then return ('procname')(args) end if; FirmModelFC(alpha, delta)[3]-FirmModelHmax(alpha, .2, delta)[3] end proc

diffr3 := proc (alpha, delta) if not [alpha, delta]::(list(numeric)) then return ('procname')(args) end if; FirmModelFC(alpha, delta)[3]-FirmModelPP(alpha, delta)[3] end proc

P1 := implicitplot(diffr1, 0 .. 1, 0 .. 1, color = gray, thickness = 1)

P2 := implicitplot(diffr2, 0 .. 1, 0 .. 1, color = black, thickness = 1)

P3 := implicitplot(diffr3, 0 .. 1, 0 .. 1, color = pink, thickness = 1)

``

``

 

 

NULL

``

NULL

P7 := inequal({diffr1(alpha, delta) > 0, diffr2(alpha, delta) > 0, diffr3(alpha, delta) > 0}, alpha = 0 .. 1, delta = 0 .. 1, color = "LightBlue")

P8 := inequal({diffr1(alpha, delta) > 0, diffr2(alpha, delta) > 0, diffr3(alpha, delta) < 0}, alpha = 0 .. 1, delta = 0 .. 1, color = yellow)``

P9 := inequal({diffr1(alpha, delta) > 0, diffr3(alpha, delta) > 0, diffr2(alpha, delta) < 0}, alpha = 0 .. 1, delta = 0 .. 1, color = cyan)

P10 := inequal({diffr1(alpha, delta) > 0, diffr2(alpha, delta) < 0, diffr3(alpha, delta) < 0}, alpha = 0 .. 1, delta = 0 .. 1, color = grey)

P11 := inequal({diffr2(alpha, delta) > 0, diffr3(alpha, delta) > 0, diffr1(alpha, delta) < 0}, alpha = 0 .. 1, delta = 0 .. 1, color = green)

P12 := inequal({diffr2(alpha, delta) > 0, diffr1(alpha, delta) < 0, diffr3(alpha, delta) < 0}, alpha = 0 .. 1, delta = 0 .. 1, color = red)

P13 := inequal({diffr3(alpha, delta) > 0, diffr1(alpha, delta) < 0, diffr2(alpha, delta) < 0}, alpha = 0 .. 1, delta = 0 .. 1, color = pink)

P14 := inequal({diffr1(alpha, delta) < 0, diffr2(alpha, delta) < 0, diffr3(alpha, delta) < 0}, alpha = 0 .. 1, delta = 0 .. 1, color = coral)

NULL

display(P7, P8, P9, P10, P11, P12, P13, P14, textplot([.2, .9, "some text"]), scaling = constrained, view = [0 .. 1, 0 .. 1], labels = [alpha, delta])

 

``


 

Download Compare_three_regions.mw

the final plot does not look reasonable. Would you please guide me?

This worksheet contains an unnamed theorem on page 202 of David Wells's book The Penguin Dictionary of Curious and Interesting Geometry.

Somehow I have uploaded both its contents and (a) link(s) to it.

What Maple code can animate and display, in turn, each of the portrayed pursuit paths?

Pursuit_problem.mw

 

Consider a target point T which moves at constant speed along a straight line, and a moving point P which at all times moves directly towards T. If P starts anywhere on the outermost ellipse, and T starts from a focus of the outer ellipse, then P always captures T at the same point, the centre of the ellipse.

 

``

 

The concentric ellipses, whose shape depends on the relative velocities of T and P, are isochrones, and the curves of pursuit are their isoclinal trajectories

 

Download Pursuit_problem

Please how can I export a matrix from maple to matlab

A := Matrix(5, 5, {(1, 1) = -.841752461600000, (1, 2) = -71.7787800100000, (1, 3) = 0.701877157500000e-2, (1, 4) = 0.672783592200000e-2, (1, 5) = 0.646005759700000e-2, (2, 1) = -.877557898000000, (2, 2) = -100., (2, 3) = 0.701617590500000e-2, (2, 4) = 0.672545092900000e-2, (2, 5) = 0.645785863600000e-2, (3, 1) = -1.00426692300000, (3, 2) = -.677765381900000, (3, 3) = 0.700840041000000e-2, (3, 4) = 0.671830608400000e-2, (3, 5) = 0.645127072200000e-2, (4, 1) = -1.31039754100000, (4, 2) = -.820833777300000, (4, 3) = 0.699547946600000e-2, (4, 4) = 0.670643167900000e-2, (4, 5) = 0.644032067900000e-2, (5, 1) = -2.17574621300000, (5, 2) = -1.15959068100000, (5, 3) = 0.697746998200000e-2, (5, 4) = 0.668987785500000e-2, (5, 5) = 0.642505292600000e-2})

I tried

ExportMatrix(matlabData, A, target=MATLAB,format=rectangular,mode=ascii);

 but didn't work. 


 

NULL

restart:

Digits:=15;

15

A:=Matrix(4,4,[[-1,2,0,0],[2,-1,2,0],[0,2,-1,2],[0,0,1,-1]],datatype=float[8],storage=sparse);

Matrix(4, 4, {(1, 1) = -1.0000000000000000, (1, 2) = 2.0000000000000000, (1, 3) = 0., (1, 4) = 0., (2, 1) = 2.0000000000000000, (2, 2) = -1.0000000000000000, (2, 3) = 2.0000000000000000, (2, 4) = 0., (3, 1) = 0., (3, 2) = 2.0000000000000000, (3, 3) = -1.0000000000000000, (3, 4) = 2.0000000000000000, (4, 1) = 0., (4, 2) = 0., (4, 3) = 1.0000000000000000, (4, 4) = -1.0000000000000000})

b:=Vector(4,[1,0.,0,0],datatype=float[8]):

 

sol:=LinearAlgebra:-LinearSolve(A,b,method=SparseDirectMKL);

Error, invalid input: LinearAlgebra:-LinearSolve expects value for keyword parameter method to be of type identical(none,SparseLU,SparseDirect,SparseDirectMKL,SparseIterative,LU,QR,solve,hybrid,Cholesky,subs,modular), but received SparseDirectMKL

NULL


 

Download buglinearsolve.mw

I was pleased to see the description of SparseDirectMKL, but it is not implemented properly yet. It is a step in the right direction, so please make this available in the near future

 

Hello everyone. after runnig this maple code i get the error :" Error, (in Optimization:-NLPSolve) non-numeric result encountered" . How can fix it?

restart;
Digits := 20;
with(LinearAlgebra);
with(linalg);
with(Optimization);
with(Student[Calculus1]);
NULL;
n := 2;
z := 1;
a := 0;
b := 1;
m := 2;
NULL;
PS1 := (j, t) -> piecewise(j = 0, 1, t^(j + s[j]));
PI1 := (j, t) -> piecewise(j = 0, 1, t^(j + q[j]));
PH1 := (j, t) -> piecewise(j = 0, 1, t^(j + h[j]));
PL1 := (j, t) -> piecewise(j = 0, 1, t^(j + l[j]));
NULL;
PS := t -> local j; Transpose(convert([seq(PS1(j, t), j = 0 .. n - 1)], Matrix));
P_I := t -> local j; Transpose(convert([seq(PI1(j, t), j = 0 .. n - 1)], Matrix));
PH := t -> local j; Transpose(convert([seq(PH1(j, t), j = 0 .. n - 1)], Matrix));
PL := t -> local j; Transpose(convert([seq(PL1(j, t), j = 0 .. n - 1)], Matrix));
Warning, (in PS) `j` is implicitly declared local
Warning, (in P_I) `j` is implicitly declared local
Warning, (in PH) `j` is implicitly declared local
Warning, (in PL) `j` is implicitly declared local


NULL;
NULL;
B1 := (i, j) -> piecewise(i = j and j = 1, 1, i = 1 and 1 < j, 0, i = j and j = 2, 1, i = 2 and 2 < j, 0, j <= i, (i + j - 2)!/(2^(j - 1)*(j - 1)!*(i - j)!));
NULL;

B := t -> Matrix(n, B1);

DS1 := (m, k) -> piecewise(m <= 1, 0, m = k, GAMMA(k + s[k - 1])/GAMMA(k + s[k - 1] - z));
DI1 := (m, k) -> piecewise(m <= 1, 0, m = k, GAMMA(k + i[k - 1])/GAMMA(k + i[k - 1] - z));
DH1 := (m, k) -> piecewise(m <= 1, 0, m = k, GAMMA(k + h[k - 1])/GAMMA(k + h[k - 1] - z));
DL1 := (m, k) -> piecewise(m <= 1, 0, m = k, GAMMA(k + l[k - 1])/GAMMA(k + l[k - 1] - z));
DS := t -> t^(-z)*Matrix(n, DS1);
DI := t -> t^(-z)*Matrix(n, DI1);
DH := t -> t^(-z)*Matrix(n, DH1);
DL := t -> t^(-z)*Matrix(n, DL1);
NULL;
NULL;

cs := i -> ps[i];
ci := j -> p[j];
ch := j -> ph[j];
cl := i -> pl[i];
CS := convert([seq(cs(i), i = 1 .. n)], Matrix);
CI := convert([seq(ci(i), i = 1 .. n)], Matrix);
CH := convert([seq(ch(i), i = 1 .. n)], Matrix);
CL := convert([seq(cl(i), i = 1 .. n)], Matrix);

NULL;
NULL;
NULL;
S := unapply(simplify(Multiply(Multiply(CS, B(t)), PS(t))[1, 1]), [t]);
I1 := unapply(simplify(Multiply(Multiply(CI, B(t)), P_I(t))[1, 1]), [t]);
H := unapply(simplify(Multiply(Multiply(CH, B(t)), PH(t))[1, 1]), [t]);
L := unapply(simplify(Multiply(Multiply(CL, B(t)), PL(t))[1, 1]), [t]);
DSS := unapply(simplify(Multiply(Multiply(Multiply(CS, B(t)), DS(t)), PS(t))[1, 1]), [t]);
DII := unapply(simplify(Multiply(Multiply(Multiply(CI, B(t)), DI(t)), P_I(t))[1, 1]), [t]);
DHH := unapply(simplify(Multiply(Multiply(Multiply(CH, B(t)), DH(t)), PH(t))[1, 1]), [t]);
DLL := unapply(simplify(Multiply(Multiply(Multiply(CL, B(t)), DL(t)), PL(t))[1, 1]), [t]);
NULL;

NULL;
RS := unapply(evalf(DSS(t) + (-0.0043217^z + 0.5944^z*Multiply(S(t), I1(t)) + (0.025^z + 0.0008^z)*S(t))), [t]);

RI := unapply(evalf(DII(t) + (-0.5944^z*Multiply(S(t), I1(t)) - 0.0056^z*Multiply(H(t), I1(t)) - 0.027^z*L(t) + (((0.025^z + 0.0008^z) + 0.025^z) + 0.5^z)*I1(t))), [t]);
RH := unapply(evalf(DHH(t) + (-0.535^z + 0.0056^z*Multiply(H(t), I1(t)) - 0.5^z*I1(t) + (0.025^z + 0.0008^z)*H(t))), [t]);
RL := unapply(evalf(DLL(t) + (-0.025^z*I1(t) + (0.025^z + 0.0008^z + 0.027^z)*L(t))), [t]);
R := unapply(evalf(RS(t)^2 + RI(t)^2 + RH(t)^2 + RL(t)^2), [t]);

NULL;

p1 := x -> (x^2 - 1)^m;
dmp1 := x -> diff(p1(x), x $ m);
NULL;
p := x -> dmp1(x)/(2^m*m!);
eq := p(x) = 0;
r := solve(eq, x);
NULL;
ss := Vector[row](m);
w := Vector[row](m);
for i to m do
    w[i] := 2/((-r[i]^2 + 1)*D(p)(r[i])^2);
    ss[i] := w[i]*evalf(R((b - a)/2*r[i] + (b + a)/2));
end do;


SS := add(ss);

Lambda := evalf((b - a)/2*SS);
;

C1 := S(0) - 43994 = 0;
C2 := I1(0) - 0.1 = 0;
C3 := H(0) = 0;
C4 := L(0) - 1 = 0;

NULL;
NLP := NLPSolve(Lambda, {C1, C2, C3, C4});
Error, (in Optimization:-NLPSolve) non-numeric result encountered

 

First 9 10 11 12 13 14 15 Last Page 11 of 2263