Download Leontief2.mw

Download odesys.mw

GraphTheory:-Graph will accept an adjacency matrix but not an incidence matrix. Here is a procedure that will convert an incidence matrix into an adjacency matrix.

IncidenceToAdjacency:= proc(J::Matrix)
uses LA= LinearAlgebra;
local
     r:= LA:-RowDimension(J),
     A:= Matrix(r,r, shape= symmetric),
     c:= LA:-ColumnDimension(J),
     j, ones
;
     for j to c do
          ones:= ArrayTools:-SearchArray(J[..,j]);
          A[ones[1], ones[2]]:= 1
     end do;
     A
end proc:

(This procedure only handles the case of an unweighted, undirected graph.)

Using your example incidence Matrix:

M:=Matrix([[1,1,0,1,0,0,0],[0,1,1,0,1,0,0],[1,0,1,0,0,1,1],[0,0,0,1,1,1,0],[0,0,0,0,0,0,1]]):
A:= IncidenceToAdjacency(M);

G:= GraphTheory:-Graph(A):
GraphTheory:-DrawGraph(G);

Here is a complete simulation for the three-door Monty Hall problem (the classic). The code is very straightforward to read. See if you can extend this to the four-door problem.

MontyHall:= proc(N::posint, {switch::truefalse:= false})
# Simulation for the 3-door Monty Hall problem. N is the number
# of trials. Returns the number of trials the player wins.
local
     wins:= 0,
     PrizeDoor, #The door with the prize behind it (1..3).
     MontysDoor,  # The door that Monty reveals (1..3).
     Players1stDoor,  # The door the player picks first.
     Players2ndDoor,  # The door the player picks second.
     Doors:= {1,2,3},
     rand3:= rand(1..3), rand2:= rand(1..2)
;
     to N do
          PrizeDoor:= rand3();
          Players1stDoor:= rand3();
          if Players1stDoor = PrizeDoor then
               MontysDoor:= (Doors minus {PrizeDoor})[rand2()]
          else
               MontysDoor:= (Doors minus {PrizeDoor,Players1stDoor})[]
          end if;
          if switch then
               Players2ndDoor:= (Doors minus {MontysDoor,Players1stDoor})[]
          else
               Players2ndDoor:= Players1stDoor
          end if;
          if Players2ndDoor = PrizeDoor then
               wins:= wins+1
          end if
     end do;
     wins
end proc:     

MontyHall(100000);
                             33282
MontyHall(100000, switch);
                             66915

You asked:

How would i label the axis and change its font and size?

By using plot options labels and labelfont. (Example shown below.)

Also, how would i move the horizontal axis down, so that the whole diagram can be seen?

By using plot options axes= box and view. (Example shown below.)

Finally, there are some odd points between 0 and 1 on the vertical axis. How do you get rid of these?

Those are very interesting and seem to be an integral result of the computation for small values of r and some x, although it may also be an effect of round-off error. I will investigate that further. To eliminate them from the plot, slightly increase the lower bound for the view of the horizontal axis. I used 0.1 as the increase amount in the example show below.

Bonus: Here's how you can generate this plot in well under 1 second by using evalhf. It's at least a factor-of-50 improvement timewise over the old way. I set this up such that it will be easy to modify for other bifurcation diagrams.

restart:
Rloop:= proc(Pts::Matrix, r_min::realcons, r_max::realcons, N::posint, M::posint)
local n, r, x, r_range:= r_max - r_min;
     for n to N do
          r:= r_min+n/N*r_range;
          x:= Pts[n,2];
          to M do  x:= x*exp(r*(1-x))  end do;
          Pts[n,1]:= r;  Pts[n,2]:= x
     end do;
     NULL
end proc:     

(N,M):= 10^~(4,2):  (x_min,x_max):= (0,1):  (r_min,r_max):= (0,4):
bifpoint:= Matrix(N,2, datatype= float[8]):
bifpoint[.., 2]:=
     < RandomTools:-Generate(list(float(range= x_min..x_max, method= uniform), N))[] >:
evalhf(Rloop(bifpoint, r_min, r_max, N, M)):
pitchf:= plots:-pointplot(bifpoint, symbol=point):
plots:-display(
     pitchf,
     axes= box, view= [r_min+0.1..r_max, -1..5],
     labels= ['r','x'], labelfont= [TIMES,BOLDITALIC,16],
     axesfont= [HELVETICA,BOLD,12]
);

