Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 169 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

Not sure...

tomleslie 13876
February 08 2020
1 1

that a "closed-form" recursion relation could be obtained,

The attached will calculate the functions b() and c() reasonably efficiently

> 

  restart:
> 

  b := proc( k::nonnegint)
             option remember:
             return(-1)^k*GAMMA(1 + 1/2*k)^2/k!;
       end proc:
  c:= proc( k::nonnegint )
            option remember:
            if   k=0
            then return b(0)^N
            else return add( (q*(N + 1) - k)*b(q)*c(k - q),
                             q = 1 .. k
                           )/(k*b(0));
            fi;
      end proc:
> 

#
# values of b()
#
  b(1)^N;
  b(2)^N;
  b(3)^N;
  b(4)^N;
#
# values of c()
#
  collect(expand(c(1)),N);
  collect(expand(c(2)),N);
  collect(expand(c(3)),N);
  collect(expand(c(4)),N);

(-(1/4)*Pi)^N

  

(1/2)^N

  

(-(3/32)*Pi)^N

  

(1/6)^N

  

-(1/4)*N*Pi

  

(1/32)*Pi^2*N^2+(-(1/32)*Pi^2+1/2)*N

  

-(1/384)*Pi^3*N^3+((1/128)*Pi^3-(1/8)*Pi)*N^2+(-(1/192)*Pi^3+(1/32)*Pi)*N

  

(1/6144)*Pi^4*N^4+(-(1/1024)*Pi^4+(1/64)*Pi^2)*N^3+((11/6144)*Pi^4-(3/128)*Pi^2+1/8)*N^2+(-(1/1024)*Pi^4+(1/128)*Pi^2+1/24)*N

 (1)

 

> 

 

NULL

Download recur.mw

As a general rule...

tomleslie 13876
February 07 2020
0 0

No-one her will respond to an "obvious" homework question unless the poster show some sign of actually having attempted the question and posted their code.

So all I am prepared t  say is that you should investigate the StateSpace() command in the DynamicSystems package. This will allow you to generate a standard state-space descriiption (the conventional A, B, C, D matrices) simply by entering your differetnial equation.

If I remember my 'control theory' you can then derive the state transition matrix from the 'A' matrix in the state-space description, although I think you will also need four inital state conditions.

A caveat. Even with appropriate initial conditions, I suspect that the state transition matrix will be incredibly complicated for symbolic parameters a__0, a__1, a__2, a__3, in your ODE. Given a set of initial conditons and numeric values for the parameters, the problem ought to be pretty trivial

Have you tried the page...

tomleslie 13876
February 06 2020
0 0

shown below?

https://www.maplesoft.com/support/downloads/index.aspx

What do you hope to achieve??!!!...

tomleslie 13876
February 06 2020
1 2

because I do not believe there is any justification for just ignoring the non-linear term in an ODE and hoping that the solution to the corresponding linear ODE will bear any relation that whihc would be obtained from the non-linear equation - so what is the point???

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

Having said that, the attached worksheet attempts to solve the ODE 'eq1' for the parameter set

a=0.322,  b=2.05e-05,  y0,  x0=5,  x1=-0.6

where y0=0.1 to 0.2 in steps of 0.001.

The dsolve(eq1...,numeric) command fails every time with a singularity warning. The x-value of the singularity varies slightly.

When  dsolve(eq1...,numeric) fails, the 'try-catch' construct in the attached does the following for all y0 values in the range 0.1 to 0.2 in steps of 0.001.

  1. Obtains the last value of 'x' for which eq1 (the non-linear ODE) generates a result
  2. Evaluates/stores  x, y(x), diff(y(x),x) at the value determined in (1) above
  3. Switches to using eq2,  (the linear ODE) and then
    1. Uses eq2 to evaluate/store x, y(x), diff(y(x),x) at the value determined in (1) above
    2. Uses eq2 to evaluate/store x, y(x), diff(y(x),x) at x=-0.6, which is the "final" value in the desired range. This is just to demonstrate that the numerical solution of eq2 has been performed over the requested range

