Maple 2024 Questions and Posts

These are Posts and Questions associated with the product, Maple 2024

I'm new to Mapple 2024, and probably miss something.
lie, e[i] i=1..4, and G are previously defined objects. Then whe I define:

cf:=(k,l,m,t,x,y,z)->seq(lie(p,e[k](t,x,y,z),e[l](t,x,y,z))*G(t,x,y,z)[p,m],p=1..4);
cf(2,4,2,t,x,y,z);

I get:

0, 0, 0, -6*(-6*z^3*(7*z - 4)*(z - 1)^2*(y^2 - 1)^3*(x^2 - 1)^2*x - 6*z^4*(1 - z)^3*(-x^2 + 1)^2*(-y^2 + 1)^3*x*(-3*z^2*(7*z - 4)*(z - 1)^2*(y^2 - 1)^3*(x^2 - 1)^3 - 7*z^3*(z - 1)^2*(y^2 - 1)^3*(x^2 - 1)^3 - 2*z^3*(7*z - 4)*(z - 1)*(y^2 - 1)^3*(x^2 - 1)^3) - (1 - z^3*(7*z - 4)*(z - 1)^2*(y^2 - 1)^3*(x^2 - 1)^3)*(-24*z^3*(1 - z)^3*(-x^2 + 1)^2*(-y^2 + 1)^3*x + 18*z^4*(1 - z)^2*(-x^2 + 1)^2*(-y^2 + 1)^3*x))*z^4*(1 - z)^3*(-x^2 + 1)^2*(-y^2 + 1)^3*x*(1 - z^3*(7*z - 4)*(z - 1)^2*(y^2 - 1)^3*(x^2 - 1)^3)

The 4th item eval as:
cf(2, 4, 2, t, x, y, z)[4] = -36*z^10*x^2*(x^2 - 1)^7*(7*z^2 - 8*z + 4)*(z - 1)^7*(y^2 - 1)^9*(-1 + z^3*(7*z - 4)*(z - 1)^2*(y^2 - 1)^3*(x^2 - 1)^3)

But when I define another function by just changing the seq operator with a sum operator:

cf2:=(k,l,m,t,x,y,z)->sum(lie(p,e[k](t,x,y,z),e[l](t,x,y,z))*G(t,x,y,z)[p,m],p=1..4);
cf2(2,4,2,t,x,y,z);

I get 0.

What do I miss ?

Good day to all of you nice people.
I'm currently attempting to plot a vector field where each component of the vector is defined by the equations S_x, S_y, and S_z, which are functions of the radial coordinate. Here is a depiction of how the vectors change with respect to r:

The next step, which I'm unsure how to do, is to plot the vectors around the z-direction, or I should say, in phi direction, to achieve something similar to this example:

Thank you a lot for your kind help. 

Here is my code:
Maple_Question.mw

In our recent project, we're diving deep into understanding the SIR model—a fundamental framework in epidemiology that helps us analyze how diseases spread through populations. The SIR model categorizes individuals into three groups: Susceptible (S), Infected (I), and Recovered (R). By tracking how people move through these categories, we can predict disease dynamics and evaluate interventions.

Key Points of the SIR Model:

  • Susceptible (S): Individuals who can catch the disease.
  • Infected (I): Those currently infected and capable of spreading the disease.
  • Recovered (R): Individuals who have recovered and developed immunity.

Vaccination Impact: One of the critical interventions in disease control is vaccination, which moves individuals directly from the susceptible to the recovered group. This simple action reduces the number of people at risk, thereby lowering the overall spread of the disease.

We're experimenting with a simple model to understand how different vaccination rates can significantly alter the dynamics of an outbreak. By simulating scenarios with varying vaccination coverage, students can observe how herd immunity plays a crucial role in controlling diseases. Our goal is to make these abstract concepts clear and relatable through practical modeling exercises.


 

In this exercise, we are going back to the simple SIR model, without births or deaths, to look at the effect of vaccination. The aim of this activity is to represent vaccination in a very simple way - we are assuming it already happened before we run our model! By changing the initial conditions, we can prepare the population so that it has received a certain coverage of vaccination.

We are starting with the transmission and recovery parameters  b = .4/daysand c = .1/days . To incorporate immunity from vaccination in the model, we assume that a proportion p of the total population starts in the recovered compartment, representing the vaccine coverage and assuming the vaccine is perfectly effective. Again, we assume the epidemic starts with a single infected case introduced into the population.​
We are going to model this scenario for a duration of 2 years, assuming that the vaccine coverage is 50%, and plot the prevalence in each compartment over time.

 

restart
with(plots)

b := .4; c := .1; n := 10^6; p := .5

deS := diff(S(t), t) = -b*S(t)*I0(t); deI := diff(I0(t), t) = b*S(t)*I0(t)-c*I0(t); deR := diff(R(t), t) = c*I0(t)

diff(R(t), t) = .1*I0(t)

(1)

F := dsolve([deS, deI, deR, S(0) = 1-p, I0(0) = 1/n, R(0) = p], [S(t), I0(t), R(t)], numeric, method = rkf45, maxfun = 100000)

odeplot(F, [[t, S(t)], [t, I0(t)], [t, R(t)]], t = 0 .. 730, colour = [blue, red, green], legend = ["S(t)", "I0(t)", "R(t)"], labels = ["Time (days)", "  Proportion\nof Population "], title = "SIR Model with vaccine coverage 50 %", size = [500, 300])

 

F(100)

[t = 100., S(t) = HFloat(0.46146837378273076), I0(t) = HFloat(0.018483974421123688), R(t) = HFloat(0.5200486517961457)]

(2)

eval(S(:-t), F(100))

HFloat(0.46146837378273076)

(3)

Reff := proc (s) options operator, arrow; b*(eval(S(:-t), F(s)))/(c*n) end proc; Reff(100)

HFloat(1.845873495130923e-6)

(4)

plot(Reff, 0 .. 730, size = [500, 300])

 

Increasing the vaccine coverage to 75%

NULL

restart
with(plots)

b := .4; c := .1; n := 10^6; p := .75

deS := diff(S(t), t) = -b*S(t)*I0(t); deI := diff(I0(t), t) = b*S(t)*I0(t)-c*I0(t); deR := diff(R(t), t) = c*I0(t)

diff(R(t), t) = .1*I0(t)

(5)

NULL

F1 := dsolve([deS, deI, deR, S(0) = 1-p, I0(0) = 1/n, R(0) = p], [S(t), I0(t), R(t)], numeric, method = rkf45, maxfun = 100000)

