tomleslie

5473 Reputation

17 Badges

9 years, 338 days

MaplePrimes Activity


These are replies submitted by tomleslie

  1. It is possible to ask Maple to generate "power series" solutions for systems of ODEs. This is done in the attached, where I have generated the relevant series up to Order(12). However as the comparison graphs demonstrate, the power series solution deviates substantially from the numerical solution obtained by Runge-Kutta methods, as T increases. The attached shows that "reasonable" agreement between the two methods is only obtained up to ~6. One could increase the order of the power series solutions (the attached uses 12) and the range of "reasonable" will be extended.
  2. The code you supply in your worksheet does not implement the method required to produce a power series solutiion for an ODE system. In fact I have absolutely no idea what your code calculates, but I guarantee that it has nothing to do with power series solution of the ODE system.

If you really want to develop your own code for producing power series solutions, I suggest you read

  1. Chaper 3 of Elementary Differention Equations by Edwards  and  Penney, or
  2. https://en.wikipedia.org/wiki/Power_series_solution_of_differential_equations (for a quick introduction).

before starting

restart; with(plots); _local(gamma)

M__h := .50; beta__o := 0.34e-1; beta__1 := 0.25e-1; mu__r := 0.4e-3; sigma := .7902; alpha := .11; psi := 0.136e-3; xi := 0.5e-1; gamma := .7; M__c := .636; mu__b := 0.5e-2

ODEs := diff(B(T), T) = M__h-beta__1*psi*B(T)*G(T)-mu__r*B(T), diff(C(T), T) = beta__1*psi*B(T)*G(T)-sigma*psi*beta__1*E(T)*DD(T)-(alpha+xi+mu__r)*C(T), diff(DD(T), T) = alpha*C(T)-(gamma+mu__r)*DD(T), diff(E(T), T) = gamma*DD(T)+sigma*psi*beta__1*E(T)*G(T)-mu__r*E(T), diff(F(T), T) = M__c-psi*beta__o*F(T)*C(T)-mu__b*F(T), diff(G(T), T) = psi*beta__o*F(T)*C(T)-mu__b*G(T); bcs := B(0) = .50, C(0) = .30, DD(0) = .21, E(0) = .14, F(0) = .70, G(0) = .45

ans := dsolve([ODEs, bcs], numeric); for j in indets([ODEs], function(name)) do odeplot(ans, [T, j], T = 0 .. 30, title = convert(j, string), titlefont = [tims, bold, 20], thickness = 2, color = red) end do

 

 

 

 

 

 

interface(rtablesize = [22, 10]); M := Matrix([`~`[lhs](ans(0)), seq(`~`[rhs](ans(j)), j = 0 .. 20)])

Matrix(%id = 18446744074590329062)

(1)

#
# Generate power series solutions of order 12
# (Maple default is 6) for all ODEs
#
  Order:=12:
  ans2:= convert~
         ( evalf~
           ( dsolve
             ( {ODEs, bcs},
               indets([ODEs], function(name)),
               'series'
             )
           ),
           polynom
         );
#
# Plot the "numeric" solution obtained earlier
# and the equivalent power series solution
# obtaine above on th same graph.
#
# Note that the range of the plots is reduced,
# because even with 12-th order polynomials,
# the power series soluton and the numerical
# solution start to deviate for T>~6
#
  for j in indets([ODEs], function(name)) do
      display
      ( [ odeplot
          ( ans,
            [T, j],
            T = 0 .. 10,
            thickness = 2,
            color = red
          ),
          plot
          ( eval(j, ans2),
            T = 0 .. 10,
            thickness = 2,
            color=blue
          )
        ],
        title = cat(convert(j,string), " (Numerical and PowerSeries compared)"),
        titlefont = [tims, bold, 20]
      )
   end do