On completion of the for-loop containing the try-catch construct, the attached

  1. Prints the values obtained at (2) above
  2. Prints the values obtained at (3)(1) above
  3. Prints the values obtained at (3)(2) above
> 

  restart;
  eq1:=(a,b)->-(2/3)*b*(1+y(x))*(1+x)*(1+a*(1+x)^3)^2*(diff(y(x), x, x))+(8/9)*b*(1+a*(1+x)^3)^2*(diff(y(x), x))^2-(1/3)*b*(1+a*(1+x))*(1+y(x))*a*(1+a*(1+x)^3)*(diff(y(x), x))+y(x)*(1+y(x))^2*(1+x)^2*a^2*(2+a*(1+x)):
  eq2:=(a,b)->-(2/3)*b*(1+a*(1+x)^3)^2*(1+x)*(diff(y(x), x, x))-(1/3*(1+a*(1+x)))*a*(1+a*(1+x)^3)*b*(diff(y(x), x))+a^2*y(x)*(2+a*(1+x))*(1+x)^2:
  s1:=(a, b, y0, x0, x1)->dsolve( { eq1(a,b), y(x0)=y0, (D(y))(x0)=-y0/(1+x0)}, type=numeric):
  s2:=(a, b, y0, x0, x1)->dsolve( { eq2(a,b), y(x0)=y0, (D(y))(x0)=-y0/(1+x0)}, type=numeric):
> 

#
# Set up some arrays for results
#
  A:=Array(0..11):
  B:=Array(0..11):
  C:=Array(0..11):
  for i from 0 by 1 to 11 do
    #
    # Generate solution modules for each value
    # of y0
    #
      y0:= 0.01+i*0.001:
      sol1:= s1(0.322, 2.05e-05, y0, 5, -0.6):
      sol2:= s2(0.322, 2.05e-05, y0, 5, -0.6):
      try #
          # Attempt to generate the solution of eq1
          # at the left hand end of the desired range.
          #
          # For all values of y0 above, this *always*
          # fails
          #
          A[i]:= sol1(-0.6)
      catch "cannot evaluate":
          #
          # Store values at the singularity obtained from eq1
          #
            A[i]:= [ `sol1`, y0, rhs~(sol1([lastexception][3])) ]:
          #
          # Store values at the singularity obtained from eq2
          #
            B[i]:= [ `sol2`, y0, rhs~(sol2([lastexception][3])) ]:
          #
          # Store values at the left-hand end of the desired
          # range obtained from eq2
          #
            C[i]:= [ `sol2`, y0, rhs~(sol2(-0.6)) ]:
      end try:
  od:
> 

  printf( "     ***********************************************************\n"):
  printf( "     *                                                         *\n"):
  printf( "     * Values at the singularity by solving eq1                *\n"):
  printf( "     *                                                         *\n"):                        
  printf( "     ***********************************************************\n"):

  printf("        y0          x(sing)           y(sing)         dy/dx(sing)\n");
  seq( printf("     %6.1e%18.6e%18.6e%18.6e\n",
              A[i][2], seq(A[i][3][j], j=1..3)
             ),
       i=0..10
     );
  printf(" \n"):
  printf( "     ***********************************************************\n"):
  printf( "     *                                                         *\n"):
  printf( "     * Values at the singularity by solving eq2                *\n"):
  printf( "     *                                                         *\n"):                        
  printf( "     ***********************************************************\n"):
  printf("       y0           x(sing)           y(sing)         dy/dx(sing)\n");
  seq( printf("     %6.1e%18.6e%18.6e%18.6e\n",
              B[i][2], seq(B[i][3][j], j=1..3)
             ),
       i=0..10
     );
  printf(" \n"):
  printf( "     ***********************************************************\n"):
  printf( "     *                                                         *\n"):
  printf( "     * Values at the left hand end (x=-0.6) by solving eq2     *\n"):
  printf( "     *                                                         *\n"):                        
  printf( "     ***********************************************************\n"):
  printf("        y0          x(-0.6)           y(-0.6)         dy/dx(-0.6)\n");
  seq( printf("     %6.1e%18.6e%18.6e%18.6e\n",
              C[i][2], seq(C[i][3][j], j=1..3)
             ),
       i=0..10
     );

     ***********************************************************
     *                                                         *
     * Values at the singularity by solving eq1                *

     *                                                         *
     ***********************************************************
        y0          x(sing)           y(sing)         dy/dx(sing)
     1.0e-02      4.039683e+00     6.470087e+012    -1.305346e+020
     1.1e-02      4.050764e+00     2.191077e+013    -8.093599e+020
     1.2e-02      4.060928e+00     8.575603e+012    -1.972576e+020
     1.3e-02      4.070320e+00     1.037853e+014    -8.269505e+021
     1.4e-02      4.079050e+00     1.047595e+013    -2.641427e+020
     1.5e-02      4.087209e+00     7.699123e+013    -5.243228e+021
     1.6e-02      4.094868e+00     1.246622e+013    -3.404310e+020
     1.7e-02      4.102086e+00     4.503155e+013    -2.329597e+021
     1.8e-02      4.108913e+00     2.215896e+013    -8.016540e+020
     1.9e-02      4.115390e+00     4.119879e+014    -6.407932e+022
     2.0e-02      4.121552e+00     1.632731e+013    -5.041440e+020
 
     ***********************************************************
     *                                                         *
     * Values at the singularity by solving eq2                *
     *                                                         *
     ***********************************************************
       y0           x(sing)           y(sing)         dy/dx(sing)
     1.0e-02      4.039683e+00     4.726846e+000    -4.064049e+001
     1.1e-02      4.050764e+00     4.728147e+000    -4.045244e+001
     1.2e-02      4.060928e+00     4.729303e+000    -4.028074e+001
     1.3e-02      4.070320e+00     4.730344e+000    -4.012283e+001
     1.4e-02      4.079050e+00     4.731282e+000    -3.997657e+001
     1.5e-02      4.087209e+00     4.732139e+000    -3.984045e+001
     1.6e-02      4.094868e+00     4.732920e+000    -3.971307e+001
     1.7e-02      4.102086e+00     4.733641e+000    -3.959344e+001
     1.8e-02      4.108913e+00     4.734304e+000    -3.948062e+001
     1.9e-02      4.115390e+00     4.734922e+000    -3.937392e+001
     2.0e-02      4.121552e+00     4.735492e+000    -3.927266e+001
 
     ***********************************************************
     *                                                         *
     * Values at the left hand end (x=-0.6) by solving eq2     *
     *                                                         *
     ***********************************************************
        y0          x(-0.6)           y(-0.6)         dy/dx(-0.6)
     1.0e-02     -6.000000e-01     3.772773e+090    -3.001057e+092
     1.1e-02     -6.000000e-01     4.150051e+090    -3.301162e+092
     1.2e-02     -6.000000e-01     4.527328e+090    -3.601268e+092
     1.3e-02     -6.000000e-01     4.904606e+090    -3.901374e+092
     1.4e-02     -6.000000e-01     5.281883e+090    -4.201480e+092
     1.5e-02     -6.000000e-01     5.659161e+090    -4.501586e+092
     1.6e-02     -6.000000e-01     6.036438e+090    -4.801691e+092
     1.7e-02     -6.000000e-01     6.413716e+090    -5.101797e+092
     1.8e-02     -6.000000e-01     6.790993e+090    -5.401903e+092
     1.9e-02     -6.000000e-01     7.168270e+090    -5.702009e+092
     2.0e-02     -6.000000e-01     7.545548e+090    -6.002114e+092

  
> 

#
# Just for illustration, plot the results of solving
# eq1 and eq2 in the range where no singularity occurs
#
# No-one could claim that these "agree"
#
  plots:-odeplot( s1(0.322, 2.05e-05, 0.01, 5, -0.6),
                  [x, y(x)],
                  x=4.05..5
                );
  plots:-odeplot( s2(0.322, 2.05e-05, 0.01, 5, -0.6),
                  [x, y(x)],
                  x=4.05..5
                );

 

 