odeplot(F1, [[t, S(t)], [t, I0(t)], [t, R(t)]], t = 0 .. 730, colour = [blue, red, green], legend = ["S(t)", "I0(t)", "R(t)"], labels = ["Time (days)", "  Proportion\nof Population "], title = "SIR Model with vaccine coverage 75%", size = [500, 300])

 

F(1100)

eval(S(:-t), F1(100))

HFloat(0.249990000844159)

(6)

Reff := proc (s) options operator, arrow; b*(eval(S(:-t), F1(s)))/(c*n) end proc; Reff(100)

HFloat(9.99960003376636e-7)

(7)

plot(Reff, 0 .. 730, size = [500, 300])

 

Does everyone in the population need to be vaccinated in order to prevent an epidemic?What do you observe if you model the infection dynamics with different values for p?

No, not everyone in the population needs to be vaccinated in order to prevent an epidemic . In this scenario, if p equals 0.75 or higher, no epidemic occurs - 75 % is the critical vaccination/herd immunity threshold . Remember,, herd immunity describes the phenomenon in which there is sufficient immunity in a population to interrupt transmission . Because of this, not everyone needs to be vaccinated to prevent an outbreak .

What proportion of the population needs to be vaccinated in order to prevent an epidemic if b = .4and c = .2/days? What if b = .6 and "c=0.1 days^(-1)?"

In the context of the SIR model, the critical proportion of the population that needs to be vaccinated in order to prevent an epidemic is often referred to as the "herd immunity threshold" or "critical vaccination coverage."

• 

Scenario 1: b = .4and c = .2/days

``

restart
with(plots)

b := .4; c := .2; n := 10^6; p := .5``

deS := diff(S(t), t) = -b*S(t)*I0(t); deI := diff(I0(t), t) = b*S(t)*I0(t)-c*I0(t); deR := diff(R(t), t) = c*I0(t)

diff(R(t), t) = .2*I0(t)

(8)

F1 := dsolve([deS, deI, deR, S(0) = 1-p, I0(0) = 1/n, R(0) = p], [S(t), I0(t), R(t)], numeric, method = rkf45, maxfun = 100000)

odeplot(F1, [[t, S(t)], [t, I0(t)], [t, R(t)]], t = 0 .. 730, colour = [blue, red, green], legend = ["S(t)", "I0(t)", "R(t)"], labels = ["Time (days)", "  Proportion\nof Population "], title = "SIR Model with vaccine coverage 50 %", size = [500, 300])

 


The required vaccination coverage is around 50% .

• 

Scenario 1: b = .6and c = .1/days

restart
with(plots)

b := .6; c := .1; n := 10^6; p := .83NULL

deS := diff(S(t), t) = -b*S(t)*I0(t); deI := diff(I0(t), t) = b*S(t)*I0(t)-c*I0(t); deR := diff(R(t), t) = c*I0(t)

diff(R(t), t) = .1*I0(t)

(9)

NULL

F1 := dsolve([deS, deI, deR, S(0) = 1-p, I0(0) = 1/n, R(0) = p], [S(t), I0(t), R(t)], numeric, method = rkf45, maxfun = 100000)

odeplot(F1, [[t, S(t)], [t, I0(t)], [t, R(t)]], t = 0 .. 730, colour = [blue, red, green], legend = ["S(t)", "I0(t)", "R(t)"], labels = ["Time (days)", "  Proportion\nof Population "], title = "SIR Model with vaccine coverage 83% ", size = [500, 300])

 

"The required vaccination coverage is around 83 `%` ."


Download SIR_simple_vaccination_example.mw

I have a complicated bivariate function f(Gamma,rho) that is a RootOf of a quartic. I know that it is strictly positive (one of the four roots at least) for Gamma=0..10 and rho in (-1,+1), with bounds excluded.

I need to find the signs of its first and second derivatives (wrt to Gamma and wrt to rho: 4 derivatives in total).

I encounter numerical issues when I plot3d the derivatives using D[]() vs. fdiff() (numerical function evaluations of the RootOf). I was hoping for the two commands to produce the same output, but they don't it seems. What's going on?

Script:

restart;
_quartic := RootOf(-8*(rho + 1)^4*_Z^4 + 12*(rho + 1)^3*Gamma*(rho - 1)*_Z^3 - 5*(rho + 1)^2*(-4/5 + Gamma^2*rho^2 + 2*(-2/5 - Gamma^2)*rho + Gamma^2)*_Z^2 - 4*(rho + 1)*Gamma*(rho^2 - 1)*_Z + Gamma^2*(rho + 1)*(rho - 1)^2);
convert(_quartic,radical):
f(Gamma,rho) := simplify(%):

RootOf((8*rho^3+24*rho^2+24*rho+8)*_Z^4+(-12*Gamma*rho^3-12*Gamma*rho^2+12*Gamma*rho+12*Gamma)*_Z^3+(5*Gamma^2*rho^3-5*Gamma^2*rho^2-5*Gamma^2*rho+5*Gamma^2-4*rho^2-8*rho-4)*_Z^2+(4*Gamma*rho^2-4*Gamma)*_Z-Gamma^2*rho^2+2*rho*Gamma^2-Gamma^2)

(1)


Synthetic representation of derivatives

der1_Gamma := diff(_quartic, Gamma):
der1_rho := diff(_quartic, rho):

Diff('f(Gamma,rho)', Gamma) = collect~(normal(eval(der1_Gamma, _quartic = 'f(Gamma,rho)')), 'f(Gamma,rho)');
Diff('f(Gamma,rho)', rho) = collect~(normal(eval(der1_rho, _quartic = 'f(Gamma,rho)')), 'f(Gamma,rho)');

der2_Gamma := diff(der1_Gamma, Gamma):
der2_rho := diff(der1_rho, rho):

Diff('f(Gamma,rho)', Gamma$2) = collect~(normal(eval(der2_Gamma, _quartic = 'f(Gamma,rho)')), 'f(Gamma,rho)');
Diff('f(Gamma,rho)', rho$2) = collect~(normal(eval(der2_rho, _quartic = 'f(Gamma,rho)')), 'f(Gamma,rho)');

