tomleslie

5023 Reputation

15 Badges

9 years, 218 days

MaplePrimes Activity


These are answers submitted by tomleslie

that a 3D-plot would be a good idea - anything wrong with producing a 2-D plot for each differetn value of the variable 'tswitch' as in the attached?

NULL

Set up functions

 

restart; F := [k[a1]*C[T]*(R-x[1](t)-x[2](t))-k[d1]*x[1](t), k[a2]*C[T]*(R-x[1](t)-x[2](t))-k[d2]*x[2](t)]; H := alpha*(x[1](t)+x[2](t))
Pars := [k[a1] = 10^(-4), k[a2] = 2*10^(-2), k[a1] = 10^(-4), k[d1] = 10^(-5), k[d2] = 10^(-5), R = 1, k[m] = 10^(-4)]

Start plotting

 

with(plots); FnOut := [op(F), F[1]+F[2]]; plts := NULL; bcs := x[1](0) = 0, x[2](0) = 0, y(0) = 0; for tswitch from 100 by 10 to 200 do Sys := (eval(`~`[`=`]([diff(x[1](t), t), diff(x[2](t), t), diff(y(t), t)], `~`[piecewise](t < tswitch, subs(C[T] = 1, FnOut), subs(C[T] = 0, FnOut))), Pars))[]; p := dsolve([Sys, bcs], numeric); plts := plts, odeplot(p, [t, y(t)], 0 .. 300, size = [800, 800]) end do; display(plts)

[k[a1]*C[T]*(R-x[1](t)-x[2](t))-k[d1]*x[1](t), k[a2]*C[T]*(R-x[1](t)-x[2](t))-k[d2]*x[2](t), k[a1]*C[T]*(R-x[1](t)-x[2](t))-k[d1]*x[1](t)+k[a2]*C[T]*(R-x[1](t)-x[2](t))-k[d2]*x[2](t)]

 

 

op(F)

k[a1]*C[T]*(R-x[1](t)-x[2](t))-k[d1]*x[1](t), k[a2]*C[T]*(R-x[1](t)-x[2](t))-k[d2]*x[2](t)

(2.1)

NULL


 

Download pwODE.mw

to Rouben's solution, which incorporate

  1. Using the Optimization package to obtain the desired diffusion coefficient: produces the same value as Rouben, ie 0.138865421347717
  2. Separating the setup and solving of the PDE from the routine used to fit it. This makes it easier to reuse the same code to generate the complete solution when the desired diffusion coefficient is found from (1) above
  3. From (2) above, one can calculate the residual errors for each of the values in the vector t_number, as well as producing plots etc, using the solution with the optimal diffusion coefficient found in (1) above

