Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 168 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

You are complicating a simple problem...

tomleslie 13876
June 07 2020
3 4

and getting bogged down by using clever(?) Maple input methods whihc you appear not to be able to understand or interpret.

So strip your problem down to the basics, and ask yourself what is wrong with the attached (apart from maybe not looking pretty!)

> 

  restart;
  eqs:= -3 <= 6*x - 1,
         6*x - 1 < 3;
  solve([eqs]);

0 <= 6*x+2, 6*x < 4

  

{-1/3 <= x, x < 2/3}

 (1)
> 

 

Download solineq.mw

 

notepad++ on Windows7...

tomleslie 13876
June 06 2020
2 3
  1. For "small" projects, I do pretty much everything in the GUI (utilising 1D text, in worksheet mode)
  2. For "bigger" projects, I write ".mpl" files in Notepad++. The latter comes with "customisation" files for several programming languages. Although the capabilities are pretty basic: not much beyond, syntax highlighting, parenthesis matching and auto-indentation. I find that it is about all I really need. I think I obtained the Maple customisation file for Notepad++ from this site
  3. Notepad++ is an editor, so (AFAIK) there is no way to set up a two-way communication process between it and Maple.
  4. My debugging process is pretty unscientific, basically
    1. Read error messages carefully - they can be informative, then
    2. Read offending code carefully
    3. Insert a few "print" statements in/around the suspect code
    4. If I still have a problem, then I fire up the debugger. I think the best introduction to the debugger is in chapter 16 of the Maple Programming guide - either enter ?programming at the Maple command prompt or download the pdf version from Maplesoft. The user interface to this tool isn't the friendliest I have ever seen, but it does everything you would expect, eg breakpoints, variable "watch", step in/out/over etc etc - see the help at ?debugger for details

Use pdsolve(..., numeric)...

tomleslie 13876
June 02 2020
2 1

You have a partial differetnial equation, so you should probably be using pdsolve.

Using pdsolve() produces the same 'analytic' solution as using dsolve(), I'm somewhat surprised the latter wrked at all!. However the analytic solution contains an infinite sum and an unevaluated integral, so obtaining a 'numeric' solution is probably a better bet. See the attached

> 

  restart;
  eq1 := diff(c(x, t), t) - 1.5*diff(c(x, t), x, x) + 0.05/(x + 0.001) = 0;
  ics := c(x, 0) = 0, c(0, t) = 10, D[1](c)(20, t) = 0;
  #sol := rhs(dsolve({eq1, ics}, c(x, t)));
  sol := pdsolve(eq1, {ics}, numeric);
  sol:-plot3d( c(x,t), x=0..20, t=0..100, labels=[x,t,c]);

diff(c(x, t), t)-1.5*(diff(diff(c(x, t), x), x))+0.5e-1/(x+0.1e-2) = 0

  

c(x, 0) = 0, c(0, t) = 10, (D[1](c))(20, t) = 0

  

_m714060160

  

 
> 

 

> 

 

Download pdsol.mw

Rule is - "last one loaded"...

tomleslie 13876
June 01 2020
1 0

As you load packages using the with() command, a package export which has the same name as a "currently operative name" will "overload" the "currently operative name"

Notice that the VectorCalculus package "overloads" several 'basic' operators, since the command with(VectorCalculus) returns as exports (my emphasis);

[`&x`, `*`, `+`, `-`, `.`, `<,>`, `<|>`, About, AddCoordinates, ArcLength, BasisFormat, Binormal, ConvertVector, CrossProduct, Curl, Curvature, D, Del, DirectionalDiff, Divergence, DotProduct, Flux, GetCoordinateParameters, GetCoordinates, GetNames, GetPVDescription, GetRootPoint, GetSpace, Gradient, Hessian, IsPositionVector, IsRootedVector, IsVectorField, Jacobian, Laplacian, LineInt, MapToBasis, VectorCalculus[Nabla], Norm, Normalize, PathInt, PlotPositionVector, PlotVector, PositionVector, PrincipalNormal, RadiusOfCurvature, RootedVector, ScalarPotential, SetCoordinateParameters, SetCoordinates, SpaceCurve, SurfaceInt, TNBFrame, TangentLine, TangentPlane, TangentVector, Torsion, Vector, VectorField, VectorPotential, VectorSpace, Wronskian, diff, eval, evalVF, int, limit, series]