> 

 

Download tryCat2.mw

Well...

tomleslie 13876
February 05 2020
1 2

You don't provide enough information for anyone here to "solve" your problem, so only very gerneral guidelines can be given.

You can use the `try-catch` construct.  If the 'try' statement sequence succeeds - no problem. If the 'try' statement sequence produces an ERROR then the 'catch' statement sequence is executed - and can do something completely different

See the attached "toy" example, which "catches" a potential divide-by-zero error, and does something  different and (vaguely?) useful when the error occurs

> 

  restart;
  res1:=Array(-5..5):
  res2:=Array(-5..5):
  f := proc(x)
            if   x=3
            then error "Naughty - you tried to divide by 0 when x= %1", x
            else 1/(x-3)
            end if
       end proc:
> 

#
# This code will fail and execution will cease
# when the error condition is met
#
# Note that the result array is not populated
# beyond the point where the error occurs
#
  for i from -5 by 1 to 5 do
      res1[i]:=f(i);
  end do:
  res1;

Error, (in f) Naughty - you tried to divide by 0 when x= 3

  

Array(%id = 18446744074327883774)

 (1)
> 

#
# You can access the error message above through the
# 'lastexception' global
#
  lastexception;

f, "Naughty - you tried to divide by 0 when x= %1", 3

 (2)
> 

#
# This code will also "fail" when the error condition
# is met, but the try-catch construct allows
# execution to continue (and do something different)
# when the error condition occurs
#
# Note that the result array is fully populated
#
  for i from -5 by 1 to 5 do
      try res2[i]:=f(i)
      catch "Naughty":
             printf( "\n     Whoops! Got an error" );
             printf( " when i=%d\n", [lastexception][-1]);
             res2[ [lastexception][-1] ]:=infinity;
      end try;
  od:
  res2;


     Whoops! Got an error when i=3

  

Array(%id = 18446744074327883894)

 (3)
> 

 

> 

 


 

Download tryCat.mw

You might want to read the help page...

tomleslie 13876
February 04 2020
3 1

to understand the difference between the commands

StringTools:-Split(); and
StringTools:-StringSplit();

whose actions are illustrated in the attached

> 

restart;
StringTools['Split']("fdsg543656fgh:576fghs:dsfg::657ufdgdsg", "::");
StringTools['Split']("fdsg543656fgh:576fghs:dsfg::657ufdgdsg", ":");
StringTools['StringSplit']("fdsg543656fgh:576fghs:dsfg::657ufdgdsg", "::");
StringTools['StringSplit']("fdsg543656fgh:576fghs:dsfg::657ufdgdsg", ":");

["fdsg543656fgh", "576fghs", "dsfg", "", "657ufdgdsg"]

  

["fdsg543656fgh", "576fghs", "dsfg", "", "657ufdgdsg"]

  

["fdsg543656fgh:576fghs:dsfg", "657ufdgdsg"]

  

["fdsg543656fgh", "576fghs", "dsfg", "", "657ufdgdsg"]

 (1)
> 

 

Download strSpl.mw

Is this what you expect...

tomleslie 13876
February 01 2020
1 1

because I'm not very familiar with the regularChains package either

> 

  restart:
  kernelopts(version);
  with(RegularChains):
 

`Maple 2019.2, X86 64 WINDOWS, Nov 26 2019, Build ID 1435526`

 (1)
> 

  R := PolynomialRing([x,a,b,c]):
  sys := [a*x^2+b*x+c]:
  N := []:
  P := [a]:
  H := [x]:
  dec := RealTriangularize(sys,N,P,H,R):
  Display(dec, R);

"[{[[a x^2+b x+c=0,],[-4 a c+b^2>0 and a>0 and c<>0,]],{[[a x+b=0,],[c=0,],[a>0 and b<>0,]],{[[2 a x+b=0,],[4 c a-b^2=0,],[c>0 and b<>0,]]]"

 (2)
> 

 

Download regch.mw

A further observation...

tomleslie 13876
January 31 2020
0 1

together with a somewhat simpler method of computing the maximum error.