See the attached (based on Rouben's code)

  restart;

#
# Parameters and procedures
#
  t_number:=<0, 10, 20, 30, 40>:
  m_number:=<11.50000000, 4.641182547, 1.273311993, 0.361845238, 0.288711649>:
  L:=2:
  getSol:= proc(d)
                local Mx0:= 0.05, cx0:= Mx0/(1-Mx0),
                      Mdb_i:= m_number[1], ct0:= Mdb_i,
                      PDE:= diff(C(x,t),t)=d*diff(C(x,t),x,x),
                      IBC:= {C(x,0)=ct0, C(0,t)=cx0, D[1](C)(L,t)=0};
                return pdsolve(PDE, IBC, numeric);
          end proc :
  Q:= proc(d)
           local pds, i, S:=0:
           if   not type(d, numeric)
           then return 'procname(_passed)'
           end if;
           pds:=getSol(d);
           for i from 1 to 5 do
               pds:-value( t=t_number[i],
                           output=listprocedure
                         );                  # get the solution at the desired time
               eval(C(x,t), %);              # extract the C(x,t) at that time
               int(%, 0..2, numeric) / L;    # compute the average of C(x,t) at that time
               S += (% - m_number[i])^2;     # accumulate the residuals
           end do;
           return sqrt(S);
      end proc:

#
# Rouben's calculation
#
# A few representative values - shows that the
# "optimum" value for the parameter 'd' is
# somewhere around 0.14
#
  Q(0.1), Q(0.2), Q(0.5), Q(1.0);
  plot(Q(d), d=0.1 .. 0.3, numpoints=5, view=0..2);

HFloat(2.0125989798253805), HFloat(1.9268252656215228), HFloat(4.344399869795477), HFloat(4.74473396546933)

 

 

#
# Use the Optimization() package to get the value
# of diffusion coefficient which produces the
# minimum overall residual. Compute the solution
# at this "optimal" diffusion coefficient
#
  oSol:= Optimization:-Minimize(Q, 0..1);
  pdsOpt:= getSol( oSol[2][1] ):
#
# Using the optimal diffusion coefficient compute
# the average residual at each of the time-values
# in the vector t_number. These values *ought* to
# match somewhat with the corresponding value in
# the vector m_number
#
  resid:= Vector( 5,
                  i-> int
                      ( eval
                        ( C(x,t),
                          pdsOpt:-value
                                  ( t=t_number[i],
                                    output=listprocedure
                                  )
                        ),
                        0..2,
                        numeric
                      )/L
                 );
#
# Get the mean square error in the average C(x,t)-value
# at each of the time-values in the vector t_number
#
  `^`~((resid-m_number),2);

[.8821785835008983, Vector(1, {(1) = .13886542134771712})]

 

Vector(5, {(1) = 11.499999999999842, (2) = 4.009882865095136, (3) = 1.7332102428767426, (4) = .7664331411209234, (5) = .3558085433772424})

 

Vector[column](%id = 18446744074389251246)

(1)

#
# Plot C(x,t) for the optimal value of the diffusion
# coefficient for each of the time-values in t_number
#
  cols:= [red, green, blue, brown, black]:
  plots:-display( [ seq
                    ( pdsOpt:-plot
                      ( t=t_number[j],
                        x=0..L,
                        color=cols[j]
                      ),
                      j=1..numelems(t_number)
                    )
                  ]
                );
#
# Plot the "general 3-D" solution
#
  pdsOpt:-plot3d( x=0..L, t=0..t_number[-1]);

 

 

 

Download diffCoeff.mw

 

 

 

In this definition, what is s(t)?

Some totally unknown, completely random function perhaps??

Because until it is defined in some way, this system cannot be solved!

While you are thinking about this, you might want to think about the parameter 's0' - now that wouldn't be s(0), would it? Because if you define s(t) using an additional differential equation, then you will need another boundary condition. (Obviously, if it is a straightforward function definition, then you won't, and 's0' will be a 'simple' parameter)

requires that the variable be 'L' be given a value.

The attached uses L=10, althogh it can easily be changed

  restart;
#
# Have to give 'L' a value to enable
# numeric solution
#
  L:=10:
  PDE := 35139.16048*(diff(w(x, t), x, x, x, x))+98646.00937*(diff(w(x, t), t, t))-2771.636*(diff(w(x, t), x, x)) = 24883.20000:
  bcs:= w(0,t)=0, D[1](w)(0,t)=0, D[1,1](w)(L,t)=0, D[1,1,1](w)(L,t)=0:
  ics:= w(x,0)=0, D[2](w)(x,0)=0:
  sol:=pdsolve(PDE, [bcs, ics], numeric):
#
# Plot w(x,t) for x=L/2 and t=0..10
#
  sol:-plot(x=L/2, t=0..10);
#
# Plot w(x,t) for x=0..L and t=5
#
  sol:-plot(t=0.5, x=0..L);
#
# 3D plot of w(x,t) for x=0..L, and t=0..10)
#
  sol:-plot3d( x=0..L, t=0..10);

 

 

 

 

 

 


 

Download solvePDE2.mw

You cannot really expect the drawing tool within Maple to have the capabilities of a serious drawing program. Maple is a system for mathemetical calculation - not drawing.

You have two choices

  1. Use an external drawing packing (Gimp/Visio/whatever). Export from that package into a format understood by the Maple drawing tool, and then use an "import" command
  2. Use the plottoools() package within Maple to produce a suitable drawing, export the result in an appropriate format. Then use an import command in the target worksheet as in (1) above

upload actual worksheet using the big green up-arrow in the Mapleprimes toolbar.

The code you present in your oriignal post cannot produce the output you present for three reasons

  1. the variable 'q' is undefined
  2. you have two first order differential equations, which obviously require two initlal conditions in order for a numeric solution to be feasible. No initial conditions are supplied

