tomleslie

"Elementwise"...

Answered: tomleslie 13876

May 17 2020

0 0

Maple commands which accept a single argument can be made to apply to all elements of a "container" (eg matrix,table array, list etc) by using the elementwise operator. So in my earlier worksheet,

parse(ImportMatrix(cat(fpath, "test.txt"), datatype=string))

would fail, because the parse() command expects its argument to be a "string" not a matrix. On the other hand

parse~(ImportMatrix(cat(fpath, "test.txt"), datatype=string))

succeeds because the parse() function is being applied tto the elements of the matrix, which are all strings.

You can achieve the same effect by 'mapping' the command parse() over the elements of a "container" with

map( parse, ImportMatrix(cat(fpath, "test.txt"), datatype=string))

Maybe...

Answered: tomleslie 13876

May 17 2020

0 0

@J F Ogilvie

If you upload the

" valuable worksheet in Maple 4.3 by 'Steve Adams' "

to this site (using the big green up-arrow in the Mapleprimes toolbar), then maybe soemone here will be able to determine its purpose, and make it executable in modern-day Maple

I doubt it...

Answered: tomleslie 13876

May 12 2020

0 0

@ogunmiloro

I think that obtaining an analytic solution to pdes which contain differentials raised to powers like 0.01 is highly unlikely. If you can find any example, in any textbook, anywhere where a pde containing a differential raised to a floating point power (like 0.01) could be solved analytically - I'd like to see it!

I have tried it for your problem - and it "fails" with an error message. Now it probably shouldn't fail with an "error" - it should perform the "graceful" fail of returning unevaluated, signifying that no analytic soltuoin can be obtained.

If you want to try it yourself then just add the code

#
# Generate solutions for three different
# values of alpha__1
#
  alphaVal:=[3.0,4.0,5.0]:
  for j from 1 by 1 to 3 do
      aVal:=alphaVal[j]:
      sol[j]:=pdsolve([eval( pdes, pars)[], bcs[]]):
  od:

to the bottom of the worsheet I poste previously

An illustration...

Answered: tomleslie 13876

May 12 2020

0 0

@ogunmiloro

This shows examples of how you can achieve your requirements

>

restart;

>	pars:= [ sigma__1 = 0.110^(-5), alpha__1 = aVal, k__1 = 1.1, k__2 = 0.01, sigma__2 = 0.110^(-5), alpha__2 = 4.0, k__3 = 0.01, k__4 = 0.12, beta__1 = 0.2, beta__2 = 0.02 ]:

>	pdes := [ diff(C(t, x), t) = sigma__1diff(C(t, x), x, x) + alpha__1C(t, x)^k__1 + alpha__1C(t, x)^k__2B(t, x)^k__3 - beta__1C(t, x), diff(B(t, x), t) = sigma__2diff(B(t, x), x, x) + alpha__2B(t, x)^k__3 + alpha__2C(t, x)^k__2B(t, x)^k__4 - beta__2B(t, x)];

[diff(C(t, x), t) = sigma__1*(diff(diff(C(t, x), x), x))+alpha__1*C(t, x)^k__1+alpha__1*C(t, x)^k__2*B(t, x)^k__3-beta__1*C(t, x), diff(B(t, x), t) = sigma__2*(diff(diff(B(t, x), x), x))+alpha__2*B(t, x)^k__3+alpha__2*C(t, x)^k__2*B(t, x)^k__4-beta__2*B(t, x)]

(1)

>	bcs := [ D[2](C)(t, 0) = 0, D[2](C)(t, 1) = 0, D[2](B)(t, 0) = 0, D[2](B)(t, 1) = 0, C(0, x) = 200.0, B(0, x) = 300 ];

(2)

>	# # Generate solutions for three different # values of alpha__1 # alphaVal:=[3.0,4.0,5.0]: for j from 1 by 1 to 3 do aVal:=alphaVal[j]: sol[j]:=pdsolve(eval( pdes, pars), bcs, numeric): od:

>	# # Plot eg C(t,x) for three different values of alpha__1 # p1 := (sol[1]):-plot3d(C(t, x), t = 0..0.5, x = 0..0.5); p2 := (sol[2]):-plot3d(C(t, x), t = 0..0.5, x = 0..0.5); p3 := (sol[3]):-plot3d(C(t, x), t = 0..0.5, x = 0..0.5);

