Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 172 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

On reflection...

tomleslie 13876
July 01 2019
0 0

would't it be easier to use the geometry pacjkage for this?

As in the attached

> 

with(geometry):
point(A, [-5, 6]):
point(B, [7, -3]):
MakeSquare(sqr, [A, B, diagonal]):
draw(sqr, axes=none);

A

  

B

  

sqr

  

 
> 

detail(sqr);;

GeometryDetail(["name of the object", sqr], ["form of the object", square2d], ["the four vertices of the square", [[-5, 6], [-7/2, -9/2], [7, -3], [11/2, 15/2]]], ["the length of the diagonal", sqrt(225)])

 (1)
> 

 


 

Download bsquare.mw

 

Not sure why you have an error...

tomleslie 13876
July 01 2019
0 1

because the following works.


 

> 

restart;
interface(version);
car_2som_opp := proc (U::list, V::list)  #construction d'un carré connaissant 2 sommets opposés
                     local dist, eqCerU, eqCerV, r, sol, X, Y;
                     dist := proc (M, N)
                                  sqrt(add((M[i]-N[i])^2, i = 1 .. 2))
                             end proc;
                     r := dist(U, V)/sqrt(2);
                     eqCerU := (x-U[1])^2+(y-U[2])^2 = r^2;
                     eqCerV := (x-V[1])^2+(y-V[2])^2 = r^2;
                     sol := solve([eqCerU, eqCerV], [x, y],allsolutions,explicit);  
                     map(allvalues,sol):
                     X := [subs(op(sol[1]), x), subs(op(sol[1]), y)];
                     Y := [subs(op(sol[2]), x), subs(op(sol[2]), y)];
                     plots:-display(plot([U, X, V, Y, U],scaling = constrained, axes = none))
                 end proc:

car_2som_opp([-5,6],[7,-3]);

`Standard Worksheet Interface, Maple 2018.2, Windows 7, November 16 2018 Build ID 1362973`

  

 
> 

 


 

Download asquare.mw

Not completely clear what you want...

tomleslie 13876
July 01 2019
2 0

and you will need at keast one boundary condition (which I have "guessed"), but maybe something in the attached will cover your requirement


 

> 

  restart;
#
# Define ODE using piecewise()
#
  ode:= diff(B[1](t),t)= piecewise
                         ( t<=1000,
                           kaC*(R-B[1](t))-kd1*B[1](t),
                           -kd1*B[1](t)
                         );
#
# Add a "guesswork" boundary condition and
# solve the ODE
#
  sol:=dsolve([ode, B[1](0)=1]);
#
# Plot the ODE solutions for values of
# kd1 fromm 7e-03 to 7e-02
#
  plot( [ seq
          ( eval
            ( rhs(sol),
              [kaC=6e-02, kd1=j, R=1]
            ),
            j=7e-03..7e-02, 1e-03
          )
        ],
        t=0..2000
      );
#
# Plot the ODE solutions for values of
# kaC from 6e-02 to 6e-01
#
  plot( [ seq
          ( eval
            ( rhs(sol),
              [kaC=j, kd1=7e-03, R=1]
            ),
            j=6e-02..6e-01, 1e-02
          )
        ],
        t=0..2000
      );
#
# Plot B[1](t) versus t and kd1, with kaC
# as a parameter (produces mulitple surfaces)
#
  colors:=[red, green, blue, brown, black, yellow]:
  plots:-display
         ( [ seq
             ( plot3d
               ( eval
                 ( rhs(sol),
                   [kaC=6e-02+(j-1)*2e-02, R=1]
                 ),
                 t=0..2000,
                 kd1=7e-03..7e-02,
                 color=colors[j]
               ),
               j=1..6
             )
           ]
         );

ode := diff(B[1](t), t) = piecewise(t <= 1000, kaC*(R-B[1](t))-kd1*B[1](t), -kd1*B[1](t))

  

B[1](t) = piecewise(t < 1000, exp(-kaC*t-kd1*t)*(1-kaC*R/(kaC+kd1))+kaC*R*exp((kaC+kd1)*t)*exp(-kaC*t-kd1*t)/(kaC+kd1), 1000 <= t, kaC*R*exp(1000*kaC+1000*kd1)*exp(-kd1*t-1000*kaC)/(kaC+kd1)+exp(-kd1*t-1000*kaC)*(1-kaC*R/(kaC+kd1)))

  

 

 

 
> 

 


 