Examine closely the vertical scales on the three error plots

> 

restart;
> 

HQ:= {diff(u(x, t), t) = 0.5*VectorCalculus:-Laplacian(u(x, t), 'cartesian'[x])};
BCs:= {u(0, t) = 0, u(Pi, t) = 0, u(x, 0) = 5*sin(3*x)};

{diff(u(x, t), t) = .5*(diff(diff(u(x, t), x), x))}

  

{u(0, t) = 0, u(Pi, t) = 0, u(x, 0) = 5*sin(3*x)}

 (1)
> 

#
# Analytic solution
#
  solA:= pdsolve([HQ[], BCs[]]);
#
# Numeric solution: default space- and time-step
#
  solN1:= pdsolve(HQ, BCs, numeric):
#
# Numeric solution: OP's values for  space- and time-step
#
  solN2:= pdsolve(HQ, BCs, numeric, timestep = 1/10000, spacestep = Pi/100):
#
# Numeric solution: yet another combination of space- and time-step
#
  solN3:= pdsolve(HQ, BCs, numeric, timestep = 1/1000, spacestep = 1/1000):

u(x, t) = 5*sin(3*x)*exp(-(9/2)*t)

 (2)
> 

#
# Plots of the discrepancy between the exact (ie analytic)
# solutions and the three numeric solutions obtained above
# with differetn space- and time-steps
#
  p1:=solN1:-plot3d(u(x, t) - rhs(solA), x = 0 .. Pi, t = 0 .. 2);
  p2:=solN2:-plot3d(u(x, t) - rhs(solA), x = 0 .. Pi, t = 0 .. 2);
  p3:=solN3:-plot3d(u(x, t) - rhs(solA), x = 0 .. Pi, t = 0 .. 2);

 

 

 
> 

####################################################
# Maximum error between the analytic solution and
# each of the three numeric solutions obtained
# above
####################################################
#
# Error with step sizes on Maple defaults
#
  max(plottools:-getdata(p1)[3]);
#
# Error with OP's step sizes (very close to that
# quoted by CarlLove!)
#
  max(plottools:-getdata(p2)[3]);
#
# Error with step sizes I originally recommended
# Approximately 1000x smaller
#
  max(plottools:-getdata(p3)[3]);

HFloat(0.0980695025205307)

  

HFloat(0.0013492645986454832)

  

HFloat(1.7123907076754818e-6)

 (3)
> 

 

 

NULL

Download pdeErr4.mw

More of an observation really...

tomleslie 13876
January 31 2020
0 0

The discrepancy between the analytic and numeric solutions is very dependent on the values of the 'step' and 'timestep' options supplied to the pdsolve/numeric command.
The two plots in the attached repesent the "error". Check the differencces in the vertical scale!


 

restart

HQ := {diff(u(x, t), t) = VectorCalculus:-Laplacian(u(x, t), 'cartesian'[x])}

{diff(u(x, t), t) = diff(diff(u(x, t), x), x)}

 (1)

BCs := {u(0, t) = 0, u(Pi, t) = 0, u(x, 0) = 5*sin(3*x)}

{u(0, t) = 0, u(Pi, t) = 0, u(x, 0) = 5*sin(3*x)}

 (2)

solA := pdsolve([HQ[], BCs[]])

u(x, t) = 5*sin(3*x)*exp(-9*t)

 (3)

``

solN1 := pdsolve(HQ, BCs, numeric)

_m612552576

 (4)

solN2 := pdsolve(HQ, BCs, numeric, timestep = 1/1000, spacestep = 1/1000)

_m610125472

 (5)

solN1:-plot3d(u(x, t)-rhs(solA), x = 0 .. Pi, t = 0 .. 2)

  

``

``

solN2:-plot3d(u(x, t)-rhs(solA), x = 0 .. Pi, t = 0 .. 2)

  

``

``

``


 

Download pdeErr3.mw

Tricky one!...

tomleslie 13876
January 25 2020
0 0

I have tried the technique of mouse-selecting the desired 2-D output expression and using the context menu sequence