>

#
# Generate a matrix of results for each of the
# three values of alpha__1
#
  for j from 1 by 1 to 3 do
      a:=sol[j]:-value( output=procedurelist);
      numRes[j]:=Matrix( [ [ `t`, `x`, `B(t,x)`, `C(t,x)`],
                           seq( seq( rhs~(a(i, j)), i=0..0.5,0.1), j=0..0.5,0.1)
                         ]
                       );
  od:
  interface(rtablesize=[50,10]):
  numRes[1];
  numRes[2];
  numRes[3];
  interface(rtablesize=[10,10]):

Matrix(%id = 18446744074370352782)

(3)

>

Download pde_prob2.mw

Makes no sense...

Answered: tomleslie 13876

May 11 2020

0 0

@666 jvbasha

NOne of this response makes much sense to me, because

Your original worksheet iteratively computes a bivariate polynomial in 'x' and 'h'. I have no idea what the significance of this bivariate polynomial actually is
Your original worksheet "converts" this bivariate polynomial in the variables 'x' and 'h' to a univariate polynomial in 'h', by the simpe expedient of setting x=1
Having generated the univariate polynomial in the variable 'h', you then plot it
All of the above is fine - and in my response all of the above now works
Now you are making comments about
1. "residual error". Neither a bivariate nor a univariate polynomial has any concept of a "residual error"
2. "solution interval". Neither a bivariate nor a univariate polynomial has any concept of a "solution interval"
3. "integral interval". Neither a bivariate nor a univariate polynomial has any concept of an "integral interval"

Now given that all you have is a bivariate polynomial and a univariate polynomial - you might want to define what you mean by

"residual error" for a bivariate/univariate polynomial
"solution interval" for a bivariate/univariate polynomial
"integral interval". for a bivariate/univariate polynomial

By the way I have also read your comment

Where the curves are a straight line(-2..2), that value (eg. h=(-1) in hh) we can take for velocity plot.

and concluded that it is completely meaningless - Given that you have have a bivariate/univariate polynomial, where are the "curves" (what curves?) a straight line? And wtf is a "velocity plot"?

Don't have this particular package insta...

Answered: tomleslie 13876

May 09 2020

0 0

@rlopez

but I have installed/used this one https://www.maplesoft.com/applications/view.aspx?SID=33406, which covers a variety of orthogonal expansions

In the attached, it appears to cover what the OP requires - as in basically generate a Fourier serie representation fpr a sawtooth

>

  restart;
#
# Define a sawtooth
#
  f:= x->x/2-round(x/2);
#
# Generate the Fourier series for the sawtooth
#
  Trig:= OrthogonalExpansions:-FourierSeries
                               ( f(x),
                                 x=-1..1,
                                 n,
                                 `Coefficients`
                               );
#
# Plot the sawtooth and its successive Fourier
# series approximations
#
  plots:-animate
         ( plot,
           [ [ f(x), Trig],
             x=-4..4,
             color= [red, blue],
             title=typeset("N=%1",round(n)),
             titlefont=[times, bold, 20]
           ],
           n=1..50,
           paraminfo=false
         );

Sum(-(-1)^i*sin(Pi*i*x)/(Pi*i), i = 1 .. n)

>

Download fourStuff.mw

As procedures...

Answered: tomleslie 13876

May 06 2020

0 0

I should have noted that you want procedures which accept any argument for number and base, so revise versions of the baove which has all three methods as procedures

>	# # The really easy way # convert(3657, base, 8);

(1)

>	# # And if you prefer "most-significant" first # (a matter of personal taste, I admit) then # you could simply reverse the above # ListTools:-Reverse(convert(3657, base, 8))

(2)

>	# # Or if you feel under some sort of "obligation" # to write "code", then you could use # q:=3657: j:=1: while q>0 do a[j]:=irem(q, 8,'q'); j:=j+1; od: convert(a, list);;

(3)

>

#
# Writing each of the above as procedures which will
# accept any number, and any base
#
  bc1:= (x, b)->convert(x, base, b):

  bc2:= (x, b)->ListTools:-Reverse(convert(x, base, b)):

  bc3:= proc(num::integer, bas::integer)
             local j:=1, q:=num, a:
             while q>0 do
                   a[j]:=irem(q, 8,'q');
                   j:=j+1;
             od:
             return convert(a, list);
       end proc:

  bc1(3657, 8);
  bc2(3657,8);
  bc3(3657,8);