What is wrong with a simple function definition, as in the attached; the function may be called from top-level or from within a package. Scoping is pretty much the same as any other function definition

  restart;
  r:= (x,y)-> Record( 'name'=x, 'age'=y):
  myModule:= module()
                    global r;
                    option package;
                    export makeRec;
                    makeRec:= proc( a, b )
                                    return r(a,b)
                              end proc;
             end module:
  with(myModule):

  s1:=r("JohnDoe1", 26):
  s2:=r("JohnDoe2", 28):
  s3:=makeRec("JohnDoe3", 30):
  s1:-name;
  s2:-age;
  s3:-age;

"JohnDoe1"

 

28

 

30

(1)

 

 


 

Download makeRec.mw

With R=0.5, the integrand is positive over some of the domain z=-R..R, y=0..R and negative over the rest. My "calibrated eyeball" thnks that, overall, the integral ought to be negative - which it is. See the attached

No idea what your plot command is doing because 'phi' is undefined


 

Integrand := parse(FileTools[Text][ReadFile]("C:/Users/TomLeslie/Desktop/1.txt")):
#
# Evaluate the integral
#
  R:=0.5:
  evalf(Int(Integrand, [z = -R .. R, y = 0 .. R]));
#
# Plot the integrand over the region specified
# by the limits of integration
#
  plot3d(Integrand, z=-R..R, y=0..R);

-0.1342931338e-2

 

 

 


 

Download posNegInt.mw

 

which doesn't use the 'background' option. (and just because I like alternatives)

Rendering is better in the worksheet than it is on this site - honest!

Differential Geometry

restart;

with(VectorCalculus):with(LinearAlgebra):with(plots):

Digits:=10:
interface(displayPrecision=10):
with(VectorCalculus):BasisFormat(false):

Frenet Frame

 

frenet := proc(p)
               local vel, speed, unitT, curv, kappa, unitN, unitB, tors:
               vel := diff(p,s):
               speed := simplify(Norm(vel)):
               unitT := simplify(vel/speed):
               curv := simplify(diff(unitT,s)):
               kappa := simplify(Norm(curv)):
               unitN := simplify(curv/kappa):
               unitB := simplify(unitT &x unitN):
               tors := diff(unitB,s)/~unitN:
              # print("Velocity",vel,"Speed",speed,"T",unitT,"Curvature",curv,"Kappa",kappa,"N",unitN,"B",unitB,"Torsion",tors);
               return([vel,curv,unitB]):
          end proc:

Testing

 

p :=unapply( <3*arcsinh(s/5), (25+s^2)^(1/2), 4*arcsinh(s/5)>, s):
ptan:=unapply(frenet(p(s))[1],s):
display( [ spacecurve
           ( p(s),
             s = 0 .. 16*Pi,
             numpoints = 1000
           ),
           display
           ( [ seq
               ( arrow
                 ( p(A),
                   ptan(A),
                   width=0.3,
                   length=4
                 ),
                 A=0..50
               )
             ],
             insequence=true
           )
         ]
       );

 

 

Download anim.mw

the same way as you would generate any other subGraph - see the attached

  restart;
  with(GraphTheory):
  with(RandomGraphs):
#
# Create some random weighted graph
# and draw it
#
  G := RandomDigraph(10, .5, weights=1..5);
  DrawGraph(G);

GRAPHLN(directed, weighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], Array(%id = 18446744074329129078), `GRAPHLN/table/1`, Matrix(%id = 18446744074329130638))

 

 

#
# Output the edges of the above graph then
# select a few of these
#
  E:=Edges(G);
  e:=E[1..5];
#
# Create a subgraph on the restricted set of edges
# then draw it
#
  s:=Subgraph(G, e);
  DrawGraph(s);
  

{[1, 2], [1, 6], [1, 9], [1, 10], [2, 3], [2, 5], [2, 6], [2, 7], [2, 8], [2, 9], [2, 10], [3, 2], [3, 6], [3, 9], [3, 10], [4, 1], [4, 2], [4, 5], [4, 7], [4, 8], [5, 4], [5, 7], [5, 8], [6, 1], [6, 2], [6, 3], [6, 4], [6, 8], [6, 9], [6, 10], [7, 1], [7, 2], [7, 4], [7, 6], [7, 8], [7, 10], [8, 1], [8, 2], [8, 3], [8, 4], [8, 5], [8, 6], [8, 9], [9, 2], [9, 5], [9, 10], [10, 3], [10, 4], [10, 7], [10, 8]}

 

{[1, 2], [1, 6], [1, 9], [1, 10], [2, 3]}

 