Diff(f(Gamma, rho), Gamma) = -((-6*rho^3-6*rho^2+6*rho+6)*f(Gamma, rho)^3+(5*Gamma*rho^3-5*Gamma*rho^2-5*Gamma*rho+5*Gamma)*f(Gamma, rho)^2+(2*rho^2-2)*f(Gamma, rho)-Gamma*rho^2+2*Gamma*rho-Gamma)/((16*rho^3+48*rho^2+48*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^3-18*Gamma*rho^2+18*Gamma*rho+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^3-5*Gamma^2*rho^2-5*Gamma^2*rho+5*Gamma^2-4*rho^2-8*rho-4)*f(Gamma, rho)+2*Gamma*rho^2-2*Gamma)

 

Diff(f(Gamma, rho), rho) = -(1/2)*((24*rho^2+48*rho+24)*f(Gamma, rho)^4+(-36*Gamma*rho^2-24*Gamma*rho+12*Gamma)*f(Gamma, rho)^3+(15*Gamma^2*rho^2-10*Gamma^2*rho-5*Gamma^2-8*rho-8)*f(Gamma, rho)^2+8*Gamma*rho*f(Gamma, rho)-2*rho*Gamma^2+2*Gamma^2)/((16*rho^3+48*rho^2+48*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^3-18*Gamma*rho^2+18*Gamma*rho+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^3-5*Gamma^2*rho^2-5*Gamma^2*rho+5*Gamma^2-4*rho^2-8*rho-4)*f(Gamma, rho)+2*Gamma*rho^2-2*Gamma)

 

Diff(f(Gamma, rho), Gamma, Gamma) = ((448*rho^8+1792*rho^7+1792*rho^6-1792*rho^5-4480*rho^4-1792*rho^3+1792*rho^2+1792*rho+448)*f(Gamma, rho)^8+(-1632*Gamma*rho^8-3264*Gamma*rho^7+3264*Gamma*rho^6+9792*Gamma*rho^5-9792*Gamma*rho^3-3264*Gamma*rho^2+3264*Gamma*rho+1632*Gamma)*f(Gamma, rho)^7+(2120*Gamma^2*rho^8-8480*Gamma^2*rho^6-208*rho^7+12720*Gamma^2*rho^4-624*rho^6-208*rho^5-8480*Gamma^2*rho^2+1040*rho^4+1040*rho^3+2120*Gamma^2-208*rho^2-624*rho-208)*f(Gamma, rho)^6+(-1200*Gamma^3*rho^8+2400*Gamma^3*rho^7+2400*Gamma^3*rho^6-7200*Gamma^3*rho^5+640*Gamma*rho^7+640*Gamma*rho^6+7200*Gamma^3*rho^3-1920*Gamma*rho^5-2400*Gamma^3*rho^2-1920*Gamma*rho^4-2400*Gamma^3*rho+1920*Gamma*rho^3+1200*Gamma^3+1920*Gamma*rho^2-640*Gamma*rho-640*Gamma)*f(Gamma, rho)^5+(250*Gamma^4*rho^8-1000*Gamma^4*rho^7+1000*Gamma^4*rho^6+1000*Gamma^4*rho^5-632*Gamma^2*rho^7-2500*Gamma^4*rho^4+632*Gamma^2*rho^6+1000*Gamma^4*rho^3+1896*Gamma^2*rho^5+1000*Gamma^4*rho^2-1896*Gamma^2*rho^4+16*rho^6-1000*Gamma^4*rho-1896*Gamma^2*rho^3+32*rho^5+250*Gamma^4+1896*Gamma^2*rho^2-16*rho^4+632*Gamma^2*rho-64*rho^3-632*Gamma^2-16*rho^2+32*rho+16)*f(Gamma, rho)^4+(240*Gamma^3*rho^7-720*Gamma^3*rho^6+240*Gamma^3*rho^5+1200*Gamma^3*rho^4-32*Gamma*rho^6-1200*Gamma^3*rho^3-240*Gamma^3*rho^2+96*Gamma*rho^4+720*Gamma^3*rho-240*Gamma^3-96*Gamma*rho^2+32*Gamma)*f(Gamma, rho)^3+(-25*Gamma^4*rho^7+125*Gamma^4*rho^6-225*Gamma^4*rho^5+125*Gamma^4*rho^4+125*Gamma^4*rho^3-225*Gamma^4*rho^2+125*Gamma^4*rho-25*Gamma^4)*f(Gamma, rho)^2+(16*Gamma^3*rho^6-64*Gamma^3*rho^5+80*Gamma^3*rho^4-80*Gamma^3*rho^2+64*Gamma^3*rho-16*Gamma^3)*f(Gamma, rho)-5*Gamma^4*rho^6+30*Gamma^4*rho^5-75*Gamma^4*rho^4+100*Gamma^4*rho^3-75*Gamma^4*rho^2+30*Gamma^4*rho-5*Gamma^4)/(((16*rho^2+32*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^2+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^2-10*Gamma^2*rho+5*Gamma^2-4*rho-4)*f(Gamma, rho)+2*Gamma*rho-2*Gamma)*((16*rho^3+48*rho^2+48*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^3-18*Gamma*rho^2+18*Gamma*rho+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^3-5*Gamma^2*rho^2-5*Gamma^2*rho+5*Gamma^2-4*rho^2-8*rho-4)*f(Gamma, rho)+2*Gamma*rho^2-2*Gamma)^2)

 