Download plotODE.mw

The command...

tomleslie 13876
July 01 2019
1 1

numbers := [random[empirical[0.5, 0.5]](N)]

uses the 'stats' package, which has been 'deprecated' since Maple version 9.5 in 2004. You cannot access these commands using 'with(Statistics)'. You have to use 'with(stats)' , as shown in the attached - which actually "works".

Note that you really shouldn't be using such deprecated packages without a very very good reason. Since all this command actually does is produce a list containing the integers 1,2 at random - there really doesn't seem to be any sensible reason for continuing to use it!

Also

  1. in the final for loop: when you are just adding a finite list of "known" numbers, it is much more efficient to use add() rather than sum()
  2. the procedure Wpath() is defined but never used - deliberate?


 

> 

restart
> 

with(stats)
> 

Wpath := proc (steps, t) local walk, i, N, ww; N := nops(steps); walk[0] := 0; for i from 0 to N-1 do walk[i+1] := walk[i]+steps[i+1]*sqrt(t/N) end do; ww := seq(plot(walk[i], t*i/N .. t*(i+1)/N), i = 0 .. N-1); plots[display]([ww]) end proc
> 

N := 400
> 

numbers := [random[empirical[.5, .5]](N)]
> 

st1 := map(proc (x) options operator, arrow; 2*x-3 end proc, numbers)
> 

list_of_k := [40, 20, 10, 5, 2, 1]
> 

for j to nops(list_of_k) do k := list_of_k[j]; st[k] := [seq((sum(st1[p], p = i*k-k+1 .. i*k))/sqrt(k), i = 1 .. N/k)] end do

40

  

[(1/10)*10^(1/2), -(1/5)*10^(1/2), (2/5)*10^(1/2), -(1/10)*10^(1/2), (2/5)*10^(1/2), -(2/5)*10^(1/2), -(3/10)*10^(1/2), 0, 0, (1/5)*10^(1/2)]

  

20

  

[-(1/5)*5^(1/2), (2/5)*5^(1/2), 0, -(2/5)*5^(1/2), (1/5)*5^(1/2), (3/5)*5^(1/2), -(3/5)*5^(1/2), (2/5)*5^(1/2), -(1/5)*5^(1/2), 5^(1/2), -(4/5)*5^(1/2), 0, 0, -(3/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), (2/5)*5^(1/2), -(2/5)*5^(1/2), (2/5)*5^(1/2), 0]

  

10

  

[-(1/5)*10^(1/2), 0, 0, (2/5)*10^(1/2), 0, 0, -(1/5)*10^(1/2), -(1/5)*10^(1/2), (1/5)*10^(1/2), 0, 0, (3/5)*10^(1/2), -(2/5)*10^(1/2), -(1/5)*10^(1/2), 0, (2/5)*10^(1/2), -(1/5)*10^(1/2), 0, (2/5)*10^(1/2), (3/5)*10^(1/2), -(3/5)*10^(1/2), -(1/5)*10^(1/2), (2/5)*10^(1/2), -(2/5)*10^(1/2), (1/5)*10^(1/2), -(1/5)*10^(1/2), -(3/5)*10^(1/2), 0, -(2/5)*10^(1/2), (1/5)*10^(1/2), (2/5)*10^(1/2), -(1/5)*10^(1/2), 0, (2/5)*10^(1/2), -(2/5)*10^(1/2), 0, (2/5)*10^(1/2), 0, (1/5)*10^(1/2), -(1/5)*10^(1/2)]

  

5

  