(4)

>

Download base2.mw

More observations...

Answered: tomleslie 13876

May 01 2020

0 0

@Carl Love

"Floating point subscripts" I can make this behaviour "go away" if I export the OP's original as a .mpl file and then read() the .mpl file in a new worksheet. I only tried this on the off-chance that there might be some "non-printing" character in the original, which would either break the export process, or be identifiable in an external editor. However the export process worked as normal: and the .mpl file showed no sign of non-ascii characters in an external editor - so pretty inconclusive.
Why does the revised worksheet oddEqs.mw I posted obtain a non-zero result? This seems to depend on the precise "preprocessing" of the equations prior to the solve() command. OP uses combine(expand(...)) where as my original worksheet omitted the combine(..) step. Since I was using solve(), to provide an exact solution, I thought the combine() step was redundant, although as the attached shows - it does make a significant difference to the solution(s) obtained. I then thought this might be some floating-point-error phenomenon (might still be??), but changing the 'Digits' setting doesn't seem to make a difference - without the combine() step, a non-zero answer is obtained: with the combine() step only a zero answer is obtained. See the attached.

>

restart:

>

  odeA1:=   0.0001995000000*a__1*omega^2*cos(omega*t)
          + 0.0001995000000*a__3*omega^2*cos(3*omega*t)
          + 0.004987500000*a__5*omega^2*cos(5*omega*t)
          - 93.86*a__1*cos(omega*t)
          - 10.42888889*a__3*cos(3*omega*t)
          - 93.86*a__5*cos(5*omega*t)
          + 1.633844046*10^18*(A__1*cos(omega*t)
          + A__3*cos(3*omega*t)/9
          + A__5*cos(5*omega*t))^2*(a__1*cos(omega*t)
          + a__3*cos(3*omega*t)/9 + a__5*cos(5*omega*t))
          + 0.02250237053*a__1*omega*sin(omega*t)
          + 0.007500790177*a__3*omega*sin(3*omega*t)
          + 0.1125118526*a__5*omega*sin(5*omega*t):
  odeA2:=   A__1*omega*sin(omega*t)
          + A__3*omega*sin(3*omega*t)/3
          + 5*A__5*omega*sin(5*omega*t)
          - 1.499250375*10^(-7)*(A__1*cos(omega*t)
          + A__3*cos(3*omega*t)/9
          + A__5*cos(5*omega*t))*(2.110712380*10^10
          - 1.633844050*10^18*(a__1*cos(omega*t)
          + a__3*cos(3*omega*t)/9
          + a__5*cos(5*omega*t))^2)
          + 1.499250375*10^(-7):
#
# With the 'combine()' step
#
  sol1:= eval
         ( solve
           ( [ seq
               ( coeff~
                 ( combine~
                   ( expand~
                     ( [odeA1,odeA2]
                     )
                   ),
                   cos(j*omega*t)
                 )[],
                 j=1..15,2
               )
             ],
             [ a__1, a__3, a__5, A__1, A__3, A__5]
           )[],
          omega=485.0135000
        );
#
# Without the 'combine()' step
#
  sol2:= eval
         ( solve
           ( [ seq
               ( coeff~
                 ( expand~
                   ( [odeA1,odeA2]
                   ),
                   cos(j*omega*t)
                 )[],
                 j=1..15,2
               )
             ],
             [ a__1, a__3, a__5, A__1, A__3, A__5]
           )[],
           omega=485.0135000
         );

(1)

>

Download oddEqs2.mw

I don't think you understood...

Answered: tomleslie 13876

April 30 2020

0 0

@erik10

my original answer.

No-one apart from you has reported a problem opening .mw files in Maple 2020 in Windows 10
Quite a lot of people seem to have a problem with broken file associations in Windows 10. This is a Windows 10 problem
So somehow your Windows installation won't allow/accept/retain/setup a file association between Maple 2020 and .mw files
Can only suggest you Google "Problems with Windows 10 file associations" and start ploughing through the many suggestions

The solution is complex...

Answered: tomleslie 13876

April 28 2020

0 0

@raj2018