Diff(f(Gamma, rho), rho, rho) = (1/4)*((21504*rho^6+129024*rho^5+322560*rho^4+430080*rho^3+322560*rho^2+129024*rho+21504)*f(Gamma, rho)^10+(-80640*Gamma*rho^6-347136*Gamma*rho^5-526080*Gamma*rho^4-245760*Gamma*rho^3+157440*Gamma*rho^2+199680*Gamma*rho+56064*Gamma)*f(Gamma, rho)^9+(127680*Gamma^2*rho^6+336512*Gamma^2*rho^5+122944*Gamma^2*rho^4-319744*Gamma^2*rho^3-18176*rho^5-246976*Gamma^2*rho^2-90880*rho^4+40576*Gamma^2*rho-181760*rho^3+53696*Gamma^2-181760*rho^2-90880*rho-18176)*f(Gamma, rho)^8+(-110016*Gamma^3*rho^6-112128*Gamma^3*rho^5+191808*Gamma^3*rho^4+172032*Gamma^3*rho^3+57344*Gamma*rho^5-105792*Gamma^3*rho^2+191488*Gamma*rho^4-59904*Gamma^3*rho+192512*Gamma*rho^3+24000*Gamma^3+2048*Gamma*rho^2-94208*Gamma*rho-37888*Gamma)*f(Gamma, rho)^7+(54960*Gamma^4*rho^6-28480*Gamma^4*rho^5-102480*Gamma^4*rho^4+56960*Gamma^4*rho^3-74176*Gamma^2*rho^5+40080*Gamma^4*rho^2-126144*Gamma^2*rho^4-28480*Gamma^4*rho+41088*Gamma^2*rho^3+7440*Gamma^4+138368*Gamma^2*rho^2+5120*rho^4+19776*Gamma^2*rho+20480*rho^3-25536*Gamma^2+30720*rho^2+20480*rho+5120)*f(Gamma, rho)^6+(-15300*Gamma^5*rho^6+30000*Gamma^5*rho^5-2100*Gamma^5*rho^4-24000*Gamma^5*rho^3+50048*Gamma^3*rho^5+14100*Gamma^5*rho^2+5632*Gamma^3*rho^4-6000*Gamma^5*rho-89856*Gamma^3*rho^3+3300*Gamma^5-1024*Gamma^3*rho^2-13056*Gamma*rho^4+39808*Gamma^3*rho-31232*Gamma*rho^3-4608*Gamma^3-15360*Gamma*rho^2+10752*Gamma*rho+7936*Gamma)*f(Gamma, rho)^5+(1875*Gamma^6*rho^6-6250*Gamma^6*rho^5+7125*Gamma^6*rho^4-3500*Gamma^6*rho^3-18316*Gamma^4*rho^5+2125*Gamma^6*rho^2+25540*Gamma^4*rho^4-2250*Gamma^6*rho+12072*Gamma^4*rho^3+875*Gamma^6-26520*Gamma^4*rho^2+12992*Gamma^2*rho^4+6244*Gamma^4*rho+10240*Gamma^2*rho^3+980*Gamma^4-15488*Gamma^2*rho^2-9728*Gamma^2*rho-512*rho^3+3008*Gamma^2-1536*rho^2-1536*rho-512)*f(Gamma, rho)^4+(3320*Gamma^5*rho^5-9240*Gamma^5*rho^4+7600*Gamma^5*rho^3-560*Gamma^5*rho^2-6144*Gamma^3*rho^4-1320*Gamma^5*rho+4992*Gamma^3*rho^3+200*Gamma^5+6784*Gamma^3*rho^2-4992*Gamma^3*rho+1024*Gamma*rho^3-640*Gamma^3+1536*Gamma*rho^2-512*Gamma)*f(Gamma, rho)^3+(-200*Gamma^6*rho^5+800*Gamma^6*rho^4-1200*Gamma^6*rho^3+800*Gamma^6*rho^2+1248*Gamma^4*rho^4-200*Gamma^6*rho-3136*Gamma^4*rho^3+1920*Gamma^4*rho^2+576*Gamma^4*rho-768*Gamma^2*rho^3-608*Gamma^4+768*Gamma^2*rho)*f(Gamma, rho)^2+(-16*Gamma^5*rho^4+192*Gamma^5*rho^3-480*Gamma^5*rho^2+448*Gamma^5*rho+256*Gamma^3*rho^3-144*Gamma^5-384*Gamma^3*rho^2+128*Gamma^3)*f(Gamma, rho)-20*Gamma^6*rho^4+80*Gamma^6*rho^3-120*Gamma^6*rho^2+80*Gamma^6*rho-32*Gamma^4*rho^3-20*Gamma^6+96*Gamma^4*rho^2-96*Gamma^4*rho+32*Gamma^4)/(((16*rho^2+32*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^2+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^2-10*Gamma^2*rho+5*Gamma^2-4*rho-4)*f(Gamma, rho)+2*Gamma*rho-2*Gamma)*((16*rho^3+48*rho^2+48*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^3-18*Gamma*rho^2+18*Gamma*rho+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^3-5*Gamma^2*rho^2-5*Gamma^2*rho+5*Gamma^2-4*rho^2-8*rho-4)*f(Gamma, rho)+2*Gamma*rho^2-2*Gamma)^2)

(2)


Signs of derivatives: fdiff (numerical function evaluations of the RootOf) vs. D[]()

restart;
with(plots):

_quartic := RootOf(-8*(rho + 1)^4*_Z^4 + 12*(rho + 1)^3*Gamma*(rho - 1)*_Z^3 - 5*(rho + 1)^2*(-4/5 + Gamma^2*rho^2 + 2*(-2/5 - Gamma^2)*rho + Gamma^2)*_Z^2 - 4*(rho + 1)*Gamma*(rho^2 - 1)*_Z + Gamma^2*(rho + 1)*(rho - 1)^2):

plot3d(_quartic, Gamma=0..10, rho=-1..+1, labels=[Gamma,rho,Lambda(Gamma,rho)],axesfont=["helvetica","roman",20],labelfont=["helvetica","roman",30]);
 

 

Define it as a f and test it for Gamma=1 and rho=0.5

f := (Gamma,rho) -> RootOf(-8*(rho + 1)^4*_Z^4 + 12*(rho + 1)^3*Gamma*(rho - 1)*_Z^3 - 5*(rho + 1)^2*(-4/5 + Gamma^2*rho^2 + 2*(-2/5 - Gamma^2)*rho + Gamma^2)*_Z^2 - 4*(rho + 1)*Gamma*(rho^2 - 1)*_Z + Gamma^2*(rho + 1)*(rho - 1)^2):
evalf(f(1.0,0.5));

HFloat(0.5110796212870378)

(3)

Value at zero:

f(0,0):
allvalues(%):
fl := select(is, [allvalues(f(0,0))], positive)[];evalf(%);

(1/2)*2^(1/2)

 

.7071067810

(4)

Value at infinity (commented out because too slow)

#limit(f(x,y), {x = infinity, y = 0}):
#fh := select(is, [allvalues(%)], positive)[];evalf(%);

Derivative at zero:

allvalues([D[1](f)(0,0)]):
Dfl := %[1][];

-1/4

(5)

Derivative at a point, evaluated, vs numerical derivative at a point:

D[1](f)(1,0.5):
evalf(%);
fdiff(f(x,y), x, {x = 1.0, y = 0.5});
fdiff(f, [1], [1.0,0.5]);

D[2](f)(1,0.5):
evalf(%);
fdiff(f(x,y), y, {x = 1.0, y = 0.5});
fdiff(f, [2], [1.0,0.5]);

HFloat(-0.05086932918910799)

 

-0.5086932919e-1

 