Copy Special -> Copy as MathML

then a simple CTRL-V (ie paste) in the target application. This sometimes works, and sometimes doesn't! Success or failure appears to depend on both the target application and the filename extension which I specify in the target application.

Since you do not specify your "target" application, I can't specify a command sequence which will definitely work with this approach.

On the other hand, if I add an execution group to your worksheet as shown in the attached, then it always(?) produces the output file testML.txt which looks OK to me (I'm not a mathML guru).

You willl have to change the target file name to something appropriate for your machine

> 

restart:
> 

n:= 2:
> 

X:= Vector(n, symbol= x);

Vector(2, {(1) = x[1], (2) = x[2]})

 (1)
> 

W:= Matrix(3,2, symbol= w);

Matrix(3, 2, {(1, 1) = w[1, 1], (1, 2) = w[1, 2], (2, 1) = w[2, 1], (2, 2) = w[2, 2], (3, 1) = w[3, 1], (3, 2) = w[3, 2]})

 (2)
> 

V:= W.X;

Vector(3, {(1) = w[1, 1]*x[1]+w[1, 2]*x[2], (2) = w[2, 1]*x[1]+w[2, 2]*x[2], (3) = w[3, 1]*x[1]+w[3, 2]*x[2]})

 (3)
> 

answer:=seq(diff~(v,W), v= V);

answer := Matrix(3, 2, {(1, 1) = x[1], (1, 2) = x[2], (2, 1) = 0, (2, 2) = 0, (3, 1) = 0, (3, 2) = 0}), Matrix(3, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = x[1], (2, 2) = x[2], (3, 1) = 0, (3, 2) = 0}), Matrix(3, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 0, (3, 1) = x[1], (3, 2) = x[2]})

 (4)
> 

fname:="C:/Users/TomLeslie/Desktop/testML.xml";
XMLTools:-PrintToFile(fname, MathML:-Export(X));
XMLTools:-PrintToFile(fname, MathML:-Export(W));
XMLTools:-PrintToFile(fname, MathML:-Export(V));
XMLTools:-PrintToFile(fname, MathML:-Export(answer[1]));
XMLTools:-PrintToFile(fname, MathML:-Export(answer[2]));
XMLTools:-PrintToFile(fname, MathML:-Export(answer[2]));
close(fname):

"C:/Users/TomLeslie/Desktop/testML.xml"

 (5)
> 

 

Download toMML.mw

One way...

tomleslie 13876
January 25 2020
1 0

is shown in the attached

> 

  restart;
  with(plots):
  animate( implicitplot,
           [ x = -3*abs(sin(y)),
             x = -3 .. 3,
             y = 0 .. t,
             thickness = 5,
             size = [150, 150],
             axes = none,
             color = red
           ],
           t=0..6
         );

 
> 

 

Download animimpl.mw

Yet another way...

tomleslie 13876
January 24 2020
1 2

using complexplot3d() is shown in the attached.

Easy to change the options for 'grid', 'transparency', 'thickness', 'colorscheme', etc etc

I haven't bothered redefining sqrt(-1) - so the attached uses the default 'I'.

> 

  restart;
  with(plots):
  with(plottools):
  Aplot:= (fn, p1)-> complexplot3d
                     ( fn,
                       z=-1-I..1+I,
                       if   p1='surface'
                       then style=surface,
                            transparency=0.8
                       else style=wireframe,
                            grid=[10,10],
                            thickness=4
                       fi,
                       colorscheme=[ "zcoloring",
                                     z->z
                                   ]
                     ):
  Bplot:= fn-> display
               ( [ Aplot(fn, 'surface'),
                   Aplot(fn, 'wireframe'),
                   reflect( Aplot(fn, 'surface'),
                            [ [0,0,0], [1,0,0], [0,1,0] ]
                          ),
                   reflect( Aplot(fn,'wireframe'),
                            [ [0,0,0], [1,0,0], [0,1,0] ]
                          )
                 ]
               ):
   Bplot( Re(sqrt(z)) );
   Bplot( sqrt(z) );

 

 