{B(T) = .5000000000+.4997992350*T-0.1003402817e-3*T^2+0.1470042520e-7*T^3-0.4186689729e-11*T^4+0.1797722185e-13*T^5-0.7652371692e-15*T^6+0.2923183302e-16*T^7-0.8844182475e-18*T^8+0.2215403975e-19*T^9-0.4589870880e-21*T^10+0.5389509965e-23*T^11, C(T) = .3000000000-0.4811931399e-1*T+0.3859529419e-2*T^2-0.2063396439e-3*T^3+0.8267055682e-5*T^4-0.2631388485e-6*T^5+0.6554216926e-8*T^6-0.5551714044e-10*T^7-0.1520371564e-10*T^8+0.2777703195e-11*T^9-0.3923000816e-12*T^10+0.4970214040e-13*T^11, DD(T) = .2100000000-.1140840000*T+0.3730565453e-1*T^2-0.8568110732e-2*T^3+0.1494601849e-2*T^4-0.2091819518e-3*T^5+0.2441368229e-4*T^6-0.2442660302e-5*T^7+0.2138541461e-6*T^8-0.1664279070e-7*T^9+0.1165691616e-8*T^10-0.7422668732e-10*T^11, E(T) = .1400000000+.1469441693*T-0.3995870043e-1*T^2+0.8709964152e-2*T^3-0.1500287681e-2*T^4+0.2093639084e-3*T^5-0.2441847491e-4*T^6+0.2442759169e-5*T^7-0.2138545263e-6*T^8+0.1664257448e-7*T^9-0.1165658906e-8*T^10+0.7422262550e-10*T^11, F(T) = .7000000000+.6324990290*T-0.1581608397e-2*T^2+0.2679492283e-5*T^3-0.6093270788e-8*T^4+0.1272029637e-9*T^5-0.4253521584e-11*T^6+0.1189725836e-12*T^7-0.2702521878e-14*T^8+0.3074973933e-16*T^9+0.3481499101e-17*T^10-0.6346661500e-18*T^11, G(T) = .4500000000-0.2249028960e-2*T+0.5983397429e-5*T^2-0.5345061622e-7*T^3+0.2810718704e-8*T^4-0.1239204116e-9*T^5+0.4250786124e-11*T^6-0.1189706297e-12*T^7+0.2702520657e-14*T^8-0.3074973865e-16*T^9-0.3481499102e-17*T^10+0.6346661500e-18*T^11}

 

 

 

 

 

 

 

 


 

Download odeRes2.mw

@Christian Wolinski 

Mine is set to [10,10] which is the default.

Until Maple2019 any matrix bigger than the interface(rtablesize) setting  wouldn't display at all (unless one reset the rtablesize). By default in Maple 2019 it now displays the matrix up to the rtablesize setting with the other entries elided, and just the little indicator at the bottom which tells you the actual matrix size - so the "15*9 Matrix" in  the output of my worksheet above

@Christian Wolinski 

One can get the output matrix directly

restart;
incog:=[theta[1, 1], theta[1, 2], theta[2, 1], theta[2, 2], theta[2, 6],
        theta[3, 0], theta[3, 3], theta[3, 4], theta[3, 5]
       ]:
eq:=[ 1, theta[1, 1]+theta[2, 2]+theta[3, 3], -theta[1, 1]-theta[2, 2],
      theta[2, 6]*theta[3, 5], -theta[1, 1]*theta[3, 3]-theta[2, 2]*theta[3, 3],
     -theta[1, 1]*theta[2, 6]*theta[3, 5]+theta[1, 2]*theta[2, 6]*theta[3, 4],
      theta[1, 1]*theta[2, 2]*theta[3, 3]-theta[1, 2]*theta[2, 1]*theta[3, 3]+theta[1, 2]*theta[2, 6]*theta[3, 0]
    ]:
myProc:= proc(P, V)
              local c, C, v;
              coeffs(P, V, 'c');
              seq([seq(`if`(type(C, dependent(v)),1,0), v = V)], C = c);
         end proc:
Matrix(map(myProc, eq, incog));

_rtable[18446744074332602366]

(1)

 


 

Download buildMat.mw

@acer 

I'm thinking about building a new workstation, using an AMD processor. (I've never gone the AMD route before - but you know - current price/performance etc!)

The fact that Maple itself appears to use the Intel MKL has now got me (slightly) worried, since (many years ago?)  this library was nicknamed "crippleAMD".

According to https://en.wikipedia.org/wiki/Math_Kernel_Library,

As of 2019, MKL, which remains the choice of many pre-compiled Mathematical applications on Windows (such as NumPy, SymPy, and MATLAB), still significantly underperforms on AMD CPUs with equivalent instruction sets

Would I really get a performance penalty if my (future) workstation uses an AMD processor??

@vv 

was actually much more interested in the discrepancy which arises between 'physics' versions

@Tegewaldt 

which command in the following you do not understand, or find "needlessly complicated"

  restart;
#
# Define a function
#
  f:= x-> 4*sin(x)+2*x^2+x+3;

proc (x) options operator, arrow; 4*sin(x)+2*x^2+x+3 end proc

(1)

#
# Return the value of this function
#
  f(p);

4*sin(p)+2*p^2+p+3

(2)

#
# Return the value of a function at some
# value of the independent variable
#
  f(2);

4*sin(2)+13

(3)

#
# Return the derivative of the function.
# Note that this derivative is itself a
# funcion
#
  D(f);

proc (x) options operator, arrow; 4*cos(x)+4*x+1 end proc

(4)

#
# Return the derivative of the funcion at a
# some value of the independent variable
#
  D(f)(2);

4*cos(2)+9

(5)