and by "complex" I don't mean "complicated", but (literally) it generates complex numbers. Now Maple can solve IVP problems whihc generate complex solutions - but pretty much everything gets more complicated, very quickly.

A trivial instance of the complication - since solutions of A(x) are complex, what does one choose as a suitable initial condition, since A(0) is now allowed (required?) to be a complex number?

The attached shows a couple of solution cases with A(x) complex, for real initial conditions.

Before I go any further with this, a number of questions need answering, such as

Are you expecting A(x) to be complex? If not, why not?
Since A(x) appears to be generally complex, what is an acceptable value for the (possibly complex) initial value A(0). Whether (or not) numerical solutions can be obtained, seems to be very dependent on this initial condition
If an ODE is providing solutions that you really (really) do not expect, this *usualll* means
1. The problem is incorrectly set up - maybe some typo somewhere - just crawl all over everything
2. You are (without realising) making some kind of assumptions - which of course Maple will not do. So think long and hard about why you are assuming certain kinds of solutions

>

  restart:
  q1:=0.4: q:=0.4: s1:=0: mu:=2: mu1:=0.7:
  mu3:=0.2: s:=0.5: M:=1.9109: mu2:=1+mu3-mu1:

  V:=unapply( M^(2)-M*(M^(2)+2*A)^(1/(2))+mu1*(M^(2)*mu+s1)-mu1*((M^(2)*mu)/(2))^(1/(2))*(((M^(2)*mu+3*s1-2*A)+sqrt((M^(2)*mu+3*s1-2*A)^(2)-12*M^(2)*mu*s1))^(1/(2))+(4*M^(2)*mu*s1)*((M^(2)*mu+3*s1-2*A)+sqrt((M^(2)*mu+3*s1-2*A)^(2)-12*M^(2)*mu*s1))^(-3/(2)))-(2*mu2)/((3*q-1))*((1-(q-1)*A)^((3*q-1)/(2*q-2))-1)-(2*mu3)/(s*(3*q1-1))*((1+s*(q1-1)*A)^((3*q1-1)/(2*q1-2))-1),A):

# ode:=diff(A(x),x)^2+V(A(x))=0;
#
# Manually rewrite ODE selecting one sign of the
# square root
#
  ode:=diff(A(x),x)=sqrt(V(A(x))):
  odeSol:=dsolve([ode, A(0)=1], numeric, complex=true):
#
# Plot the real and imaginary parts
#
  plots:-odeplot( odeSol, [ [x, Re(A(x))],[x, Im(A(x))] ], x=0..10, color=[red,blue]);

>	# # Now use the other sign of the square roo # ode:=diff(A(x),x)=-sqrt(V(A(x))): odeSol:=dsolve([ode, A(0)=-1], numeric, complex=true): # # Plot the real and imaginary parts # plots:-odeplot( odeSol, [ [x, Re(A(x))],[x, Im(A(x))] ], x=-10..0, color=[red,blue]);

>

Download odeComplex.mw

In the attached...

Answered: tomleslie 13876

April 28 2020

0 0

@AHSAN

I have written three (very simple) rootfinding procedures, using three conceptually different iteration processes, namely

regula falsi
bisection
newton's method

>

  restart;
#
# Procedure which computes 'root' using regula falsi,
# returns root, residual, and number of iterations
#
  regFals:= proc( g, X__A, X__B, err)
                  local j, p, X__1:=X__A, X__2:=X__B, iters:=0:
                  uses Student[Precalculus]:
                #
                # Iterate until tolerance is achieved
                #
                  for j from 1 do
                      iters:=iters+1:
                    #
                    # Get the x-intercept of the line between
                    # the two points given by {X__1, f(X__1)]
                    # and [X__2, f(X__2)]
                    #
                      p:=evalf( [ Line
                                  ( [X__1, g(X__1)],
                                    [X__2, g(X__2)]
                                  )
                                ][4]
                              );
                    #
                    # Exit if required tolerance has be achieved
                    #
                      if   abs(g(p))<err
                      then break
                    #
                    # Otherwise update the range end-points
                    #
                      elif g(p)>0
                      then X__2:=p:
                      else X__1:=p;
                      fi;
                 od:
                 return root=p, residual=g(p), numberOfIterations=iters;
            end proc:
#
# Procedure which computes root using the bisection method
#
  bisection:= proc( g, X__A, X__B, err )
                    local j, p, X__1:=X__A, X__2:=X__B, iters:=0:
                  #
                  # Iterate until tolerance is achieved
                  #
                    for j from 1 do
                        iters:= iters+1:
                      #
                      # Get updated x-value
                      #
                        p:= (X__1+X__2)/2.0:
                      #
                      # Exit if required tolerance has be achieved
                      #
                        if   abs(g(p)) < err
                        then break
                      #
                      # Otherwise update the range end-points
                      #
                        elif g(p)>0
                        then X__2:=p:
                        else X__1:=p;
                        fi;
                   od:
                   return root=p, residual=g(p), numberOfIterations=iters;
              end proc:
#
# Procedure wich computes root using Newton's method
#
  newton:= proc( g, X__A, X__B, err )
                 local j, p, X__1:=X__A, X__2:=X__B, iters:=0:
               #
               # Iterate until tolerance is achieved
               #
                 for j from 1 to 50 do
                     iters:= iters+1:
                   #
                   # Get updated x-value
                   #
                     p:= evalf(X__1-eval(g(X__1)/diff(g(x),x), x=X__1)):
                   #
                   # Exit if required tolerance has be achieved
                   #
                     if   abs(g(p)) < err
                     then break
                   #
                   # Otherwise update the
                   #
                     else X__1:=p;
                     fi;
                od:
                return root=p, residual=g(p), numberOfIterations=iters;
          end proc:

>	# # Toy example # f:=x->x*exp(x)-4: X__1:=0: X__2:=2: tol:=1.0e-06: plot( f(x), x=0..2); # # Try each of these three methods # regFals(f, X__1, X__2, tol); bisection(f, X__1, X__2, tol); newton(f, X__1, X__2, tol);

(1)

>

Download rtFinf.mw

I believe...

Answered: tomleslie 13876

April 26 2020

0 0

@AHSAN

I have already fixed this elsewhere in this thread

Fixed a few syntax errors...

Answered: tomleslie 13876

April 26 2020

0 0

@AHSAN

and restricted the search range for the root to a region where the function returns real (rather than complex) values, produces the following

>

(1)

>

F := (-(288*(lambda+2*k*(1/3)-2/3))*sqrt(2)*lambda^2*arctan((1/2)*sqrt(4*lambda/(1-k)-2)*sqrt(2))/(1-k)^2-(144*(lambda+2*k*(1/3)-2/3))*Pi*lambda^2*sqrt(2)/(1-k)^2-(36*((lambda+2*k*(1/3)-2/3)*(4*lambda/(1-k)-2)+10*lambda*(1/3)+4*k*(1/3)-4/3))*sqrt(4*lambda/(1-k)-2))*(1-k)^2/(256*lambda^2)

0.3164062500e-2*(-355.5555555*(lambda-.6000000000)*2^(1/2)*lambda^2*arctan((1/2)*(4.444444444*lambda-2)^(1/2)*2^(1/2))-558.5053607*(lambda-.6000000000)*lambda^2*2^(1/2)-36*((lambda-.6000000000)*(4.444444444*lambda-2)+(10/3)*lambda-1.200000000)*(4.444444444*lambda-2)^(1/2))/lambda^2

(2)

>

proc (lambda) options operator, arrow; 0.3164062500e-2*(-355.5555555*(lambda-.6000000000)*2^(1/2)*lambda^2*arctan((1/2)*(4.444444444*lambda-2)^(1/2)*2^(1/2))-558.5053607*(lambda-.6000000000)*lambda^2*2^(1/2)-36*((lambda-.6000000000)*(4.444444444*lambda-2)+(10/3)*lambda-1.200000000)*(4.444444444*lambda-2)^(1/2))/lambda^2 end proc

(3)

>

seq(evalf(F(j)), j = .1 .. 1, .1)

(4)

>

tol := 0.10e-5; `λ__1` := .5; `λ__2` := 1.0

(5)

>

for j do p := evalf([Line([`λ__1`, F(`λ__1`)], [`λ__2`, F(`λ__2`)])][4]); if abs(evalf(F(p))) < tol then break elif evalf(F(p)) > 0 then `λ__2` := p else `λ__1` := p end if end do; root = p, residual = evalf(F(p))