> 

 

> 

 

Download comploit.mw

Need one more BC...

tomleslie 13876
January 23 2020
1 0

otherwise the solution still contains an integration constant.

See the attached

> 

  restart;
  ode:=diff(y(x),x,x)-y(x)=-4*sin(x)^3+9*sin(x);
#
# Solve the ode, subject to the bc
#
  ans:=dsolve([ode, y(0)=y(2*Pi)]);
#
# "Simplify" this solution.
#
# Note that another boundary conditon is necessary
# in order to evaluate the integration constant _C1
#
  ans:= op(1, ans)=simplify(add(op([2,1..2], ans)))+add(op([2,3..4], ans));

diff(diff(y(x), x), x)-y(x) = -4*sin(x)^3+9*sin(x)

  

y(x) = -exp(-x)*_C1*(exp(2*Pi)-1)/(exp(-2*Pi)-exp(2*Pi))+exp(x)*_C1*(exp(-2*Pi)-1)/(exp(-2*Pi)-exp(2*Pi))-(1/10)*sin(3*x)-3*sin(x)

  

y(x) = (exp(2*x)+exp(2*Pi))*_C1*exp(-x)/(exp(2*Pi)+1)-(1/10)*sin(3*x)-3*sin(x)

 (1)
> 

#
# Explicit solutions can be obtained if one sets
# y(0)=y(2*Pi)=some_numeric_value - as in (for
# example)
#
  ans2:=dsolve([ode, y(0)=0, y(2*Pi)=0]);
  ans2:=dsolve([ode, y(0)=1, y(2*Pi)=1]);
  ans2:=dsolve([ode, y(0)=2, y(2*Pi)=2]);

y(x) = -(1/10)*sin(3*x)-3*sin(x)

  

y(x) = -(exp(2*Pi)-1)*exp(-x)/(exp(-2*Pi)-exp(2*Pi))+exp(x)*(exp(-2*Pi)-1)/(exp(-2*Pi)-exp(2*Pi))-(1/10)*sin(3*x)-3*sin(x)

  

y(x) = -2*(exp(2*Pi)-1)*exp(-x)/(exp(-2*Pi)-exp(2*Pi))+2*exp(x)*(exp(-2*Pi)-1)/(exp(-2*Pi)-exp(2*Pi))-(1/10)*sin(3*x)-3*sin(x)

 (2)
> 

 


 

Download odeSol.mw

Don't use plot3d!...

tomleslie 13876
January 22 2020
0 2

whhc is trying to plot a "surface". You are only getting a line because the way you are defining the "surface", it has zero width in one dimension. (Nothing depends on your plot variable 't')

To plot  curves in 3D, your best bet is probably to use spacecurve(). You can color these in a variety of ways, just by setting the HSV or RGB triple  ot be dependent on 'i' and 'j'.

One possibility is shown in the attached

> 

  restart:
  with(plots):
  display
  ( seq
    ( seq
      ( spacecurve
        ( [i/sin(u), u, j],
          u = 0.00001 .. 3/2,
          color=ColorTools:-Color("HSV", [i*j/100, 0.5, 0.5]),
          thickness=5,
          view = [0 .. 10, 0 .. Pi/2, 0 .. 10]
        ),
        i=1..10
      ),
      j=1..10
    )
  );

 
> 

 


 

Download hueplot.mw

Consider a different approach!...

tomleslie 13876
January 20 2020
1 0
  1. Your original setup of seven loops will execute 20000000 times
  2. The conditions for executing the solve() command ( ie a > c and igcd(a, b, c, d, t, m, n) = 1 and abs(b)+abs(d)-n <> 0 ) will be satisfied 853750 times (out of the possible 20000000)
  3. How long might it take for a single 'solve()' command to execute. Make the (rash?) assumption that it might be ~0.1 seconds, so that the total execution time for the 'solve()' operations will be 85375 seconds or about 23.7hours
  4. You have two choices
    1. Wait a long time
    2. Consider whether a different approach to the problem might be feasible
First 84 85 86 87 88 89 90 Last Page 86 of 207