GRAPHLN(directed, weighted, [1, 2, 3, 6, 9, 10], Array(%id = 18446744074347199174), `GRAPHLN/table/5`, Matrix(%id = 18446744074347200374))

 

 

 


 

Download subG.mw

so I'm going to "fix" all of your typos by pure guesswork - in which case I get the attached

restart;
w:=(x)->sin(lambda*x)+b*cos(lambda*x)-sinh(lambda*x)-b*cosh(lambda*x);
b := -(sin(lambda*L)+sinh(lambda*L))/(cos(lambda*L)+cosh(lambda*L));

proc (x) options operator, arrow; sin(lambda*x)+b*cos(lambda*x)-sinh(lambda*x)-b*cosh(lambda*x) end proc

 

-(sin(lambda*L)+sinh(lambda*L))/(cos(lambda*L)+cosh(lambda*L))

(1)

expr:=eval(int(w(x)^2,x=0..L), L=10);
plot( expr, lambda=-1..1)

(1/8)*(20*lambda*exp(40*lambda)+80*sin(10*lambda)*lambda*exp(30*lambda)-80*cos(10*lambda)^2*lambda*exp(20*lambda)-6*sin(10*lambda)*cos(10*lambda)*exp(40*lambda)+6*exp(40*lambda)*cos(10*lambda)^2+40*lambda*exp(20*lambda)-80*sin(10*lambda)*lambda*exp(10*lambda)+exp(40*lambda)-10*sin(10*lambda)*exp(30*lambda)+14*cos(10*lambda)*exp(30*lambda)-12*sin(10*lambda)*cos(10*lambda)*exp(20*lambda)-4*exp(20*lambda)*sin(10*lambda)*cosh(10*lambda)-4*exp(20*lambda)*sinh(10*lambda)*cos(10*lambda)-4*exp(20*lambda)*sinh(10*lambda)*cosh(10*lambda)+20*lambda-10*sin(10*lambda)*exp(10*lambda)-14*cos(10*lambda)*exp(10*lambda)-6*sin(10*lambda)*cos(10*lambda)-6*cos(10*lambda)^2-1)*exp(-20*lambda)/(lambda*(cos(10*lambda)^2+2*cos(10*lambda)*cosh(10*lambda)+cosh(10*lambda)^2))

 

 

fsolve( expr=1, lambda=-0.5..0);
fsolve( expr=1, lambda= 0..0.5);

-.1007127230

 

.1007127230

(2)

 


 

Download solveEq2.mw

with a large numeber of points (eg 1000) you are going to run out of labels in the dataframe quite quickly - in fact after about 26 entries. Generating appropriate labels for your 'dataframe' will probably become more cumbersome than generating the data - so why bother??? For your requirement a 2X1000(ish) matrix would suffice.

If your datapoints have integer values, then Maple will calculate exact distances which will (almost certainly) contain square root symbols. Trust me, Excel isn't going to like these! It is probably a good idea to "force" a floating point conversion within the ExcelTools:-Export() command.

The attached illustrates the technique

  restart;
  with(Student:-Calculus1):
#
# Specify the "reference" points
#
  refpts:=[ [1,1],
            [2,4]
          ];
#
# Generate 'n' random points whose x-coordinates
# are in the range -xRange..xRange and y-coordinates
# are in the range -yRange..yRange.
#
# Specific values chose below are just for illustration
#
  n:=100: xRange:=10: yRange:=20:
  xR:=rand(-xRange..xRange):
  yR:=rand(-yRange..yRange):
  pts:=[seq([xR(), yR()],j=1..n)]:
#
# Produce a matrix
#
#  1. first row is the distance from refpts[1] to
#     each of the random points generated above
#  2. second row is the distance from refpts[2] to
#     each of the random points generated above
#
  d:= Matrix
      ( 2,
        n,
        (i,j)-> Distance
                ( refpts[i],
                  pts[j]
                )
      );
#
# Export the distance matrix to Excel. Note that if
# integers have been used in the above calculations
# then the computed distance value is likely to contain
# a square root symbol, whihc Excel isn't going to like!
# Best to force a floating point evaluation at this point
#
# OP will also have to change the file path to something
# appropriate for his/her machine
#
  ExcelTools:-Export( evalf~(d),
                      "C:/Users/TomLeslie/Desktop/test.xlsx",
                       1, "B2"
                    );