-0.5086932919e-1

 

HFloat(-0.05166477232109392)

 

-0.5166477232e-1

 

-0.5166477232e-1

(6)

Can make a function out of fdiff

fDfG := (Gamma,rho) -> fdiff(f, [1], [Gamma,rho]);
fDfr := (Gamma,rho) -> fdiff(f, [2], [Gamma,rho]);

proc (Gamma, rho) options operator, arrow; fdiff(f, [1], [Gamma, rho]) end proc

 

proc (Gamma, rho) options operator, arrow; fdiff(f, [2], [Gamma, rho]) end proc

(7)

Check for numerical values close to thresholds:

Digits := 15:
evalf('D[1]'(f)(0.1e-8,0.5));fdiff(f, [1], [0.1e-8,0.5]);
evalf('D[1]'(f)(0.1e-7,0.5));fdiff(f, [1], [0.1e-7,0.5]);
evalf('D[1]'(f)(0.1e-5,0.5));fdiff(f, [1], [0.1e-5,0.5]);
evalf('D[1]'(f)(0.00001,0.5));fdiff(f, [1], [0.00001,0.5]);
evalf('D[1]'(f)(0.001,0.5));fdiff(f, [1], [0.001,0.5]);


evalf('D[2]'(f)(1,-0.99));fdiff(f, [2], [1,-0.99]);
evalf('D[2]'(f)(1,-0.97));fdiff(f, [2], [1,-0.97]);
evalf('D[2]'(f)(1,-0.1));fdiff(f, [2], [1,-0.1]);
evalf('D[2]'(f)(1,0.98));fdiff(f, [2], [1,0.98]);
evalf('D[2]'(f)(1,-0.99));fdiff(f, [2], [1,-0.99]);

57735026.8022959

 

57735026.8022959

 

-0.833333329724894e-1

 

-0.833333329724894e-1

 

-0.833332972489415e-1

 

-0.833332972489415e-1

 

-0.833329724894151e-1

 

-0.833329724894151e-1

 

-0.832972489466445e-1

 

-0.832972489466445e-1

 

-223.615892086941

 

-223.615892086941

 

-43.0236130145893

 

-43.0236130145893

 

-.212392503268663

 

-.212392503268663

 

-0.127828146340716e-2

 

-0.127828146340716e-2

 

-223.615892086941

 

-223.615892086941

(8)

Compare with D (vertical range here to prevent effect of large values from fdiff near zero):

d1G := plot3d([D[1](f), fDfG], 0..10, -0.95..+0.95, view=-0.3..0, color = [red, blue]);
d1r := plot3d([D[2](f), fDfr], 0..10, -0.95..+0.95, color = [red, blue]);

 

 

 

Second derivatives:

evalf('D[1,1]'(f)(1.0,0.5));
fdiff(f, [1, 1], [1.0,0.5]);

evalf('D[2,2]'(f)(1.0,0.5));
fdiff(f, [2, 2], [1.0,0.5]);

fD2fG := (Gamma,rho) -> fdiff(f, [1, 1], [Gamma]);
fD2fr := (Gamma,rho) -> fdiff(f, [2, 2], [Gamma]);

0.266607527050519e-1

 

0.266607527050519e-1

 

.151600577769391

 

.151600577769391

 

proc (Gamma, rho) options operator, arrow; fdiff(f, [1, 1], [Gamma]) end proc

 

proc (Gamma, rho) options operator, arrow; fdiff(f, [2, 2], [Gamma]) end proc

(9)

d2G:= plot3d([D[1,1](f), fD2fG], 0..10, -0.9..+0.9, color = [red, blue]);
d2r:= plot3d([D[2,2](f), fD2fr], 0..10, -0.9..+0.9, color = [red, blue]);
 

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

 

 

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

 

 

d1d2G := plot3d([fDfG, fD2fG], 0.1e-6 .. 10, -0.98 .. +0.98, axesfont=["helvetica","roman",20],labelfont=["helvetica","roman",30], size=[1000,1000]);
d1d2r := plot3d([fDfr, fD2fr], 0.1e-6 .. 10, -0.98 .. +0.98, axesfont=["helvetica","roman",20],labelfont=["helvetica","roman",30], size=[1000,1000]);

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

 

 

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

 

 
 

NULL

Download signs_derivatves_bivariate.mw

When I calculate the difference between two dates with the same day using the DateDifference Calendar function, I don't get an integer number of months, but results with extra days, hours, minutes, seconds and milliseconds. How can I calulate the actual number of months?

I appreciate any help.

restart

with(Calendar)

d1 := Date(2024, 1, 5, 0, 0, 0)

_m4674925856

(1)

d2 := Date(2024, 5, 5, 0, 0, 0)

_m4881744064

(2)

DateDifference(d1, d2, 'units' = 'mixed')

4*Units:-Unit(mo)+2*Units:-Unit(d)+20*Units:-Unit(h)+3*Units:-Unit(min)+48*Units:-Unit(s)+800*Units:-Unit(ms)

(3)

NULL

NULL

 

NULL

Download DateDifference.mw

I never used Maple's StringTools:-RegSubs before.

I have a Latex string generated by Maple that has this form

"......    \\textrm{ ......   }  .... "

Where the dots mean anything, including space characters.  I need to change the part of the string  "\\textrm{ ...... }" to this

"......  \\begin{minipage}{\\linewidth}\\textrm{  ......   }\\end{minipage}  .... "

It is easy to change the opening part. Here is an example

s:="A& B& \\textrm{The fundamental matrix has } &C ";
s:=StringTools:-SubstituteAll(s,"\\textrm{","\\begin{minipage}{\\linewidth}\\textrm{");

Gives

The problem is how to replace the closing "}"  with "}\\end{minipage}".   I can't just replace "}" with "}\\end{minipage}" using SubstituteAll because I need to only change the "}" that closes the corresponding "\\textrm{" and not any other "}" that could be in the string.

So I need to use pattern matching or regular expression substitutions. The examples in help are not easy to understand.

 

Basically I need regular expression that matches 

                        "\\textrm{ WILDCARD * }" 

and change it to 

                        "\\begin{minipage}{\\linewidth}\\textrm{ WILDCARD * }\\end{minipage}" 

Any one know how to use Maple's egular expression to do this substitution using the above example?

Maple 2024 on windows 10.

Update

Thanks for the replies. I ended up writing small function that simply scan the string and do the replacement. 

Can this be better done in Maple ? , see worksheet.