[-(3/5)*5^(1/2), (1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (3/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), -(1/5)*5^(1/2), -(1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), (1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), 5^(1/2), (1/5)*5^(1/2), -5^(1/2), -(1/5)*5^(1/2), -(1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), (1/5)*5^(1/2), (3/5)*5^(1/2), (1/5)*5^(1/2), -(3/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (3/5)*5^(1/2), (1/5)*5^(1/2), (3/5)*5^(1/2), (3/5)*5^(1/2), -(1/5)*5^(1/2), -5^(1/2), (1/5)*5^(1/2), -(3/5)*5^(1/2), -(1/5)*5^(1/2), 5^(1/2), -(1/5)*5^(1/2), -(3/5)*5^(1/2), (1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), -(1/5)*5^(1/2), -5^(1/2), -(1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), -(3/5)*5^(1/2), (1/5)*5^(1/2), (1/5)*5^(1/2), 5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), -(3/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), 5^(1/2), -5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), (1/5)*5^(1/2), (3/5)*5^(1/2), (3/5)*5^(1/2), -(3/5)*5^(1/2), -(3/5)*5^(1/2), 5^(1/2), -(3/5)*5^(1/2), (1/5)*5^(1/2)]

  

2

  

[-2^(1/2), 0, 0, 0, 0, 2^(1/2), -2^(1/2), 0, 0, 0, 2^(1/2), -2^(1/2), 0, -2^(1/2), 2^(1/2), 0, 2^(1/2), 0, 0, 2^(1/2), 2^(1/2), -2^(1/2), -2^(1/2), 2^(1/2), 0, 0, 0, 0, 2^(1/2), -2^(1/2), -2^(1/2), 0, 2^(1/2), -2^(1/2), 0, 0, 0, 0, -2^(1/2), 0, 2^(1/2), -2^(1/2), 2^(1/2), 0, 0, 0, 2^(1/2), 0, 0, -2^(1/2), 2^(1/2), 0, -2^(1/2), 2^(1/2), -2^(1/2), 2^(1/2), 0, 0, 2^(1/2), 2^(1/2), 0, 2^(1/2), -2^(1/2), -2^(1/2), -2^(1/2), 0, -2^(1/2), 0, 0, 0, -2^(1/2), 0, 0, 2^(1/2), 0, 0, 2^(1/2), -2^(1/2), 2^(1/2), 2^(1/2), 0, 2^(1/2), -2^(1/2), -2^(1/2), 0, 2^(1/2), 0, -2^(1/2), 0, 0, 2^(1/2), 0, 2^(1/2), 0, 0, 2^(1/2), 0, 2^(1/2), 2^(1/2), 0, 0, -2^(1/2), 0, -2^(1/2), -2^(1/2), 0, 2^(1/2), -2^(1/2), -2^(1/2), 0, -2^(1/2), 2^(1/2), 0, 2^(1/2), 2^(1/2), -2^(1/2), 2^(1/2), -2^(1/2), -2^(1/2), 0, 0, 0, 0, 2^(1/2), 0, 0, -2^(1/2), 2^(1/2), -2^(1/2), 0, -2^(1/2), -2^(1/2), -2^(1/2), 2^(1/2), -2^(1/2), -2^(1/2), 0, 2^(1/2), 2^(1/2), -2^(1/2), 0, -2^(1/2), 0, -2^(1/2), 0, 0, 0, 0, 2^(1/2), 0, 2^(1/2), 2^(1/2), 2^(1/2), 0, -2^(1/2), 0, 0, 0, -2^(1/2), 0, 0, 0, 0, 0, 0, -2^(1/2), 2^(1/2), 0, 2^(1/2), 2^(1/2), -2^(1/2), -2^(1/2), 0, -2^(1/2), 2^(1/2), -2^(1/2), 0, 2^(1/2), -2^(1/2), 2^(1/2), 0, 2^(1/2), -2^(1/2), 2^(1/2), 2^(1/2), 0, 2^(1/2), 0, -2^(1/2), 0, -2^(1/2), -2^(1/2), 2^(1/2), 2^(1/2), 2^(1/2), -2^(1/2), 0, -2^(1/2), 2^(1/2), 0]

  

1

  

[-1, -1, -1, 1, -1, 1, -1, 1, -1, 1, 1, 1, -1, -1, 1, -1, 1, -1, 1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, 1, -1, 1, 1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, -1, 1, -1, 1, 1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, -1, 1, -1, 1, -1, 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, -1, 1, 1, -1, 1, 1, -1, -1, 1, 1, 1, -1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1, 1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, -1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, -1, 1, 1, -1, 1, 1, -1, -1, 1, 1, 1, 1, -1, 1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 1, -1, 1, -1, -1, -1, 1, 1, -1, 1, 1, 1, -1, 1, 1, 1, -1, -1, 1, 1, 1, -1, 1, 1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, -1, -1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, -1, -1, -1, 1, 1, -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, 1, -1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, -1, 1, 1, 1, -1, 1, 1, -1, -1, 1, 1, 1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, -1]

 (1)