#
# Define a "new" function which contains the
# derivative of the original function
#
  g:= (x, y)-> D(f)(x)*y;

proc (x, y) options operator, arrow; (D(f))(x)*y end proc

(6)

#
# Return the value of the "new" function
#
  g(p,q);

(4*cos(p)+4*p+1)*q

(7)

#
# Return the value of the "new" function for
# a few values of the independent variables
#
  g(p,2);
  g(3,q);
  g(3,2);

8*cos(p)+8*p+2

 

(4*cos(3)+13)*q

 

8*cos(3)+26

(8)

 

Download basicDiff.mw

@Preben Alsholm 

On 64-bit Windows 7, I just updated to Physics 428.

The pop-up showed the usual progress bar, and the associated annotation reported "downloading"," installing" then "Done!".

Subsequent Physics:-Version() check in Maple gave the message

`The "Physics Updates" version in the MapleCloud is 428 and is the same as the version installed in this computer, created 2019, September 22, 7:8 hours, found in the directory C:\\Users\\TomLeslie\\maple\\toolbox\\2019\\Physics Updates\\lib\\`

as expected.

So I had no issues with this update

@minhhieuh2003 

was introducd in Maple 2015, and updated in Maple 2016. Exactly what version are you using?

@Jayaprakash J 

so you will have to specify the requirement more clearly.

Possibly??? the command

P:=unapply(doDiff(f,x,n),n);

@Melvin Brown 

This is the same error message which I obtained in a earlier post - and I suggested why this error message might occur. Try reading my earlier post

@Subhan331 

Questions on this site only get an answer in Maple code. If you want an answer in Matlab code, then post on

https://uk.mathworks.com/matlabcentral/answers

If you have a Maple problem in Maple code, then post the code here

"fminsearch" is a Matlab command to "Find minimum of unconstrained multivariable function using derivative-free method"

Why are you posting Matlab problems on a Maple forum??

"Exactly what is the problem???"

Your answer that " I think time-series result ( i mean graph) is incorrect?" isn't making it!!

On what mathematical basis do you think the solution provided by Maple is incorrect??

When I look at the 3D result provided by Maple, ie C(x,t) as a function of x and t, the only thing that looks  *slightly* suspicious is the region around C(1,0). The value at C(1,0) is forced to be 1 by the supplied ICs/BCs, but this value changes very rapidly (one might almost say "discontinuously") as (x,t) move even slightly from (1,0). This would make me double-check whether the relevant boundary/initial condition is "correct"

Otherwise I see no reason to doubt the result which Maple provides
 

The code in the image you supply is "cheating" becuase the expressions you supply for 'c' and 'd' mean two completely different things, and are certainly not equivalent

See the attached

a := (ui-uf)/lambda

(ui-uf)/lambda

(1)

b := delta*a

delta*(ui-uf)/lambda

(2)

c := expand(a*b)

delta*uf^2/lambda^2-2*delta*ui*uf/lambda^2+delta*ui^2/lambda^2

(3)

d := ui*delta*ui(1/lambda^2)-ui*delta*uf(1/lambda^2)-uf*delta*ui(1/lambda^2)+uf*delta*uf(1/lambda^2)

ui*delta*ui(1/lambda^2)-ui*delta*uf(1/lambda^2)-uf*delta*ui(1/lambda^2)+uf*delta*uf(1/lambda^2)

(4)

#
# Convert the above to 1-D input for illustration
#
# Note that the code contains the variable 'names'
# ui and uf, but also the 'functions' ui() and uf()
# with arguments (1/lambda^2).
#
# Maple consider the name 'ui' and the function 'ui()'
# to be two completely different entities
#
# The function definitions occur because there is
# no multiplication (explicit or implicit) in a
# subexpression like ui(1/lambda^2)
#
# Check what happens when these multiplications are
# added
#
  d := ui*delta*ui(1/lambda^2) - ui*delta*uf(1/lambda^2) - uf*delta*ui(1/lambda^2) + uf*delta*uf(1/lambda^2);
  d := ui*delta*ui*(1/lambda^2) - ui*delta*uf*(1/lambda^2) - uf*delta*ui*(1/lambda^2) + uf*delta*uf*(1/lambda^2);

ui*delta*ui(1/lambda^2)-ui*delta*uf(1/lambda^2)-uf*delta*ui(1/lambda^2)+uf*delta*uf(1/lambda^2)

 

delta*uf^2/lambda^2-2*delta*ui*uf/lambda^2+delta*ui^2/lambda^2

(5)

 


 

Download notEquiv.mw

 

@acer 

"You could try testing whether StringTools:-Search returns a result greater than 1."

I think I'd check whether StringTools:-Search returns a result greater than 0

3 4 5 6 7 8 9 Last Page 5 of 150