You can use the events option to dsolve to detect when the derivative is 0. If you post the equations, I'll help to set it up. Otherwise, you can read about it at ?dsolve,events . (Warning: This help page is quite esoteric.)

While it's true that you can't use solve with the output of dsolve(..., numeric), you can use fsolve.

I believe that your error "Empty script base" is related to the 2D input. I have never seen it before. When I convert your expression to 1D input with Convert tool on the Format menu, it looks very weird to me:

E~:=(k__x,k__y,k__z)~->~&+-Ɣsqrt((1)+(4~cos(3/(2)k__ya__cc)cos((((sqrt(3)))/(2))k__xa__cc))+((4)((cos^(2)(((sqrt(3))/(2))k__xa__cc))^()))6";

The tildes (~) and quotation marks make no sense to me.

Here is how I would type your expression in 1D input:

E:= (k__x, k__y, k__z)->
    &+- gamma*sqrt(1+4*cos(3/2*k__y*a__cc)*cos(sqrt(3)/2*k__x*a__cc) +
               4*cos(sqrt(3)/2*k__x*a__cc)^2);

To substitute the values for gamma and a__cc, do

E:= subs([gamma= 3*1.6e-19, a__cc= 0.142e-9], eval(E));

Names created with parse are always global. So you'll need to make PEA global instead of local for your code to work.

If you can describe what you're actually trying to accomplish, I can probably suggest a better way to do it.

Maple allows the upper bound of a plot range to be infinity:

plots([n^3, (n+1)^2], n= 3..infinity, legend= [n^3, (n+1)^2]);

or

plot(n^3 - (n+1)^2, n= 3..infinity);

Maple can prove some simple inequalities "on its own":

is(n^3 > (n+1)^2) assuming n >= 3;
                                   true

Here is a procedure that might help you. It adds a Vector to an Array iff the Vector does not already appear in the Array. It returns true if the Vector was added and false otherwise.

AddIfNew:= proc(A::{Array,Vector}, V::rtable)
local n:= numelems(A);
     if not ormap(LinearAlgebra:-Equal, A, V) then
          A(n+1):= V;
          return true
     end if;
     false
end proc:

Example of its use:

A:= Array([]):
for k to 26 do
     AddIfNew(A, LinearAlgebra:-RandomVector(2, generator= rand(1..5)))
end do:
numelems(A);
                               17
convert(A,list);

It is extremely difficult to get fsolve to produce all the solutions to a system of equations and to know that they are truly all of the solutions. Often it is impossible. In this case, the equations are solvable exactly, so use solve:

seq(eval(y,E), E in solve({a[1],a[2]}, {x,y}, explicit));

To get just the positives, produce the above list and then filter it with select and is like I showed you before:

Sols:= seq(eval(y,E), E in solve({a[1],a[2]}, {x,y}, explicit)):
select(y-> is(y>0), Sols);

If ex is your expression of square roots, try

expand(combine(ex, symbolic));

f:= unapply(Trace(M), x);

BuildMatrix:= (ex::algebraic)->
     Matrix(10, 9, (i,j)-> eval(ex, {C= i, K= j+1}))
:
BuildMatrix(2*K+2*C-3);

@john125 

The issue of dealing with multiple solutions from fsolve only arises when you are solving a single univariate polynomial. In all other cases fsolve returns at most one solution per invocation.

Let p(x) be the polynomial that you are solving. To eliminate duplicates, enclose the fsolve is set-constructor braces:

{fsolve(p(x), x)};

or

{fsolve}(p(x), x);

To only receive nonnegative solutions, include an interval for solutions in the call:

fsolve(p(x), x, 0..infinity);

To also exclude 0, thus making it just positive solutions, do

fsolve(p(x), x, avoid= {x=0}, 0..infinity);

Either of these latter two ideas can be combined with the braces to eliminate duplicates:

{fsolve(p(x), x, avoid= {x=0}, 0..infinity)};

In your second followup, you ask about solving multiple equations for just one of the variables. That is not possible, although you can easily make it so that you only see the solution for one variable. Let's say you just want to see x. Then do

eval(x, fsolve({a,b}, {x,y}));