> 

``


 

Download walk.mw

 

 

What exactly is wrong?...

tomleslie 13876
July 01 2019
1 3

Somehow, the "code" in your file Matrix.mw appears to be "Maple Output" which is odd

If I convert everything to 1-D Maple Input, then everything seems to execute correctly - see the attached

> 

  restart;
  interface(rtablesize = 10);
> 

  with(LinearAlgebra);
> 

  A := 8:
  B := 5:
  q := 0.4:
  p := 0.2:
  r := 1 - p - q:
  dimP := A + B + 1:
> 

  P := Matrix(dimP, dimP, [0 $ dimP*dimP]):
  P[1, 1] := 1:
  P[1, 2] := 0:
  P[dimP, dimP] := 1:
  P[dimP, dimP - 1] := 0:
  for i from 2 to dimP - 1 do
      P[i, i - 1] := q;
      P[i, i] := r;
      P[i, i + 1] := p;
  end do:
> 

  p0 := Matrix(dimP, 1, [0 $ dimP]):
  p0[A + 1, 1] := 1:
  pV[0] := p0:
  PT := Transpose(P):
> 

  for n to 200 do
      pV[n] := PT . (pV[n - 1]):
  end do:
> 

  map(x -> evalf(x, 3), Transpose(pV[5]));

_rtable[18446744074332538990]

 (1)
> 

 

NULL

Download mat2.mw

Difficult to tell...

tomleslie 13876
June 29 2019
0 0

what you are getting wrong unless you upload the relevant worksheet using the big green up-arrow in the toolbar.

The attached will produce the plot you want

> 

plot(x^2-9, x=-10..10);

 
> 

 

Download simplePlot.mw

One way...

tomleslie 13876
June 27 2019
1 8

is to get rid of the integral term, by differentiating the pde with respect to 't'.