refpts := [[1, 1], [2, 4]]

 

_rtable[18446744074399240054]

(1)

 


Download toXL.mw

On my machine this produces the file

test.xlsx

as in the attached

restart;
expr1:=(x^2-1)*y/x^2+(x^2-1)/x^2;
thaw(algsubs((y+1)=freeze((y+1)), expr1));

(x^2-1)*y/x^2+(x^2-1)/x^2

 

(y+1)*(x^2-1)/x^2

(1)

 

Download frth2.mw

Although, I do admit that working out exactly which sub-expression to "freeze", can sometimes be a bit of a "black art"

  1. Please upload a worksheet using the big green up-arrow in the MaplePrimes toolbar
  2. By using a lot of guesswork, I managed to come up with the attached - but this did involve butchery on the code you originally presented
  3. Main syntax-type changes I made were
    1. Using square brackets ( ie [] not () ) when one wants to access an indexaxable quantity
    2. Realising that when one uses Statistics:-Mean() to compute the "mean" of a Matrix, the result is returned as a vector. Hence the quantity 'MeanY' in your worksheet  can only be used with a single index - so MeanY[1] rather than MeanY[1,1] everywhere it is used

For what it is worth, and bearing in mind the outright butchery I performed on the code you posted, the attached now actually runs without error. Whether or not it bears any relation to what you want, I have no idea


 

restart:               

interface(rtablesize = 50); with(LinearAlgebra): with(plots):with(Statistics): with(ArrayTools):
Y := Matrix(12, 1, [2, 3, 2, 7, 6, 8, 10, 7, 8, 12, 11, 14]);
MeanY := Mean(Y);
n := ArrayNumElems(Y);
type(Y, Matrix);
X := Matrix(12, 3, [1, 0, 2, 1, 2, 6, 1, 2, 7, 1, 2, 5, 1, 4, 9, 1, 4, 8, 1, 4, 7, 1, 6, 10, 1, 6, 11, 1, 6, 9, 1, 8, 15, 1, 8, 13]);
 k := NumElems(X)/n;
Xprimed := Transpose(X);
XprimedX := Xprimed . X;
XprimedXinverse := MatrixInverse(XprimedX);
XprimedY := Xprimed . Y;
Betahat := XprimedXinverse . XprimedY;
Betahat := convert(Betahat, float);
fit := plot3d(5.38+3.01*X1-1.29*X2, X1 = 0 .. 8, X2 = 2 .. 15);
data := pointplot3d({[0, 2, 2], [2, 5, 7], [2, 6, 3], [2, 7, 2], [4, 7, 10], [4, 8, 8], [4, 9, 6], [6, 9, 12], [6, 10, 7], [6, 11, 8], [8, 13, 14], [8, 15, 11]}, axes = normal, color = red, symbol = soliddiamond, symbolsize = 40, scaling = unconstrained, title = `\` Data vs Best Fit ' `);

display({data, fit});
Yprimed := Transpose(Y);
 YprimedY := Yprimed . Y;
 BetahatPrime := Transpose(Betahat);
SSE := YprimedY-BetahatPrime . XprimedY;
SampleVariance := (YprimedY-BetahatPrime . XprimedY)/(n-k-1);
SSR := Transpose(Betahat) . Transpose(X) . Y-n*MeanY(1, 1)^2;

SST := SSR+SSE;
 F := SSR[1, 1]*(n-k-1)/(k*SSE);

RsquaredNonCentered := Determinant(SSR)/Determinant(SST);
 R := sqrt(%);
RsquaredCentered:= ((Transpose(Betahat) . Transpose(X) . Y) - n*MeanY[1]^2) .(1/((Yprimed . Y) - n*MeanY[1]^2));
F := SSR[1, 1]*(n-k-1)/(k*SSE[1, 1]);
 SSE; SSR; SSR+SSE;
 Total := Yprimed . Y-n*MeanY[1]^2;

[10, 10]

 

Vector(12, {(1) = 2, (2) = 3, (3) = 2, (4) = 7, (5) = 6, (6) = 8, (7) = 10, (8) = 7, (9) = 8, (10) = 12, (11) = 11, (12) = 14})

 

Vector(1, {(1) = 7.5})

 

n := 12

 

true

 

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

 

k := 3

 

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

 