If you examine what the `.` command from the VectorCalculus() actually does, you can use

showstat( VectorCalculus:-`.`);

whose first few lines are

VectorCalculus:-`.` := proc()
local s, i, const, a, body, DiffOp;
   1   s := [_passed];
   2   DiffOp := {identical('VectorCalculus:-Del'), identical('VectorCalculus
         :-Nabla')};
   3   if not type({_passed},('set')({op(DiffOp), 'Matrix', 'scalar', ':-
         Vector'})) then
   4       return :-`.`(_passed)
       end if;

........a lot more stuff

and notice that after checking the passed arguments, line 4 of the above code returns the "top-level" command `.`  ie :-`.` in prefix form).

This is quite normal, and you can see a similar process for other overloaded operators by examining the code for eg,

showstat( VectorCalculus:-`*`);
showstat( VectorCalculus:-`+`);

#################################################################

If you wish to continue your normal practice of prefixing commands in order to indicate which package you are using, then you can continue to do so - although I think(?) you will have to use prefix rather than infix notation.

Suggest you read the attached (quite carefully)

> 

  restart;
  alias( VC=VectorCalculus):
  f:=y*sin(x)+y^2;
#
# `.` from the vector calculus package
#
  VC:-`.`( VC:-Nabla, VC:-Nabla(f,[x,y]) );
  VC:-DotProduct(VC:-Nabla, VC:-Nabla(f,[x,y]) );
#
# Using top-level '.' in infix and prefix form
#
  VC:-Nabla . VC:-Nabla(f,[x,y]);
  `.`( VC:-Nabla, VC:-Nabla(f,[x,y]) );

f := y*sin(x)+y^2

  

2-y*sin(x)

  

2-y*sin(x)

  

VectorCalculus:-Nabla.(Vector(2, {(1) = y*cos(x), (2) = sin(x)+2*y}))

  

VectorCalculus:-Nabla.(Vector(2, {(1) = y*cos(x), (2) = sin(x)+2*y}))

 (1)
> 

  restart;
  with(VectorCalculus):
  f:=y*sin(x)+y^2;
#
# Infix method
#
  Nabla . Nabla(f,[x,y]);
#
# Prefix method
#
`.`(Nabla, Nabla(f,[x,y]));
#
# Accessing 'top-level' command '.'
#
  :-`.`(Nabla, Nabla(f,[x,y]));

f := y*sin(x)+y^2

  

2-y*sin(x)

  

2-y*sin(x)

  

VectorCalculus:-Nabla.(Vector(2, {(1) = y*cos(x), (2) = sin(x)+2*y}))

 (2)
> 

 


 

Download overload.mw

 

I'm with mmcdara on this one...

tomleslie 13876
May 30 2020
0 0

cos (s)he has posted a correct answer.

Unfortunatley for you it is unliklet that anyoen is going to download a package (Gym-pakke ) to do simple linear fits, becuase

  1. it's in Danish - for me English is fine, French is OK and German, maybe - but Danish? Sorry it isn't going to happen
  2. who needs some add-on packag to do a simple linear fit? This is trivial to do with standard Maple commands as mmcdara has already shown

See also the attached

> 

  restart:
  with(plots):
#
# The data
#
  X:=<0,    10,  20,  30,  40,  50,  60,  70>:
  Y:=<690, 716, 742, 770, 798, 825, 850, 877>:
#
# The linear fit
#
  f:=Statistics:-LinearFit(a+b*x, X, Y, x);
#
# Plots of the data and the 'fit'
#
  plots:-display( [ pointplot(X, Y, symbol=solidcircle, symbolsize=20, color=red),
                    plot(f, x=X[1]..X[-1])
                  ]
                );

HFloat(689.4999999999998)+HFloat(2.6857142857142895)*x

  

 
> 

 

> 

 

> 

 

Download linFit.mw

Main problem...

tomleslie 13876
May 29 2020
0 0

was the use of the parameter 'a', for which there is no numerical value - so I highlighted its use in the attached and (arbitrarliy) gave it the value 1.

I fixed a few other typos

> 

restart
> 

with(plots)
> 

inf := 50
> 

lambda := L(t)+k*A(t)
> 

odes := diff(S(t), t) = Lambda-(`&beta;__1`+`&beta;__2`)*S(t)*lambda-mu*S(t)+r*R(t), diff(E(t), t) = `&beta;__1`*S(t)*lambda-(alpha+mu+d__1+(1-p)*epsilon+p*eta)*E(t), diff(L(t), t) = p*eta*E(t)+e*A(t)-(`&sigma;__2`+mu+d__3+y__2)*L(t), diff(A(t), t) = (1-p)*epsilon*E(t)-(`&sigma;__1`+mu+d__2+y__1+e)*A(t), diff(R(t), t) = y__2*L(t)+y__1*A(t)-(mu+r)*R(t), diff(V(t), t) = `&beta;__2`*S(t)*lambda+alpha*E(t)+`&sigma;__1`*A(t)+`&sigma;__2`*L(t)-a*V(t)
> 

ICs := S(0) = S__0, E(0) = E__0, L(0) = L__0, A(0) = A__0, R(0) = R__0, V(0) = V__0
> 

parameters := [S__0 = 1000000, E__0 = 1, L__0 = 0, A__0 = 0, R__0 = 0, V__0 = 0, Lambda = .22, mu = .22, `&beta;__1` = 0.3e-1, `&beta;__2` = 0.3e-1, r = .5, alpha = .6, p = .1, k = 0.4e-1, y__1 = .8, y__2 = .8, d__1 = 0.2e-1, d__2 = 0.6e-1, d__3 = .2, `&sigma;__1` = 0.2e-1, `&sigma;__2` = 0.2e-1, epsilon = 0.2e-1, eta = 0.2e-1, e = 0.1e-1, a = 1]
> 

S1 := dsolve(eval([odes, ICs], parameters), numeric)
> 

P1 := odeplot(S1, [[t, S(t)]], 0 .. inf, linestyle = 1, color = blue, labels = ["time", "Population"], axes = "boxed"); P2 := odeplot(S1, [[t, E(t)]], 0 .. inf, linestyle = 2, color = brown, labels = ["time", "Population"], axes = "boxed"); P3 := odeplot(S1, [[t, L(t)]], 0 .. inf, linestyle = 3, color = red, labels = ["time", "Population"], axes = "boxed"); P4 := odeplot(S1, [[t, A(t)]], 0 .. inf, linestyle = 4, color = green, labels = ["time", "Population"], axes = "boxed"); P5 := odeplot(S1, [[t, R(t)]], 0 .. inf, linestyle = 5, color = purple, labels = ["time", "Population"], axes = "boxed")

 

 

 

 

 
> 

``


 

Download anotherCOVID.mw

Have I seen this before?...

tomleslie 13876
May 27 2020
0 0

Or at least something very like it?

The attached performs the calculations and animations - alhtough combining the pre-launch and post-launch animations is a bit uninformative becuase of the large discrepancy in both space and time in the two regions

> 

  restart;
  with(plots):
  params:= [alpha=Pi/6, l=40, V__B=4.5, d=3, g=9.81]:
#
# Pre-launch solution and animation
#
  ode1:= diff(x__1(t),t,t)=P-m*g*sin(alpha):
  ics:=  x__1(0)=0, D(x__1)(0)=0:
  sol1:=dsolve( [ ode1, ics ] );
  s1:=eval( sol1, [ t=tau, x__1(t)=l]);
  s2:=eval( diff( sol1, t ), [ diff(x__1(t),t)=V__B, t=tau])/tau:
  sol2:=algsubs( rhs(s2)=lhs(s2), [s1, sol1]);
  tau__val:=eval(isolate( sol2[1], tau), params);
  sol3:=eval(sol2[2], [ params[], tau__val]):
  animate( pointplot,
           [ eval([(rhs(sol3)-l)*sin(alpha),(rhs(sol3)-l)*cos(alpha)], params),
             symbol=solidcircle,
             symbolsize=20
           ],
           t=0..rhs(tau__val),
           frames=20,
           trace=20
        );

x__1(t) = (1/2)*(P-m*g*sin(alpha))*t^2

  

l = (1/2)*(P-m*g*sin(alpha))*tau^2

  

[l = (1/2)*tau*V__B, x__1(t) = (1/2)*t^2*V__B/tau]

  

tau = 17.77777778

  

 
> 

#
# Post-launch solutioon and animation
#
  ode2:= diff(x(t),t,t)=0,
         diff(y(t),t,t)=-g:
  ics:= x(0)=0,
        y(0)=0,
        D(x)(0)=V__B*cos(alpha),
        D(y)(0)=V__B*sin(alpha):
  sol4:= dsolve( [ ode2, ics ] ):
  sol5:=eval(eliminate( eval( sol4, [ x(t)=d, y(t)=-h, t=T] ), T), params):
  h__val=evalf(rhs(isolate( sol5[2][], h)));
  animate( pointplot,
           [ eval([rhs~(sol4)[]], params),
             symbol=solidcircle,
             symbolsize=20
           ],
           t=0..rhs(sol5[1][]),
           frames=20,
           trace=20
        );

h__val = 1.174615859

  

 
> 

#
# It is possible to combine the two animations, but
#
# 1. the vertical scales are so different, (~40 pre-
#    launch, and 1 post-launch), and
# 2. the duration of the pre-launch and post-launch
#    periods are so different ( ~18secs pre-launch
#    and ~0.8secs post launch
#
# means that the combined animation doesn't "look"
# very good
#
  animate
  ( pointplot,
    [ piecewise
      ( t<rhs(tau__val),
        eval([(rhs(sol3)-l)*sin(alpha),(rhs(sol3)-l)*cos(alpha)], params),
        eval(subs( t=t-rhs(tau__val), eval([rhs~(sol4)[]], params)), params)
      ),
      symbol=solidcircle,
      symbolsize=20
    ],
    t=0..rhs(tau__val)+rhs(sol5[1][]),
    frames=100,
    trace=100
  );

 
> 

 

Download launch.mw

 

Use the big green up-arrow...

tomleslie 13876
May 27 2020
1 5

in the Mapleprimes toolbar to upload your worksheets

As you presented it I had to do a lot of "hacking" just to figure out what tou *might* want - and I *may* have got it wrong.

So all I can say is that the attached worksheet *might* perform the calculation(s) you want!!

> 

  restart;
  with(plots):
  hz2 := phi1(t)+arcsin((h+r1*sin(phi1(t)))/r2);
  deq := diff(phi1(t), t, t)-hz2 = 0;
#
# Two sets of initial conditions
#
  ics1 := phi1(0) = -arcsin(h/r1), (D(phi1))(0) = 0;
  ics2 := phi1(0) = 0, D(phi1)(0) = -0.63;

phi1(t)+arcsin((h+r1*sin(phi1(t)))/r2)

  

diff(diff(phi1(t), t), t)-phi1(t)-arcsin((h+r1*sin(phi1(t)))/r2) = 0

  

phi1(0) = -arcsin(h/r1), (D(phi1))(0) = 0

  

phi1(0) = 0, (D(phi1))(0) = -.63

 (1)
> 

#
# Solve/plot ode with initial conditions ic1
#
  sol1 := dsolve
          ( eval
            ( [ deq, ics1 ],
              [ r1 = 8, r2 = 8, r3 = 1, h = 5,
                m2 = 1, m3 = 20, theta3 = 20
              ]
            ),
            phi1(t),
            numeric
          ):
   odeplot( sol1,
            [ t, phi1(t) ],
            t=0..2
          );

 
> 

#
# Solve/plot ode with initial conditions ic2
#
   sol2:= dsolve
          ( eval
            ( [ deq, ics2 ],
              [ r1 = 8, r2 = 8, r3 = 1, h = 5,
                m2 = 1, m3 = 20, theta3 = 20
              ]
            ),
            phi1(t),
            numeric
          ):
   odeplot( sol2,
            [ t, phi1(t) ],
            t=0..2
          );

 
> 

 

Download solveODE.mw

Overkill?...

tomleslie 13876
May 26 2020
0 0

Using the DETools:-DEplot() may be a bit "overkill" if you just want the solution curves. See the attached for an alternative

> 

  restart;
  ode1 := diff(P(t), t) = 0.2*P(t) - 300:
  sol1:=dsolve([ ode1, P(0)=1300]);
  sol2:=dsolve([ ode1, P(0)=1800]);
  plot(rhs~([sol1, sol2]), t=0..20);

P(t) = 1500-200*exp((1/5)*t)

  

P(t) = 1500+300*exp((1/5)*t)

  

 
> 

s

s

 (1)
> 

 


 

Download odeplt.mw

Always tricky...

tomleslie 13876
May 26 2020
0 1

to figure out precisely what a poster needs when they describe a plot - but maybe as shown in the attached?

BTW - the "gridlines" will not appear in an actual Maple worksheet (unless you set the option gridlines=true). I have given up wonderiing why they always appear in any plots in worksheets uploaded to this site

> 

  l1:= plottools:-line( [-Pi/4,0], [-Pi/4, sec(-Pi/4)], color=black, thickness=5):
  l2:= plottools:-line( [ Pi/4,0], [ Pi/4, sec( Pi/4)], color=black, thickness=5):
  p:=  plots:-shadebetween(sec(x),0,x=-Pi/4..Pi/4, color=red):
  plots:-display([l1,l2,p], view=[-Pi/3..Pi/3, 0..1.1*sec(Pi/4)]);

 
> 

 


 

Download secplt.mw

Butchery...

tomleslie 13876
May 26 2020
1 0

The attached is a seriously "butchered" version of your code. I removed stuff wihch wasn't doing anything, combined some loops for efficiency, etc, etc - but the resulting version is still embarrassingly bad.

However it will perform the required calculations for any supplied value of the parameter 'alpha', returning a "2D" graph (technically a 3D spacecurve) and a 3D graph.

These spacecurves and surface plots can then be combined in any way you want. Several possibilities are shown in the attached.

> 

  restart:
  getPlots:= proc( alpha1, col1, col2)
                   local k:=2, M:=2, K:=2^(k-1),
                         N:=K*M, beta:=1, i, X, T, r, Phi_mxm,
                         key:=4, rho:=1, mmm:=1, sigma:=1,
                         lambda:=(m,alpha1,beta)->sqrt((2*m+alpha1+beta+1)*GAMMA(2*m+alpha1+beta+1)*m!/(2^(alpha1+beta+1)*GAMMA(m+alpha1+1)*GAMMA(m+beta+1))),
                         psii:=(n,m,x)-> piecewise( (n-1)/K <= x and x <= n/K,
                                                     2^(k/2)*lambda(m,alpha1,beta)*simplify(JacobiP(m,alpha1,beta,2^(k)*x-2*n+1)),
                                                     0
                                                  ),
                         w:=(n,x)->(1+x)^(1/gama-1),
                         omega:=(n,x)->w(n,2^(k)*x-2*n+1),
                         psi:=t->Matrix(N,1,[seq(seq(psii(i,j,t),j=0..M-1),i=1...K)] ),
                         PP:=alpha-> Phi_mxm.P(k,M,alpha).MatrixInverse(Phi_mxm),
                         u_exact:=(x,t)-> 1/(1+exp(sqrt(key/6)*x-(5/6)*key*t))^2,
                         f1:=unapply(u_exact(x,0),x),
                         g1:=unapply(u_exact(0,t),t),
                         g2:=unapply(u_exact(1,t),t),
                         U:=Matrix( N,N,symbol=u),
                         wwxx, wwt,ww,MC, PDEson, j, coz, UU,
                         P:= proc( k,M,nn )
                                   local PB,m,i,p,j,xi;
                                   m:=M*(2^(k-1)):
                                   xi:=(i,n)->((i+1)^(n+1)-2*i^(n+1)+(i-1)^(n+1));
                                   PB:=Matrix(m,m):
                                   for i from 1 to m do
                                       p:=0:
                                       for j from 1 to m do
                                           if   i=j
                                           then PB[i,j]:=1;
                                           fi:
                                           if   i<j
                                           then p:=p+1:
                                                PB[i,j]:= xi(p,nn);
                                           fi:
                                       end do;
                                       p:=1:
                                   end do:
                                   PB:=1/m^nn/GAMMA(nn+2)*PB;
                                   return PB;
                             end proc:
                   uses plots, LinearAlgebra;
                   for i from 1 to N do
                       X[i]:=evalf((2*i-1)/((2^k)*M)):
                       T[i]:=evalf((2*i-1)/((2^k)*M)):
                       r[i]:=evalf(psi(T[i])):
                   end do:

                   Phi_mxm:=Matrix([seq(r[i],i=1...N)]);
                   wwxx:= diff(f1(x),x,x)+(psi(x)^+.U.PP(alpha1).psi(t))(1):
                   wwt:= diff(g1(t),t)+x*(diff(g2(t),t)-diff(g1(t),t)-(psi(1)^+.PP(2)^+.U.psi(t))+psi(x)^+.PP(2)^+.U.psi(t))(1) :
                   ww:= f1(x)+g1(t)-g1(0)+x*(g2(t)-g2(0)-g1(t)+g1(0)-(psi(1)^+.PP(2)^+.U.PP(alpha1).psi(t)))(1)+(psi(x)^+.PP(2)^+.U.PP(alpha1).psi(t))(1) :
                   MC:=unapply(wwt -rho*wwxx -key*ww.(1-(1/mmm)*ww^sigma ),x,t):
                   PDEson:=Matrix(N, N, (i,j)-> simplify( evalf( MC(X[i],T[j] ) ))=0):
                   coz:=fsolve( { seq( seq(PDEson(i,j),i=1..N ),j=1..N )}):
                   U:=subs(op(coz),U):
                   UU:=(x,t)->evalf(f1(x)+g1(t)-g1(0)+x*(g2(t)-g2(0)-g1(t)+g1(0)-(psi(1)^+.PP(2)^+.U.PP(alpha1).psi(t)))(1)+(psi(x)^+.PP(2)^+.U.PP(alpha1).psi(t))(1)):
                   return plot3d(UU(x,t),x=0..1,t=0..1,color=col1, transparency=0.5),
                          plots:-spacecurve([0.5, t, UU(0.5,t)],t=0..1, color=col2, thickness=3):
             end proc:
> 

  p1,p2:=getPlots(1, red, blue):
  p3,p4:=getPlots(0.3, green, orange):
  plots:-display([p1,p2]);
  plots:-display([p3,p4]);
  plots:-display([p1,p3]);
  plots:-display([p2,p4]);

 

 

 

 
> 

 


 

Download grph23v2.mw

 

use a spacecurve...

tomleslie 13876
May 26 2020
1 1

as shown in the final group of the attached


 

> 

 

> 

restart:
with(orthopoly):
with(LinearAlgebra):
with(plots):
interface(rtablesize=20):
> 

k:=2:
M:=2:
K:=2^(k-1):
N:=K*M:

 

 
> 

alpha:=1:


alpha1:=alpha: 
beta:=1:
lambda:=(m,alpha1,beta)->sqrt((2*m+alpha1+beta+1)*GAMMA(2*m+alpha1+beta+1)*m!/(2^(alpha1+beta+1)*GAMMA(m+alpha1+1)*GAMMA(m+beta+1)));
psii:=(n,m,x)->piecewise((n-1)/K <= x and x <= n/K,  2^(k/2) *lambda(m,alpha1,beta)*simplify(JacobiP(m,alpha1,beta,2^(k)*x-2*n+1)), 0);
local i,j: psi:=(x)->Array([seq(seq(psii(i,j,x),j=0..M-1),i=1..K)] ):

w:=(n,x)->(1+x)^(1/gama-1):
omega:=(n,x)->w(n,2^(k)*x-2*n+1):

lambda := proc (m, alpha1, beta) options operator, arrow; sqrt((2*m+alpha1+beta+1)*GAMMA(2*m+alpha1+beta+1)*factorial(m)/(2^(alpha1+beta+1)*GAMMA(m+alpha1+1)*GAMMA(m+beta+1))) end proc

 

proc (n, m, x) options operator, arrow; piecewise((n-1)/K <= x and x <= n/K, 2^((1/2)*k)*lambda(m, alpha1, beta)*simplify(JacobiP(m, alpha1, beta, 2^k*x-2*n+1)), 0) end proc

(1.1)
> 


psi:=t->Matrix(N,1,[seq(seq(psii(i,j,t),j=0..M-1),i=1...K)] ):
 
 

> 


for i from 1 to N do
X[i]:=evalf((2*i-1)/((2^k)*M)):
end do:

for i from 1 to N do
T[i]:=evalf((2*i-1)/((2^k)*M)):
end do:

for i from 1 to N do
r[i]:=evalf(psi(T[i])):
end do:
Phi_mxm:=Matrix([seq(r[i],i=1...N)]);
 

Matrix(4, 4, {(1, 1) = 1.732050808, (1, 2) = 1.732050808, (1, 3) = 0., (1, 4) = 0., (2, 1) = -3.872983346, (2, 2) = 3.872983346, (2, 3) = 0., (2, 4) = 0., (3, 1) = 0., (3, 2) = 0., (3, 3) = 1.732050808, (3, 4) = 1.732050808, (4, 1) = 0., (4, 2) = 0., (4, 3) = -3.872983346, (4, 4) = 3.872983346})

(1.2)
> 

 

P:=proc(k,M,nn) local PB,m,i,p,j,xi;
m:=M*(2^(k-1)):
xi:=(i,n)->((i+1)^(n+1)-2*i^(n+1)+(i-1)^(n+1));
PB:=Matrix(m,m):
for i from 1 to m do
p:=0:
for j from 1 to m do
if i=j then PB[i,j]:=1;
 fi:
if i<j then p:=p+1:
PB[i,j]:= xi(p,nn);
 fi:
end do;
p:=1:
end do:
PB:=1/m^nn/GAMMA(nn+2)*PB;

return PB;
end proc:
PP:=alpha->Phi_mxm.P(k,M,alpha).MatrixInverse(Phi_mxm);

proc (alpha) options operator, arrow; `.`(Phi_mxm, P(k, M, alpha), LinearAlgebra:-MatrixInverse(Phi_mxm)) end proc

(1.3)
> 


key:=4:
rho:=1:
mmm:=1:
sigma:=1:

u_exact:=(x,t)-> 1/(1+exp(sqrt(key/6)*x-(5/6)*key*t))^2 ;
f1:=unapply(u_exact(x,0),x):
g1:=unapply(u_exact(0,t),t);
g2:=unapply(u_exact(1,t),t);

 
> 

U:=Matrix( N,N,symbol=u);
wwxx:= diff(f1(x),x,x)+(psi(x)^+.U.PP(alpha).psi(t))(1):
wwt:= diff(g1(t),t)+x*(diff(g2(t),t)-diff(g1(t),t)-(psi(1)^+.PP(2)^+.U.psi(t))+psi(x)^+.PP(2)^+.U.psi(t))(1) :
ww:= f1(x)+g1(t)-g1(0)+x*(g2(t)-g2(0)-g1(t)+g1(0)-(psi(1)^+.PP(2)^+.U.PP(alpha).psi(t)))(1)+(psi(x)^+.PP(2)^+.U.PP(alpha).psi(t))(1) :

MC:=unapply(wwt -rho*wwxx -key*ww .(1-(1/mmm)*ww^sigma ),x,t):

u_exact := proc (x, t) options operator, arrow; 1/(1+exp(sqrt((1/6)*key)*x-(5/6)*key*t))^2 end proc

  

g1 := proc (t) options operator, arrow; 1/(1+exp(-(10/3)*t))^2 end proc

  

g2 := proc (t) options operator, arrow; 1/(1+exp((1/3)*sqrt(6)-(10/3)*t))^2 end proc

  

Matrix(%id = 18446744074368602590)

 (1)
> 

PDEson:=Matrix(N,N):

for i from 1 to N do
for j from 1 to N do
PDEson(i,j):=simplify(
evalf(

MC(X[i],T[j])


)
)=0:
end do:
end do:
> 

coz:=fsolve({seq(seq(PDEson(i,j),i=1..N ),j=1..N )}):
U:=subs(op(coz),U):
 
> 


UU:=(x,t)->evalf(f1(x)+g1(t)-g1(0)+x*(g2(t)-g2(0)-g1(t)+g1(0)-(psi(1)^+.PP(2)^+.U.PP(alpha).psi(t)))(1)+(psi(x)^+.PP(2)^+.U.PP(alpha).psi(t))(1)):
> 


with(plots):
graphic_3D:=plot3d(UU(x,t),x=0..1,t=0..1,color=gray, transparency=0.5);

 
> 

graphic_2D:=plots:-spacecurve([0.5, t, UU(0.5,t)],t=0..1, color=green, thickness=3):
plots:-display([graphic_2D, graphic_3D]);

 
> 

 

 

 
> 

 

> 

> 

> 


 

Download grph23.mw

Ypu have two choices...

tomleslie 13876
May 25 2020
1 1
  1. either use evalc() which wii essentially simplify the expression based on the assumption that all unknown names are 'real, or
  2. use Maple's assume() command, to specify that both 'x' and 'lambda' are real

Both mehod's are shown in the attached

> 

restart;
conjugate(1-I*x*lambda^2)+1-I*x*lambda^2;
evalc(%);

2+I*conjugate(x*lambda^2)-I*x*lambda^2

  

2

 (1)
> 

restart:
interface(showassumed=0):
assume(lambda::real, x::real);
conjugate(1-I*x*lambda^2)+1-I*x*lambda^2;

2

 (2)
> 

 

Download comp.mw

 

Maple help facilty...

tomleslie 13876
May 25 2020
0 1

is worth reading!!!

Excerpts from the help at ?fsolve/details states

  • For a single polynomial equation of one variable, the fsolve command computes all real (non-complex) roots.

  • For a general equation or system of equations or procedures, the fsolve command computes a single real root.

Which part of these statements don't you understand?

 

 

Apply elementwise...

tomleslie 13876
May 25 2020
1 1

You can apply a function f() to a list 'L' of values, by using the elementwise operator '~', as in  f~(L).

See the attached for a couple of ways

> 

  restart:
  fprime_expr:=x^sin(x)*(cos(x)*ln(x)+sin(x)/x);
  RootFinding:-Analytic(diff(fprime_expr,x,x)=0, re=0..10, im=-1..1);
  L:=sort(select(type,[%],realcons));

  f:=x->x^sin(x);
#
# Apply 'f' elementwise to L
#
  f~(L);

x^sin(x)*(cos(x)*ln(x)+sin(x)/x)

  

4.78206109907178, 3.48122397295652, 2.16436178885414, .748344075670115-.240667154995574*I, .748344075670115+.240667154995574*I, 9.04887132811185, 6.77657190535075, 7.93089684792655

  

[2.16436178885414, 3.48122397295652, 4.78206109907178, 6.77657190535075, 7.93089684792655, 9.04887132811185]

  

proc (x) options operator, arrow; x^sin(x) end proc

  

[1.896584783, .6599753075, .2099102798, 2.475003072, 7.882490202, 2.244817923]

 (1)
> 

  restart:
  fprime_expr:=x^sin(x)*(cos(x)*ln(x)+sin(x)/x):
  getRoots:=proc( g::procedure)
                  local x__old:=0,
                        x__new:=[],
                        rt:
                  while true do
                        rt:= RootFinding:-NextZero(g, x__old):
                        if   rt <10
                        then x__new:=[x__new[], rt];
                             x__old:=x__new[-1];
                        else break;
                        fi;
                  od:
                  return x__new
            end proc:
  L:=getRoots( unapply( diff(fprime_expr,x,x), x) );
  f:=x->x^sin(x);
#
# Apply 'f' elementwise to L
#
  f~(L);   

[2.164361788, 3.481223972, 4.782061099, 6.776571905, 7.930896847, 9.048871328]

  

proc (x) options operator, arrow; x^sin(x) end proc

  

[1.896584783, .6599753083, .2099102798, 2.475003072, 7.882490202, 2.244817923]

 (2)
> 

 


 

Download fList.mw

First 74 75 76 77 78 79 80 Last Page 76 of 207