In this exercise, we are going back to the simple SIR model, without births or deaths, to look at the effect of vaccination. The aim of this activity is to represent vaccination in a very simple way - we are assuming it already happened before we run our model! By changing the initial conditions, we can prepare the population so that it has received a certain coverage of vaccination.

We are starting with the transmission and recovery parameters b = 0.4 days-1 and  c=0.1 days-1. To incorporate immunity from vaccination in the model, we assume that a proportion p of the total population starts in the recovered compartment, representing the vaccine coverage and assuming the vaccine is perfectly effective. Again, we assume the epidemic starts with a single infected case introduced into the population.​
We are going to model this scenario for a duration of 2 years, assuming that the vaccine coverage is 50%, and plot the prevalence in each compartment over time.

I have used the following code to solve this problem but I face difficulties when I am trying to define and plot the effective reproduction number.

Any suggestions? Thank you in advance!
with(Physics, diff);
with(LinearAlgebra);

with(plots);
with(DEtools);

b := 0.4;
c := 0.1;
n := 10^6;
p := 0.5;
deS := diff(S(t), t) = -b*S(t)*I0(t);
deI := diff(I0(t), t) = b*S(t)*I0(t) - c*I0(t);
deR := diff(R(t), t) = c*I0(t);

F := dsolve([deS, deI, deR, S(0) = 1 - p, I0(0) = 1/n, R(0) = p], [S(t), I0(t), R(t)], numeric, method = rkf45, maxfun = 100000)

odeplot(F, [[t, S(t)], [t, I0(t)], [t, R(t)]], t = 0 .. 730, colour = [blue, red, green], legend = ["S(t)", "I0(t)", "R(t)"], labels = ["Time (days)", "Proportion of Population"], title = "SIR Model with vaccination")


#define the effective reproduction number
Reff := t -> b*1/c*S(t)*1/n;
Reff(100);

SI_MODEL.mw
I have tried the following code to simulate the following SI Model. But when I run it, it gives me several mistakes and most importantly does not calculate the Jacobian Matrix and the eigenvalues. 

Any help would be deeply appreciated


Many of the most interesting dynamics in nature have to do with interactions between organisms. These interactions are often subtle, indirect, and difficult to detect. Interactions in which one organism consumes all or part of another. This includes predator-prey, herbivore-plant, and parasite-host interactions. These linkages are the prime movers of energy through food chains. They are an important factor in the ecology of populations, determining the mortality of prey and the birth of new predators.

Mathematical models and logic suggest that a coupled system of predator and prey should cycle: predators increase when prey is abundant, prey are driven to low numbers by predation, the predators decline, and the prey recover, ad infinitum. One such model that simulates predator-prey interactions is the Lotka - Volterra Model.

 


We will discuss a behavior of 2 D - subsystems: SI - model.


SI Model (without predator)

 

We will study a mathematical model that appears in eco - eco-epidemiology.


The SI model describes the interactions between S-prey with density x and I-prey with density y under the assumption z≡0


Break the population into compartments.

• 

Prey (S)

• 

Infected prey (I)


SI Equations :

(D(x))(t) = r*x(t)-b*x(t)^2-c*x(t)*y(t)-beta*x(t)*y[t]/(a+x(t))

(D(y))(t) = -mu*y(t)+beta*x(t)*y(t)/(a+x(t))

– 

 r stands for the intrinsic growth rate of S-prey.

– 

b defines the intra-class competition in S-prey.

– 

c characterizes the inter-class competition between S-prey and I-prey.

– 

μ stands for the mortality of the I-prey

– 

β considered as the bifurcation parameter.

• 

Note: The dynamics of the system depend on parameter β

 

Modeling with Differential Equations

 

• 

We will fix the parameters in the study to have the following values r = 1, b = 1, c = 1/100, μ = 4/10, a=1/2.75.

• 

The initial values will be x(0)=0.2, y(0)=0.05


As mentioned above the dynamics of the system depend on parameter β. So we will consider two cases when β=0.75 and β=1


We will solve the system for the above values of parameters and initial conditions and for a time interval [0,200] , then we will create its plots

 

Case 1

  restart

beta := .75; r := 1; b := 1; c := 1/100; mu := 4*(1/10); a := 1/2.75; f := r*x(t)-b*x(t)^2-c*x(t)*y(t)-(3/4)*x(t)*y(t)/(a+x(t)); g := -mu*y(t)+(3/4)*x(t)*y(t)/(a+x(t)); deq1 := diff(x(t), t) = f; deq2 := diff(y(t), t) = g; equilibrio := solve({f = 0, g = 0}, {x(t), y(t)}); jacobian := simplify(Matrix(2, 2, [[diff(f, x(t)), diff(f, y(t))], [diff(g, x(t)), diff(g, y(t))]])); A := simplify(eval(jacobian, equilibrio[1])); B := simplify(eval(jacobian, equilibrio[2])); Q := simplify(eval(jacobian, equilibrio[3])); eigen1 := Eigenvalues(A); eigen2 := Eigenvalues(B); eigen3 := Eigenvalues(Q); soln := dsolve({deq1, deq2, x(0) = .2, y(0) = 0.5e-1}, numeric, output = listprocedure); plots:-display(plots:-odeplot(soln, [t, x(t)], 0 .. 200, color = blue, legend = ["x(t)"]), plots:-odeplot(soln, [t, y(t)], 0 .. 200, color = red, legend = ["y(t)"]), labels = ["t", "Population"], title = "Population Dynamics")

 

NULL

 

 

Case 2

restart``

beta := 1; r := 1; b := 1; c := 1/100; Mu := 4*(1/10); a := 1/2.75; f := r*x(t)-b*x(t)^2-c*x(t)*y(t)-x(t)*y(t)/(a+x(t)); g := -Mu*y(t)+x(t)*y(t)/(a+x(t)); deq1 := diff(x(t), t) = f; deq2 := diff(y(t), t) = g; equilibrio := solve({f = 0, g = 0}, {x(t), y(t)}); jacobian := Matrix(2, 2, [[diff(f, x(t)), diff(f, y(t))], [diff(g, x(t)), diff(g, y(t))]], simplify); J_at_equilibrio := [eval(jacobian, equilibrio)]; eigen1 := Eigenvalues(A); eigen2 := Eigenvalues(B); eigen3 := Eigenvalues(Q); soln := dsolve({deq1, deq2, x(0) = .2, y(0) = 0.5e-1}, numeric); plots:-display(plots:-odeplot(soln, [t, x(t)], 0 .. 200, color = blue, legend = ["x(t)"]), plots:-odeplot(soln, [t, y(t)], 0 .. 200, color = red, legend = ["y(t)"]), labels = ["t", "Population"], title = "Population Dynamics with Beta = 1")

 