The drawback with this approach is hte at it increases the degree of the PDE (with respect to 't') so an additional initial condition is needed. I tried a few possibilities for D[2](u)(x,0). These didn't make a huge difference to the overall "form" of the solution, whihc always gets *seriously* huge for t>~1.
You might want to experiment with the following, and in particular with the "extra" initial condition D[2](u)(x,0=1, which I inserted. (Changing the rhs of this initial condition doesn't make a huge difference to the overall form!)

> 

interface(rtablesize=10):
f:=(x^3+t^2*x^2-t^2*x+4*t*x-2*t+1)*exp(x*t)-x+1;
pde:=diff(u(x,t),t)+diff(u(x,t),x,x)+u(x,t)+int(u(x,tau),tau=0..t)=f;
pde:=diff(pde,t);
bounds:= u(0,t)=0, u(1,t)=0, u(x,0)=x*(x-1), D[2](u)(x,0)=1;
psol:=pdsolve(pde, [bounds], numeric);
psol:-plot3d( u(x,t), x=0..1, t=0..0.1);

(t^2*x^2-t^2*x+x^3+4*t*x-2*t+1)*exp(t*x)-x+1

  

diff(u(x, t), t)+diff(diff(u(x, t), x), x)+u(x, t)+int(u(x, tau), tau = 0 .. t) = (t^2*x^2-t^2*x+x^3+4*t*x-2*t+1)*exp(t*x)-x+1

  

diff(diff(u(x, t), t), t)+diff(diff(diff(u(x, t), t), x), x)+diff(u(x, t), t)+u(x, t) = (2*t*x^2-2*t*x+4*x-2)*exp(t*x)+(t^2*x^2-t^2*x+x^3+4*t*x-2*t+1)*x*exp(t*x)

  

u(0, t) = 0, u(1, t) = 0, u(x, 0) = x*(x-1), (D[2](u))(x, 0) = 1

  

_m906819456

  

 
> 

 

Download pdeSol2.mw

Something like...

tomleslie 13876
June 27 2019
0 0

the attached -maybe?

> 

  restart;
  interface(rtablesize=12);
  M2:=Matrix(2, 12, (i,j)->`if`( i=1,
                                 `if`(type(j, odd), U[(j+1)/2],0),
                                 `if`(type(j, even), U[j/2], 0)
                               )
            );
  Matrix( [ [seq( D[1](M2[1,j])(x,y) + 0*M2[2,j],          j =1..12)],
            [seq( 0*(M2[1,j])(x,y)   + D[2](M2[2,j])(x,y), j =1..12)],
            [seq( D[2](M2[1,j])(x,y) + D[1](M2[2,j])(x,y), j =1..12)]
          ]
        ):
  convert~([%], diff);

[10, 10]

  

Matrix(2, 12, {(1, 1) = U[1], (1, 2) = 0, (1, 3) = U[2], (1, 4) = 0, (1, 5) = U[3], (1, 6) = 0, (1, 7) = U[4], (1, 8) = 0, (1, 9) = U[5], (1, 10) = 0, (1, 11) = U[6], (1, 12) = 0, (2, 1) = 0, (2, 2) = U[1], (2, 3) = 0, (2, 4) = U[2], (2, 5) = 0, (2, 6) = U[3], (2, 7) = 0, (2, 8) = U[4], (2, 9) = 0, (2, 10) = U[5], (2, 11) = 0, (2, 12) = U[6]})

  

[Matrix(%id = 18446744074423018070)]

 (1)
> 

 


 

Download doDiffs.mw

Try the attached...

tomleslie 13876
June 27 2019
1 0

which appears to give the same answer

> 

restart
> 

N := 2; A := -N; B := N
> 

q := .3; p := .5; sa := .9; sb := .1; r := 1-p-q
> 

dimP := 2*N+1
> 

P := Matrix(dimP, dimP)
> 

P[1, 1] := sa; P[1, 2] := 1-sa; P[dimP, dimP] := sb; P[dimP, dimP-1] := 1-sb
> 

for i from 2 to dimP-1 do P[i, i-1] := q; P[i, i] := r; P[i, i+1] := p end do
> 

P

Matrix(%id = 18446744074421179446)

 (1)
> 

# change this part code to the modernversion
 with(linalg)
> 

J := diag(`$`(1, dimP))

array( 1 .. 5, 1 .. 5, [( 5, 5 ) = (1), ( 3, 3 ) = (1), ( 4, 4 ) = (1), ( 2, 2 ) = (1), ( 1, 1 ) = (1) ] )

 (2)
> 

d := matrix(dimP, 1, [`$`(1, dimP)])

array( 1 .. 5, 1 .. 1, [( 3, 1 ) = (1), ( 4, 1 ) = (1), ( 5, 1 ) = (1), ( 1, 1 ) = (1), ( 2, 1 ) = (1) ] )

 (3)
> 

b := matrix(dimP+1, 1, [`$`(0, dimP), 1])

array( 1 .. 6, 1 .. 1, [( 3, 1 ) = (0), ( 4, 1 ) = (0), ( 5, 1 ) = (0), ( 1, 1 ) = (0), ( 6, 1 ) = (1), ( 2, 1 ) = (0) ] )

 (4)
> 

augment(P-J, d); A := transpose(augment(P-J, d))

Matrix(5, 6, {(1, 1) = -.1, (1, 2) = .1, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 1, (2, 1) = .3, (2, 2) = -.8, (2, 3) = .5, (2, 4) = 0, (2, 5) = 0, (2, 6) = 1, (3, 1) = 0, (3, 2) = .3, (3, 3) = -.8, (3, 4) = .5, (3, 5) = 0, (3, 6) = 1, (4, 1) = 0, (4, 2) = 0, (4, 3) = .3, (4, 4) = -.8, (4, 5) = .5, (4, 6) = 1, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = .9, (5, 5) = -.9, (5, 6) = 1})

  

array( 1 .. 6, 1 .. 5, [( 5, 4 ) = (.5), ( 4, 3 ) = (.5), ( 5, 5 ) = (-.9), ( 4, 5 ) = (.9), ( 6, 5 ) = (1), ( 1, 4 ) = (0), ( 3, 3 ) = (-.8), ( 4, 2 ) = (0), ( 3, 1 ) = (0), ( 2, 4 ) = (0), ( 4, 4 ) = (-.8), ( 5, 2 ) = (0), ( 1, 3 ) = (0), ( 1, 5 ) = (0), ( 3, 4 ) = (.3), ( 6, 4 ) = (1), ( 5, 3 ) = (0), ( 2, 2 ) = (-.8), ( 4, 1 ) = (0), ( 1, 2 ) = (.3), ( 5, 1 ) = (0), ( 3, 2 ) = (.5), ( 2, 5 ) = (0), ( 6, 3 ) = (1), ( 1, 1 ) = (-.1), ( 2, 3 ) = (.3), ( 3, 5 ) = (0), ( 6, 2 ) = (1), ( 6, 1 ) = (1), ( 2, 1 ) = (.1) ] )

 (5)
> 

linsolve(A, b)

array( 1 .. 5, 1 .. 1, [( 3, 1 ) = (.1668726823), ( 4, 1 ) = (.2781211372), ( 5, 1 ) = (.1545117429), ( 1, 1 ) = (.3003708282), ( 2, 1 ) = (.1001236094) ] )

 (6)
> 

#
# Unload the linalg package and load the
# LinearAlgebra packages
#
  unwith(linalg):
  with(LinearAlgebra):
#
# Same calculation using LinearAlgebra()
# commands
#
  J:= DiagonalMatrix([1 $ dimP]):
  d:= Matrix(dimP, 1, [1 $ dimP]):
  b:= Matrix(dimP + 1, 1, [0 $ dimP, 1]):
  A:= Transpose(<P-J|d>):
  LinearSolve(A,b);
 

Vector(5, {(1) = .30037082818294225, (2) = .10012360939431399, (3) = .16687268232385657, (4) = .27812113720642756, (5) = .15451174289245975})

 (7)
> 

 


 

Download modernize.mw

some complex number problem...

tomleslie 13876
June 27 2019
1 2
  1. As defined 'A' is complex
  2. So everything subsequently based on 'A' will also be complex (??)
  3. Doesn't seem to cause an issue for 'H' and 'B' -maybe because imaginary part are very small???
  4. Can't fieldplot 'E', but can fieldplot 'Re(E)'. It is noticeable that for 'E', imaginary parts are "bigger" than real parts

See the attached

> 

  restart;
> 

  with(plots):
  with(LinearAlgebra):
  with(VectorCalculus):
  SetCoordinates('cartesian'[x,y,z]):
> 

  mu_0:=4*Pi*10^(-7): I_0:=2400: f:=2500: omega:=2*Pi*f:
  c:=3*10^8: k:=omega/c: l:=3*10^(-2): epsilon_0:=1/(mu_0*c^2):
> 

  J:=Vector[column]([ 0 ,
                      0 ,
                      I_0/s ]);

Vector(3, {(1) = 0, (2) = 0, (3) = 2400/s})

 (1)
> 

# IMPORTANT: R is constant for the calculation of A
#
  R:=sqrt(x^2+y^2+z^2);
#
# Note 'A' will be complex!
#
  A := VectorField(mu_0/(4*Pi)*int(J*exp(-I*k*R)*s/R, p = -l/2 .. l/2));

R := sqrt(x^2+y^2+z^2)

  

Vector(3, {(1) = 0, (2) = 0, (3) = (9/1250000)*exp(-((1/60000)*I)*Pi*sqrt(x^2+y^2+z^2))/sqrt(x^2+y^2+z^2)})

 (2)
> 

  B:=Curl(A):
#
# B ought(?) to be complex, because 'A' is - yet
# it plots??? Is fieldplot() perhaps ignoring
# "small" imaginary parts??
#
# No difference between plotting 'B' and 'Re(B)'
#
  B_plot:=fieldplot3d(B,     x=-1..1,y=-1..1,z=-1..1,fieldstrength=log,arrows=SLIM);
  B_plot:=fieldplot3d(Re(B), x=-1..1,y=-1..1,z=-1..1,fieldstrength=log,arrows=SLIM)

 

 
> 

#
# H is complex, because B is. Maybe fieldplot ia
# still ignoring "small" imaginary parts??
#
# No difference between plotting 'H' and 'Re(H)'
#
  H:=(1/mu_0)*B:
  H_plot:=fieldplot3d(H, x=-1..1,y=-1..1,z=-1..1,fieldstrength=log,arrows=SLIM);
  H_plot:=fieldplot3d(Re(H), x=-1..1,y=-1..1,z=-1..1,fieldstrength=log,arrows=SLIM);

 

 
> 

#
# E is complex, because 'H' is. Now there
# is a significant difference between plotting
# 'E' and 'Re(E)'
#
  E:=1/(omega*epsilon_0*I)*Curl(H):
  E_plot:=fieldplot3d(E, x=1..500,y=1..500,z=1..500,fieldstrength=log,arrows=SLIM);
  E_plot:=fieldplot3d(Re(E), x=1..500,y=1..500,z=1..500,fieldstrength=log,arrows=SLIM);

 

 
> 

#
# Check relative magnitudes of real and imaginary
# parts in B, H, and E
#
  evalf(eval(B, [x=1,y=1,z=1]));
  evalf(eval(H, [x=1,y=1,z=1]));
  evalf(eval(E, [x=1,y=1,z=1]));

Vector(3, {(1) = -0.1385640652e-5+0.5000000000e-18*I, (2) = 0.1385640652e-5-0.5000000000e-18*I, (3) = 0.})

  

Vector(3, {(1) = -1.102657795+0.3978873575e-12*I, (2) = 1.102657795-0.3978873575e-12*I, (3) = 0.})

  

Vector(3, {(1) = -0.8908680955e-6-7939136.102*I, (2) = -0.8908680955e-6-7939136.102*I, (3) = -0.3947841759e-5-0.4353118458e-1*I})

 (3)
> 

 


 

Download fields.mw

If...

tomleslie 13876
June 27 2019
1 3

you check the help at ?solve/details you will find that the solve() command will ignore assumptions which cannot be written as polynomial inequalities. Hence you cannot use 'assumptions' to obtain only real solutions.

For real solutions, the recommendation is to use the RealDomain() package. Alternatively you can 'filter' the roots to select only the 'real' ones.

Both approaches are shown in the attached
 

> 

restart;
with(RealDomain):
fprime := x-> x^6-(3/2)*x^5+2*x^4+(5/2)*x^3-7*x^2+2:
f := unapply(simplify(int(fprime(x), x)), x):
g := unapply(expand(f(x^2+2*x)), x):
sol := solve(And(diff(g(x),x)=0,diff(g(x),x,x)<>0),x);
evalf(sol);

-1, -1+(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 4))^(1/2), -1-(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 4))^(1/2), -1+(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 2))^(1/2), -1-(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 2))^(1/2), -1+(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 1))^(1/2), -1-(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 1))^(1/2)

  