Matrix(3, 3, {(1, 1) = 12, (1, 2) = 52, (1, 3) = 102, (2, 1) = 52, (2, 2) = 296, (2, 3) = 536, (3, 1) = 102, (3, 2) = 536, (3, 3) = 1004})

 

Matrix(3, 3, {(1, 1) = 309/317, (1, 2) = 77/317, (1, 3) = -145/634, (2, 1) = 77/317, (2, 2) = 411/2536, (2, 3) = -141/1268, (3, 1) = -145/634, (3, 2) = -141/1268, (3, 3) = 53/634})

 

Vector(3, {(1) = 90, (2) = 482, (3) = 872})

 

Vector(3, {(1) = 1704/317, (2) = 3819/1268, (3) = -815/634})

 

Vector(3, {(1) = 5.375394322, (2) = 3.011829653, (3) = -1.285488959})

 

 

 

 

Vector[row](12, {(1) = 2, (2) = 3, (3) = 2, (4) = 7, (5) = 6, (6) = 8, (7) = 10, (8) = 7, (9) = 8, (10) = 12, (11) = 11, (12) = 14})

 

Vector(1, {(1) = 840})

 

Vector[row](3, {(1) = 5.375394322, (2) = 3.011829653, (3) = -1.285488959})

 

Vector(1, {(1) = 25.458990522000022})

 

Vector(1, {(1) = 3.1823738152500027})

 

Vector(1, {(1) = 139.54100947799998})

 

Vector(1, {(1) = 165.0})

 

Vector(1, {(1) = 14.616029582415962})

 

RsquaredNonCentered := .845703087745454

 

R := .919621165342259

 

Vector(1, {(1) = .8457030877454544})

 

F := 14.6160295824160

 

Vector(1, {(1) = 25.458990522000022})

 

Vector(1, {(1) = 139.54100947799998})

 

Vector(1, {(1) = 165.0})

 

Matrix(%id = 18446744074376644238)

(1)

 


 

Download matProb.mw

 

Just put the desired values of Mh in a list and add another loop, which will permit a sequence of plots to be obtained. Thes can then be displayed together. Just make sure that the number of entries in the list 'colors' is the same as the number of entries in the list 'Mh'

Code changes are highlighted in the attached

restart:
with(plots):
m:=10:H:=1:b:=0.02:a:=0.05:V:=Array(0..m): V[0]:=1-exp(-t):
plts:=NULL:
Mh:=[1,2,3,4]:
colors:=[red, blue, green, black]:
for zz from 1 to 4 by 1 do
    for k from 1 to m do
        if   k=1
        then chi:=0:
        else chi:=1: # OP's code here didn't make sense
        fi:
        p:=0:
        for j from 0 to k-1 do
            p:=p+(V[k-1-j]*diff(V[j],t$2)-diff(V[k-1-j],t)*diff(V[j],t)-a*(2*diff(V[k-1-j],t)*diff(V[j],t$3)-diff(V[k-1-j],t$2)*diff(V[j],t$2)-V[k-1-j]*diff(V[j],t$4))):
        od:
        p:=(p+diff(V[k-1],t$3)-b*(diff(V[k-1],t$2)+t*diff(V[k-1],t$3))-Mh[zz]*diff(V[k-1],t))*h*H:
        p:=factor(p):
        V[k]:=(-int(p,t)+0.5*exp(t)*int(exp(-t)*p,t)+0.5*exp(-t)*int(exp(t)*p,t)+chi*V[k-1]+C1+C3*exp(-t)):
        v:=unapply(V[k],t):
        V[k]:=frontend(expand,[V[k]]):
        V[k]:=subs(C3=solve(eval(subs(t=0,diff(V[k],t))),C3),V[k]):
        V[k]:=frontend(expand,[V[k]]):
        V[k]:=subs(C1=solve(eval(subs(t=0,-V[k]-diff(V[k],t))),C1),V[k]):
    od:
    appr:=0:
    for k from 0 to m do
        appr:=appr+V[k]:
    od:
    u_appr:=unapply(appr,(h,t)):
    u_appr_1:=unapply(diff(u_appr(h,t),t),(h,t)):
    evalf(u_appr_1(-0.4,t)):
    plts:=plts, plot([u_appr_1(-0.4,t)],t=0..4,0..1.2,color=colors[zz],axes=frame):
od:
display([plts]);

 

 

 


 

Download odeProb.mw

2 3 4 5 6 7 8 Last Page 4 of 113