NULL


 

Download SI_MODEL.mw

 

 

The CauchyRiemann procedure (for older version  of Maple )doesn't work quite right in Maple 2024 .
Also ran the procedure through the AI for so-called code improvement and now it shows what the code stands for 
The output according to the original procedure would look like on the screenshot, but running original procedure does not give this output ? 
I also want to extend the procedure with a plot of the complex function. 
That differentiability of complex functions is not obvious even if the cauchy-riemann equation is satisfied ?

 

restart

"maple.ini in users"

(1)

NULL

CauchyRiemann:=proc(expr::algebraic) # original procedure
  local x, y, u, v, u_x, u_y, v_x, v_y, flag1, flag2;

  u:=evalc(Re(eval(expr, z=x+I*y)));
  v:=evalc(Im(eval(expr, z=x+I*y)));

  u_x:=diff(u,x);
  u_y:=diff(u,y);
  v_x:=diff(v,x);
  v_y:=diff(v,y);

  print('f(z)'=expr);
  printf("\n");
  
  print('u(x,y)'=u);
  print('u[x](x,y)'=u_x);
  print('u[y](x,y)'=u_y);
  printf("\n");

  print('v(x,y)'=v);
  print('v[x](x,y)'=v_x);
  print('v[y](x,y)'=v_y);
  printf("\n");

  if u_x=v_y then
    print('u[x]=v[y]');
    print(u_x=v_y);
    flag1:=true;
  else
    print('u[x]<>v[y]');
    print(u_x<>v_y);
    flag1:=false;
  end if;

  if u_y=-v_x then
    print('u[y]=-v[x]');
    print(u_y=-v_x);
    flag2:=true;
  else
    print('u[y]<>-v[x]');
    print(u_y<>-v_x);
    flag2:=false;
  end if;
  
printf("\n");
if flag1=true and flag2=true then
   print(`Fullfill the Cauchy-Riemann Equations`);
   print(`The derivative is:`='u[x]+I*v[y]');
   print('diff(f(z),z)'=u_x+I*v_y);
else
   print(`Cauchy-Riemann ?`);
end if

end proc:

f(z):=1/(z+2):
CauchyRiemann(f(z))

f(z) = 1/(z+2)

 

 

 

u(x, y) = (x+2)/(y^2+(x+2)^2)

 

u[x](x, y) = 1/(y^2+(x+2)^2)-(x+2)*(2*x+4)/(y^2+(x+2)^2)^2

 

u[y](x, y) = -2*(x+2)*y/(y^2+(x+2)^2)^2

 

 

 

v(x, y) = -y/(y^2+(x+2)^2)

 

v[x](x, y) = y*(2*x+4)/(y^2+(x+2)^2)^2

 

v[y](x, y) = -1/(y^2+(x+2)^2)+2*y^2/(y^2+(x+2)^2)^2

 

 

 

u[x] <> v[y]

 

1/(y^2+(x+2)^2)-(x+2)*(2*x+4)/(y^2+(x+2)^2)^2 <> -1/(y^2+(x+2)^2)+2*y^2/(y^2+(x+2)^2)^2

 

u[y] <> -v[x]

 

-2*(x+2)*y/(y^2+(x+2)^2)^2 <> -y*(2*x+4)/(y^2+(x+2)^2)^2

 

 

 

`Cauchy-Riemann ?`

(2)

 

Also ran the procedure through the AI for so-called code improvement and now it shows what the code stands for

restart;

# Improved and corrected version of the CauchyRiemann procedure :ASKED AI 
CauchyRiemann := proc(expr::algebraic)
    local x, y, u, v, u_x, u_y, v_x, v_y, CR1, CR2;

    # Assign real and imaginary parts of the function
    u := evalc(Re(eval(expr, z = x + I*y)));
    v := evalc(Im(eval(expr, z = x + I*y)));

    # Calculate partial derivatives
    u_x := diff(u, x);
    u_y := diff(u, y);
    v_x := diff(v, x);
    v_y := diff(v, y);

    # Properly format and print function details
    printf("f(z) = %a\n", expr);
    printf("u(x, y) = %a, u_x = %a, u_y = %a\n", u, u_x, u_y);
    printf("v(x, y) = %a, v_x = %a, v_y = %a\n", v, v_x, v_y);

    # Evaluate and print Cauchy-Riemann equations
    CR1 := u_x = v_y;
    CR2 := u_y = -v_x;
    printf("\nCauchy-Riemann Equations:\n");
    printf("u_x = v_y: %a\n", CR1);
    printf("u_y = -v_x: %a\n", CR2);

    # Check both equations
    if CR1 and CR2 then
        printf("The function is analytic (holomorphic) at this point.\n");
        printf("The derivative f'(z) is %a + I*%a\n", u_x, v_y);
    else
        printf("The function does not satisfy the Cauchy-Riemann equations and is not analytic.\n");
    end if;
end proc;

# Test the procedure with a specific function
f := z -> 1/(z + 2);
CauchyRiemann(f(z));

"maple.ini in users"

 