-1., -.3051050611, -1.694894939, .497747328, -2.497747328, .283407770, -2.283407770

 (1)
> 

restart;
fprime := x-> x^6-(3/2)*x^5+2*x^4+(5/2)*x^3-7*x^2+2:
f := unapply(simplify(int(fprime(x), x)), x):
g := unapply(expand(f(x^2+2*x)), x):
sol := solve(And(diff(g(x),x)=0,diff(g(x),x,x)<>0),x):
realSol:=seq( `if`( type(evalf(j), realcons), j, NULL), j in sol);
evalf(realSol);

-1, -1+(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 1))^(1/2), -1-(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 1))^(1/2), -1+(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 2))^(1/2), -1-(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 2))^(1/2), -1+(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 4))^(1/2), -1-(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 4))^(1/2)

  

-1., .283407770, -2.283407770, .497747328, -2.497747328, -.3051050611, -1.694894939

 (2)
> 

 


 

Download realSols.mw

 

Depends...

tomleslie 13876
June 26 2019
1 1

just how "fussy" you want to get about the positioning of the labels, adn exactly how they move

The attached is pretty basic, but works

Download rollCube.mw

Like this?...

tomleslie 13876
June 26 2019
2 1
> 

  with(GraphTheory):
> 

  truss:=Graph( { {1, 2}, {1, 3}, {2, 3},
                  {2, 4}, {3 ,4}, {3, 5},
                  {4, 5}, {4, 6}, {5, 6}
                }
              ):
  SetVertexPositions( truss,
                      [ [0, 1], [0, 0], [1, 1],
                        [1, 0], [2, 1], [2, 0]
                      ]
                    ):

  DrawGraph(truss);

 