proc (expr::algebraic) local x, y, u, v, u_x, u_y, v_x, v_y, CR1, CR2; u := evalc(Re(eval(expr, z = x+I*y))); v := evalc(Im(eval(expr, z = x+I*y))); u_x := diff(u, x); u_y := diff(u, y); v_x := diff(v, x); v_y := diff(v, y); printf("f(z) = %a
", expr); printf("u(x, y) = %a, u_x = %a, u_y = %a
", u, u_x, u_y); printf("v(x, y) = %a, v_x = %a, v_y = %a
", v, v_x, v_y); CR1 := u_x = v_y; CR2 := u_y = -v_x; printf("
Cauchy-Riemann Equations:
"); printf("u_x = v_y: %a
", CR1); printf("u_y = -v_x: %a
", CR2); if CR1 and CR2 then printf("The function is analytic (holomorphic) at this point.
"); printf("The derivative f'(z) is %a + I*%a
", u_x, v_y) else printf("The function does not satisfy the Cauchy-Riemann equations and is not analytic.
") end if end proc

 

proc (z) options operator, arrow; 1/(z+2) end proc

 

f(z) = 1/(z+2)
u(x, y) = (x+2)/(y^2+(x+2)^2), u_x = 1/(y^2+(x+2)^2)-(x+2)/(y^2+(x+2)^2)^2*(2*x+4), u_y = -2*(x+2)/(y^2+(x+2)^2)^2*y
v(x, y) = -y/(y^2+(x+2)^2), v_x = y/(y^2+(x+2)^2)^2*(2*x+4), v_y = -1/(y^2+(x+2)^2)+2*y^2/(y^2+(x+2)^2)^2

Cauchy-Riemann Equations:
u_x = v_y: 1/(y^2+(x+2)^2)-(x+2)/(y^2+(x+2)^2)^2*(2*x+4) = -1/(y^2+(x+2)^2)+2*y^2/(y^2+(x+2)^2)^2
u_y = -v_x: -2*(x+2)/(y^2+(x+2)^2)^2*y = -y/(y^2+(x+2)^2)^2*(2*x+4)
The function does not satisfy the Cauchy-Riemann equations and is not analytic.

 

NULL

Download CAUCHY_RIEMANN_-FORUM_VRAAG.mw

I want to solve or try to solve this equation 

PDE := diff(G(a, H, phi, PI), a)(aH) + diff(G(a, H, phi, PI), H)(k/a^2 - kappa^2/2*PI^2/a^6) + diff(G(a, H, phi, PI), phi)(PI/a^3) = diff(G(a, H, phi, PI), PI)(a^3*diff(V(phi), phi))

with pdsolve(PDE, G)

and maple answer me the next

Error, (in pdsolve/info) first argument does not have a differentiated function with name G

I nw in maple, maybe I´m make a mistake, but I can't find what

I have  4 worksheets with derived equations. So I export the equations and  possibly some procedures (but they can be handled seperately if needed)  from each worksheet as a .mpl file. 

I want to combine the .mpl files  together without using copy/paste. Then I can open that single file in the VS code editor.
There may be other ways to achieve this so I am open to suggestions.

I chased down a problem to factoring a square that has sqrt in the coefficients. All numbers are real,
The code is inside a procedure in a package. Iso I could do with something robust.

expand((sqrt(A+B)*x+sqrt(7-K)*y)^2)
     2      2            (1/2)          (1/2)        2      2
  A x  + B x  + 2 (A + B)      x (7 - K)      y - K y  + 7 y 

factor(%) 

 

Hello, everyone,

I am new to Maple and I am trying to get use of it.

I tried to plot the following linear systems in different ways. I realized that the Student Linear Algebra is not as flexible as Linear Algebra. My question is the following. Is there any other way to create a plot without defining the implicit plots?

with(Student[LinearAlgebra])

A := Matrix([[1, 1], [12, 16]]); b := Vector([10, 136]); sol := LinearSolve(A, b)

Matrix(2, 2, {(1, 1) = 1, (1, 2) = 1, (2, 1) = 12, (2, 2) = 16})

 

Vector(2, {(1) = 10, (2) = 136})

 

Vector[column](%id = 36893488153382714652)

(1)

LinearSystemPlot({x+y = 10, 12*x+16*y = 136}, axes = normal)

 

restart

 

 

with(Student[LinearAlgebra])

A := Matrix([[2, -1, 1], [0, 1, 3], [0, 0, 1]]); b := Vector([-5, 7, 2]); sol := LinearSolve(A, b)

Matrix(3, 3, {(1, 1) = 2, (1, 2) = -1, (1, 3) = 1, (2, 1) = 0, (2, 2) = 1, (2, 3) = 3, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1})

 

Vector(3, {(1) = -5, (2) = 7, (3) = 2})

 

LinearSolve(Matrix(%id = 36893488151878044716), Vector[column](%id = 36893488151878030628))

(2)

LinearSystemPlot({w = 2, y+3*w = 7, 2*x-y+w = -5}, axes = normal)

 

restart

with(plots); with(LinearAlgebra); A := Matrix([[1, 1], [12, 16]]); b := Vector([10, 136]); sol := LinearSolve(A, b); eq1 := x+y = 10; eq2 := 12*x+16*y = 136; plot1 := implicitplot(eq1, x = -5 .. 10, y = -50 .. 50, color = "red", thickness = 2, labels = ["x", "y"]); plot2 := implicitplot(eq2, x = -5 .. 10, y = -50 .. 50, color = "blue", thickness = 2); display(plot1, plot2, title = "Plot of Linear System", legend = [x+y = 10, 12*x+16*y = 136])

Matrix(2, 2, {(1, 1) = 1, (1, 2) = 1, (2, 1) = 12, (2, 2) = 16})

 

Vector(2, {(1) = 10, (2) = 136})

 

Vector[column](%id = 36893488153330030820)

 

x+y = 10

 

12*x+16*y = 136

 

 

A := Matrix([[2, -1, 1], [0, 1, 3], [0, 0, 1]]); b := Vector([-5, 7, 2]); solution := LinearSolve(A, b); eq1 := 2*x-y+z = -5; eq2 := y+3*z = 7; eq3 := z = 2; plot1 := implicitplot3d(eq1, x = -10 .. 10, y = -10 .. 10, z = -10 .. 10, color = "red", style = surface); plot2 := implicitplot3d(eq2, x = -10 .. 10, y = -10 .. 10, z = -10 .. 10, color = "blue", style = surface); plot3 := implicitplot3d(eq3, x = -10 .. 10, y = -10 .. 10, z = -10 .. 10, color = "green", style = surface); display(plot1, plot2, plot3, title = "3D Plot of Linear System", axes = boxed)

Matrix(3, 3, {(1, 1) = 2, (1, 2) = -1, (1, 3) = 1, (2, 1) = 0, (2, 2) = 1, (2, 3) = 3, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1})

 

Vector(3, {(1) = -5, (2) = 7, (3) = 2})

 

Vector[column](%id = 36893488152610632156)

 

 

NULL


 

Download linear_systems.mw

 

OEIS A034828 and OEIS A000292 (which give the Wiener index for the cycle graph and the path graph respectively) mention that 

the Wiener index of the cycle of length 19 is 855 and 
the Wiener index of the path with 19 edges is 1330

However, 

GraphTheory:-WienerIndex(GraphTheory:-CycleGraph(19));
 = 
                               38

GraphTheory:-WienerIndex(GraphTheory:-PathGraph(20));
 = 
                               38

So what happened here? 

First 36 37 38 39 40 41 42 Page 38 of 43