> 

 


Download truss.mw

You can adjust vertexColors, edgeColors, etc, etc - if yu really want to

Complex values...

tomleslie 13876
June 25 2019
0 1

If I replace the plot() command in your code with

seq( evalf(eval(c1,[gamma2=0.02, gamma1=j])), j=0.02..0.1,0.02)

just to test what values might be obtained, then I get

-0.05316538521 - 1.597731966*I,
-0.05642551576 - 1.600665914*I,
-0.05990502879 - 1.603524111*I,
-0.06361988171 - 1.606284242*I,
-0.06758719718 - 1.608919583*I

In other words every value is complex. Is this desired or expected?

If it is "expected", how were you planning to plot it?

My $0.02...

tomleslie 13876
June 25 2019
0 0

Fairly obviously, for this example, if you evaluate 'sol' under the condition 'f(x)=0', you are going to get divide-by-zero errors.

There are a few wayst to get around this problem - one of them given by vv - and another as shown in the attached

> 

restart;
> 

depvars:= (pde::{algebraic, `=`(algebraic), {set,list}(algebraic, `=`(algebraic))})->
   indets(
      indets(convert(pde, diff), specfunc(diff)),
      And(typefunc(name, name), Not(typefunc({mathfunc, identical(diff)})))
   );

is_solution_trivial:= (pde, sol::{`=`, set(`=`),`function`})->
   evalb(eval(`if`(sol::set, sol, {sol}), depvars(pde)=~ 0) = {0=0});

proc (pde::{algebraic, `=`(algebraic), ({list, set})(algebraic, `=`(algebraic))}) options operator, arrow; indets(indets(convert(pde, diff), specfunc(diff)), And(typefunc(name, name), Not(typefunc({mathfunc, identical(diff)})))) end proc

  

proc (pde, sol::{`=`, function, set(`=`)}) options operator, arrow; evalb(eval(`if`(sol::set, sol, {sol}), `~`[:-`=`](depvars(pde), ` $`, 0)) = {0 = 0}) end proc

 (1)
> 

pde := diff(w(x,y),x)-( diff(f(x),x)*y^2 -f(x)*g(x)*y+ g(x))*diff(w(x,y),y) = 0;

diff(w(x, y), x)-((diff(f(x), x))*y^2-f(x)*g(x)*y+g(x))*(diff(w(x, y), y)) = 0

 (2)
> 

sol:=pdsolve(pde,w(x,y));

w(x, y) = _F1((y*f(x)*(Int((diff(f(x), x))*exp(Int(f(x)*g(x), x))/f(x)^2, x))-f(x)*exp(Int((f(x)^2*g(x)-2*(diff(f(x), x)))/f(x), x))-(Int((diff(f(x), x))*exp(Int(f(x)*g(x), x))/f(x)^2, x)))/(y*f(x)-1))

 (3)
> 

is_solution_trivial(pde,sol)

 

Error, (in eval/Int) numeric exception: division by zero

  
> 

#
# Check output of depvars() - returns the names
# of 'dependent' functions which have derivatives
# in the pde (or pde system)
#
# Does not return the names of other 'dependent'
# functions (eg g(x) in this example) which do
# not have derivatives in the pde
#
  depvars(pde);

{f(x), w(x, y)}

 (4)
> 

#
# The is_solution_trivial() function
#
#  1. sets up the equations depvars(pde)=~0, So,
#     in this example, one obtains
#  2. evaluates its second argument (in this
#     case 'sol') under the condition depvars(pde)=~0
#     ie {f(x) = 0, w(x, y) = 0}
#  3. Thus (for this example), the whole process is
#     equivalent to
#
#     eval(sol, {f(x) = 0, w(x, y) = 0} )
#
#     Since the term 1/f(x) occurs multiple times
#     in sol, this becomes 1/0 and the division-
#     by-zero error is inevitable.
#
# So the following is also guaranteed to give the
# same error
#
  eval(sol, {f(x) = 0, w(x, y) = 0} );

Error, (in eval/Int) numeric exception: division by zero

  
> 

#
# One could circumvent this issue using something
# like
#
  if   not has(sol, 1/~depvars(pde))
  then is_solution_trivial(pde,sol);
  fi;
> 

 

Download div0.mw

First 99 100 101 102 103 104 105 Last Page 101 of 207