tomleslie

Diififult to work out...

Answered: tomleslie 13876

July 29 2020

0 0

exactly what you are trying to achieve

If I make a random gues, that you want to compute F[0..10], Theta[0..10] and G[0..10], then the attached works

>

restart:

>	F[0]:=0: F[1]:=1: F[2]:=a: Theta[0]:=b: Theta[1]:=-1: G[0]:=beta: G[1]:=K[s]*beta: delta(-1):=0: delta(0):=1: for k from 1 to 10 do delta(k):=0: od: A:=1: B:=1: C:=1:

>

  for k from 0 by 1 to 8 do;
      F[k+3]:= eval
               ( solve
                 ( (k+1)*(k+2)*(k+3)*F[k+3]
                   +2*nu*add
                         ( delta(k-m-1)*(m+1)*(m+2)*(m+3)*F[m+3],
                           m = 0 .. k
                         )
                   +2*nu*(k+1)*(k+2)*F[k+2]
                   +A*add
                      ( F[k-m]*(m+1)*(m+2)*F[m+2],
                        m = 0 .. k
                      )
                   -A*add( (k-m+1)*(m+1)*F[k-m+1]*F[m+1],
                           m = 0 .. k
                         )
                   -B*M*(k+1)*F[k+1],
                   F[k+3]
                 ),
                 [ Pr=6.2, M=2, nu=0.3, phi=0.05, lambda=1, Sc=1, Ks=1]
               );
      Theta[k+2]:= eval
                   ( solve
                     ( C*(k+1)*(k+2)*Theta[k+2]
                       +2*nu*C*add
                               ( delta(k-m-1)*(m+1)*(m+2)*Theta[m+2],
                                 m = 0 .. k
                               )
                       +2*nu*C*(k+1)*Theta[k+1]
                       +Pr*add
                           ( F[k-m]*(m+1)*Theta[m+1],
                             m = 0 .. k
                           )
                       -Pr*add
                           ( Theta[k-m]*(m+1)*F[m+1],
                             m = 0 .. k
                           ),
                       Theta[k+2]
                     ),
                     [Pr=6.2, M=2, nu=0.3, phi=0.05, lambda=1, Sc=1, Ks=1]
                   );
      G[k+2]:= eval
               ( solve
                 ( (k+1)*(k+2)*G[k+2]
                   +2*nu*add
                         ( delta(k-m-1)*(m+1)*(m+2)*G[m+2],
                           m = 0 .. k
                         )
                   +2*nu*(k+1)*G[k+1]
                   +Sc*add
                       ( F[k-m]*(m+1)*G[m+1],
                         m = 0 .. k
                       )
                   -lambda*Sc*G[k]
                   +2*lambda*Sc*add
                                ( G[k-m]*G[m],
                                  m = 0 .. k
                                )
                   -lambda*Sc*add
                              ( add
                                ( G[t]*G[m]*G[k-m-t],
                                  m = 0 .. k-t
                                ),
                                t = 0 .. k
                              ),
                   G[k+2]
                 ),
                [Pr=6.2, M=2, nu=0.3, phi=0.05, lambda=1, Sc=1, Ks=1]);
  od:

>	eval(F); eval(Theta); eval(G);

(1)

>

Download DTM2.mw

No analytic solution???...

Answered: tomleslie 13876

July 27 2020

0 1

Althugh I couldn't prove it either way, I seriously doubt that this syetem of ODEs can be solve "analytically". Which means that you are left with the option of a "numerical" solution.

In order to obtain a numerical solution, all parameters have to be assigned numerical values and initial condiotns for the system have to be specifiied

I have absolutely no idea what the numerical values of either the parameters or the initial conditions should be -so in the attached, I just gave all of these (essentially) arbitrary values - yU shoud adjust/modify as necessary

MAking all of these assumptions, I do actually get a soltuon. If I just thihnk in terms of "susceptible", "colonized" and "infwcted", it looks (vaguely) plausble.

>

restart:

>

#
# Define the ODEs for the system
#
  odeSys:= diff(S(t),t) = z*(1 - P[1] - P[2]) + w*S(t) - x*Q[1]*C(t)*S(t)/N - x*Q[1]*Iota(t)*S(t)/N - v*S(t),
           diff(C(t),t) = z*P[1] + x*Q[1]*C(t)*S(t)/N + x*Q[2]*Iota(t)*S(t)/N - w*C(t) - k*o*C(t) - v*C(t),
           diff(Iota(t),t) = C(t)*k*o - Iota(t)*v + z*P[2];

diff(S(t), t) = z*(1-P[1]-P[2])+w*S(t)-x*Q[1]*C(t)*S(t)/N-x*Q[1]*Iota(t)*S(t)/N-v*S(t), diff(C(t), t) = z*P[1]+x*Q[1]*C(t)*S(t)/N+x*Q[2]*Iota(t)*S(t)/N-w*C(t)-k*o*C(t)-v*C(t), diff(Iota(t), t) = k*o*C(t)-Iota(t)*v+z*P[2]

(1)

>

#
# I suspect that the above ODE system *cannot* be solved
# "analytically" but will have to be solved "numerically",
# in which case it will be necessary to define numerical
# values for all parameters in the system
#
# I have absolutely no idea what these parameters
# should be - so the following just sets every single
# one of them to 0.5 - just for convenience.
#
# If OP has any idea what these parameters actually
# are then they should be entered as
#
# pars:= [ N=valueForN,
#          k=valueFork,
#          ....
#          and so on
#        ]~
#
  pars:= [seq( j=0.5,j in`union`(indets~([odeSys], name)[]) minus {t})];

(2)

>	# # Define the initial conditions for the system # # Again I have no idea what these conditions should # be, so the following is entirely arbitrary # ICs:= S(0) = 100, C(0) = 1, Iota(0) = 1;

(3)

>	# # Solve the system numerically # soln:=dsolve([eval([odeSys], pars)[], ICs], numeric);

(4)

>	# # Plot these numerical solutions # plots:-odeplot( soln, [ [t, S(t)], [t, C(t)], [t, Iota(t)] ], t = 0 .. 5, color=[red, green, blue] );

>

Download odeSol.mw

That is ne "ugly" equation...

Answered: tomleslie 13876

July 26 2020

0 0

The attached worksheet show how to get a couple of answers which are *probably* correct. Even increasing Maple's Digits setting doesn't help much when dealing with this expression sincee you have used 10-digit floating point coefficients in many places. To get "better" (or more) answers, I think all the coefficients in the definition of the expression would need substantially increased precision (20-25 digits maybe) or prefereably be expressed as rational numbers (if at all possible)

>

(1)

>

(2)

>

(3)

>

(4)

>

(5)

>

plot(F(lambda), lambda = -.2 .. .2); plot(F(1000*lambda), lambda = -.2 .. .2, numpoints = 10000, view = [-.2*(1/1000) .. .2*(1/1000), -0.1e-1 .. 0.1e-1]); NewtonsMethod(F(lambda), lambda = 0); NewtonsMethod(F(lambda), lambda = 0.65e-1)

(6)

>

Download newt.mw

The answer is not 76 +2/25...

Answered: tomleslie 13876

July 25 2020

1 4

As Rouben has already pointed out , the "green" region in your diagram has an area of 952/25 - or if you prefer 38+2/25.

Note that the 'integer part of this number is exactly half of the value you *seem* to think is correct - this may help you find the error in your calculation

Whilst Rouben's answer is correct, and is essentially an application of "similar" triangles, you appear not to believe it. So try a different way of doing the same calculation using the geom3d() package as in the attached worksheet - which also demonstrates that the area of the rectangle SJIK is 1152/25, and the area of SJIK minus the area of ABCD is 952/25

>	restart: with(geom3d): # # The following will issue a warning because (by default) # the name 'I' is used by Maple to represent the square # root of -1. For the purposes of this worksheet it just # represents the name of a 'point' - so ignore the warnng! # local I:

Warning, The imaginary unit, I, has been renamed _I

>	# # Image supplied by the OP #

>

#
# Specify some of the points in the above image
#
  point(L, [ 0, 0, 12]):
  point(E, [-2, 0, 7]):
  point(H, [ 2, 0, 7]):
  point(F, [-2, 2, 7]):
  point(G, [ 2, 2, 7]):
#
# Define the lines through the points L,E,S and
# L, G, I
#
  line( LES, [L, E]):
  line( LGI, [L, G]):

#
# Define the z=0 plane
#
  plane( p, z=0, [x,y,z]):
#
# Intersection of the lines LES, LGI,
# with the plane z=0, which defines the points
# S and I
#
  intersection( S, LES, p):
  intersection( I, LGI, p):

>	# # Using the coordinates of the points 'S' and 'I' # compute the area of the rectangle SJIK # area1:=( xcoord(I)-xcoord(S) )*( ycoord(I)-ycoord(S)); # # Now subtract the area (=8) of the rectangle ABCD # to obtain the "green" area # area2:=area1-8

(1)

>

Download areaCalc.mw

try the geom3d package...

Answered: tomleslie 13876

July 24 2020

0 4

has several useful commands, such as

point() which will alow you to define a "point"

plane() which will allow you to define a "plane", and

distance() which will allow you to determine the "distance" between a variety of objects (such as a point and a plane)

which of these is causing a problem?

Use the option 'listprocedure'...

Answered: tomleslie 13876

July 24 2020

2 0

in the dsolve() command.

This will enable individual procedures for the dependent variables to be assigned to any names you want and used as simple function calls.

See the "toy" example in the attached

In future please use the big green up-arrow in the Mapleprimes toolbar to upload an executable worksheet which illustrates your problem. No-one is going to retype your example from a "picture" - which is why the "toy" example in the attached, just uses the first entry from the dsolve/numeric help page

>

restart;

>	dsys := {diff(x(t),t)=y(t),diff(y(t),t)=-x(t),x(0)=1,y(0)=0};

>	# # Basic method to get numeric results from ODE solution # dsn1 := dsolve(dsys,numeric); dsn1(2);

(1)

>

#
# Less basic method, which outputs individual procedures
# for the dependent variables, which can then be assigned
# and used.
#
  dsn2:= dsolve(dsys,numeric, output=listprocedure);
#
# Assign the procedures which compute x(t) and y(t) to the
# names xt and yt
#
  xt:=eval(x(t), dsn2);
  yt:=eval(y(t), dsn2);
#
# So the following will compute x(2) and y(2)
#
  xt(2);
  yt(2);

(2)

>	# # Or you can do this all in one command as # xtt,ytt:= eval( [x(t),y(t)], dsolve(dsys,numeric, output=listprocedure))[]; xtt(2);ytt(2);

(3)

>

Download solDE2.mw

Several typos...

Answered: tomleslie 13876

July 23 2020

3 2

mainly missing function arguments - like writing 'w' when you mean w(x,t).

Also one too many boundary conditons. v(x,t) is never differentiated wrt to x, so you cannot have a boundary condition of the form v(some_value_of_x,t)=some_function_of_t - so I just commented out the invalid bundary condition

Everything now seems to be working, although v(x,t) doesn't do anything interesting. See the attached

>

  restart;
  with(plots):
  pde1 := { diff(u(x, t), t) + Gamma*diff(u(x, t), x) = mu(w(x,t)- u(x,t)) - 0.5*beta*v(x,t)*(w(x,t) - u(x,t)) - beta*u(x,t)*v(x,t),
            diff(w(x, t), t) - Gamma*diff(w(x, t), x) = mu(u(x,t) - w(x,t)) - 0.5*beta*v(x,t)*(u(x,t) - w(x,t)) - beta*w(x,t)*v(x,t),
            diff(v(x,t),t)=beta*u(x, t)*v(x, t)-alpha*v(x, t)
          };

  IBC1 := { u(0, t) = 0,
            u(x, 0) = exp(-200*(x - 0.4)^2),
            w(1, t) = 0,
            w(x, 0) = exp(-200*(x - 0.6)^2),
            v(x,0)=0
            #v(0,t)=sin(t)
          };
  mu := 0.0009; Gamma := 10; alpha:=1; beta:=0.1;

  pds1 := pdsolve(pde1, IBC1, numeric, time = t, range = 0 .. 1, spacestep = 0.01);

{diff(u(x, t), t)+Gamma*(diff(u(x, t), x)) = mu(w(x, t)-u(x, t))-.5*beta*v(x, t)*(w(x, t)-u(x, t))-beta*u(x, t)*v(x, t), diff(w(x, t), t)-Gamma*(diff(w(x, t), x)) = mu(u(x, t)-w(x, t))-.5*beta*v(x, t)*(u(x, t)-w(x, t))-beta*w(x, t)*v(x, t), diff(v(x, t), t) = beta*u(x, t)*v(x, t)-alpha*v(x, t)}

${u(0, t) = 0, u(x, 0) = exp(-200*(x-.4)^2), v(x, 0) = 0, w(1, t) = 0, w(x, 0) = exp(-200*(x-.6)^2)}$

(1)

>	pds1:-plot3d( u(x,t), x=0..1, t=0..1); pds1:-plot3d( v(x,t), x=0..1, t=0..1); pds1:-plot3d( w(x,t), x=0..1, t=0..1);

>

Download pdesol.mw

Use the ImportMatrix() command...

Answered: tomleslie 13876

July 23 2020

0 0

with the source=MATLAB option, as in

M1:=ImportMatrix( "full_path_to_the file_you want"", source=MATLAB, output=matrices);

see the help at ?ImportMatrix

If you still have problems, upload the file containg the matrix here. This site won't accept files with a '.mat' extension, so you may have to zip the file first

Don't think...

Answered: tomleslie 13876

July 23 2020

1 2

an algebraic solution for this inequality is possible, allthough it is pretty trivial to generate a numeric one. See the attached

BTW, since your ineqality is expressed in terms of a single variable, namely 't', Im not sure what the the command plots:-inequal() is actually "doing"

>	restart: with(plots):

>	k:= 1; c:= 5; sigma:= .85; N=10; t0:= 0.0065; x:= t->4*exp(-t);

(1)

>	ineq:=sigmak/(2c+k)*abs(x(t))<abs(x(t)-x(t0));

(2)

>

#
# Plot both sides of the inequality
#
  plot( [rhs(ineq), lhs(ineq)],
        t=-1..1,
        color=[red, blue]
      );
#
# Shade the region where the inequality is
# satisfied
#
  shadebetween( rhs(ineq),
                lhs(ineq),
                t=-1..1,
                positiveonly,
                color=red,
                changefill=[color=blue]
              );

>	# # Generate the t-values which bound the region # where the inequality satisfied # # t < t1 or # t > t2 # t1, t2:= fsolve( rhs(ineq)=lhs(ineq), t=-1..0), fsolve( rhs(ineq)=lhs(ineq), t=0..1);

(3)

>

Download ineqs.mw

Hmmmm...

Answered: tomleslie 13876

July 22 2020

0 0

Well, it isn't too difficult to

solve for N solving ODE11 and ODE12 simultaneously
solve for t_2 solving ODE11 and ODE12 simultaneouly
find t__3 and t__4

but the condition

if (N<=t_4 and M>=t_3) then

is never satisfied, so it all seems a bit pointless. However at least this avoids dealing with the subsequent issue of

if (d1< d2)

whixh can never be evaluated since d2 is defined nowhere in your worksheet.

See the attached where I have generated all SOL1, t3 and t4 values as tables - just for illustration. This would probably be unnecessary if subsequent conditional clauses ever evaluated to anything meaningful

Down-stream credit and time-dependent demand rate,

(1)

Finite production rate proportional to the demand rate

(2)

3. Mathematical model for three echelon supply chain considering deteriorating items under quadratic demand

3.1. Manufacturer’s inventory model

The manufacturer’s inventory level is governed by the following differential equations:

(3)

(4)

$i__m1 := unapply(rhs(dsolve({ode1, i__m1(0) = 0})), t)$

proc (t) options operator, arrow; ((-c*t^2*theta__m^2+b*t*theta__m^2+2*c*t*theta__m-b*theta__m+theta__m^2-2*c)*exp(theta__m*t)*a*N^alpha*(lambda-1)/theta__m^3-(-b*theta__m+theta__m^2-2*c)*a*N^alpha*(lambda-1)/theta__m^3)*exp(-theta__m*t) end proc

(5)

$i__m2 := unapply(rhs(dsolve({ode2, i__m2(0) = 0})), t)$

proc (t) options operator, arrow; (-(-c*t^2*theta__m^2+b*t*theta__m^2+2*c*t*theta__m-b*theta__m+theta__m^2-2*c)*N^alpha*a*exp(theta__m*t)/theta__m^3+(-b*theta__m+theta__m^2-2*c)*N^alpha*a/theta__m^3)*exp(-theta__m*t) end proc

(6)

Taking value of from article we get:

t__1 := t__2*(2*b*t__2^2*`θ__m`-2*c*t__2^2+3*b*t__2+3*t__2*`θ__m`+6)/(6*(lambda-1))

t__2*(2*b*t__2^2*theta__m-2*c*t__2^2+3*b*t__2+3*t__2*theta__m+6)/(6*lambda-6)+t__2

(7)

Set-up cost

Holding cost

(8)

HC__m := (int(h__m(t)*i__m1(t), t = 0 .. t__1))*(int(h__m(t)*i__m2(t), t = 0 .. t__2))/T

Cost due to deterioration

DC__m := C__m*`θ__m`*(int(i__m1(t), t = 0 .. t__1)+int(i__m2(t), t = 0 .. t__2))/T

Opportunity interest loss in m- shipments by offering credit period M to the distributer

OIL__m := P__m*i__om*m*(int(R(N, t), t = 0 .. M))/T

Therefore, the total cost per unit time () for manufacturer is:

3.2. Distributor’s inventory model

T

he rate of change of inventory for distributor is governed by the following differential equation:

$i__d := unapply(rhs(dsolve({ode3, i__d(t__3) = 0})), t)$

proc (t) options operator, arrow; (-(-c*t^2*theta__d^2+b*t*theta__d^2+2*c*t*theta__d-b*theta__d+theta__d^2-2*c)*a*N^alpha*exp(theta__d*t)/theta__d^3+(-c*t__3^2*theta__d^2+b*t__3*theta__d^2+2*c*t__3*theta__d-b*theta__d+theta__d^2-2*c)*a*N^alpha*exp(theta__d*t__3)/theta__d^3)*exp(-theta__d*t) end proc

(9)

-(-b*theta__d+theta__d^2-2*c)*a*N^alpha/theta__d^3+(-c*t__3^2*theta__d^2+b*t__3*theta__d^2+2*c*t__3*theta__d-b*theta__d+theta__d^2-2*c)*a*N^alpha*exp(theta__d*t__3)/theta__d^3

(10)

Ordering cost

Holding cost

(11)

HC__d := m*(int(h__d(t)*i__d(t), t = 0 .. t__3))/T

Cost due to deterioration

DC__d := P__m*`θ__m`*m*(int(i__d(t), t = 0 .. t__3))/T

Opportunity interest loss () in mn shipments due to offer of credit period incurred by the retailer:

OIL__d := P__d*i__OD*m*n*(int(R(N, t), t = 0 .. N))/T

Case 1D: M ≤

the interest earned () is:

IE__d1 := P__d*i__e*m*(int(R(N, t)*t, t = 0 .. M))/T

the interest charged () is:

IC__d1 := P__m*i__c*m*(int(i__d(t), t = M .. t__3))/T

Case 2D: M >

IE__d2 := P__d*i__e*m*(int(R(N, t)*t, t = 0 .. t__3)+M-t__3+int(R(N, t), t = 0 .. t__3))/T

3.3. Retailer’s inventory model

the rate of change of inventory at any time t is governed by the following differential equation:

$i__r := unapply(rhs(dsolve({ode4, i__r(t__4) = 0})), t)$

proc (t) options operator, arrow; (-(-c*t^2*theta__r^2+b*t*theta__r^2+2*c*t*theta__r-b*theta__r+theta__r^2-2*c)*a*N^alpha*exp(theta__r*t)/theta__r^3+(-c*t__4^2*theta__r^2+b*t__4*theta__r^2+2*c*t__4*theta__r-b*theta__r+theta__r^2-2*c)*a*N^alpha*exp(theta__r*t__4)/theta__r^3)*exp(-theta__r*t) end proc

(12)

-(-b*theta__r+theta__r^2-2*c)*a*N^alpha/theta__r^3+(-c*t__4^2*theta__r^2+b*t__4*theta__r^2+2*c*t__4*theta__r-b*theta__r+theta__r^2-2*c)*a*N^alpha*exp(theta__r*t__4)/theta__r^3

(13)

Ordering cost () :

Holding cost ( ) :

(14)

HC__r := m*n*(int(h__r(t)*i__r(t), t = 0 .. t__4))/T

Cost due to deterioration () :

DC__r := P__d*`θ__r`*m*n*(int(i__r(t), t = 0 .. t__4))/T

Case 1R: N<=t4

IE__r1 := P__r*i__e*m*n*(int(R(N, t)*t, t = 0 .. N))/T

IC__r1 := P__d*i__c*m*n*(int(i__r(t), t = N .. t__4))/T

Case 2R: N > :

IE__r2 := P__r*i__e*m*n*(int(R(N, t)*t, t = 0 .. t__4)+N-t__4+int(R(N, t), t = 0 .. t__4))/T

params := [lambda = 3, a = 300, b = .15, c = .25, A__m = 300, A__d = 150, A__r = 50, C__m = 4, P__m = 8, P__d = 10, P__r = 12, `θ__m` = .15, `θ__d` = .12, `θ__r` = 0.5e-1, h__m = .2, h__d = .3, h__r = .5, i__m = .1, i__d = .1, i__r = .1, i__om = .1, i__OD = .15, i__c = .3, i__e = .2, M = 2, alpha = 0.2e-1]

ODE12 := diff(P1_TCS1, N) = 0

$set_m := {1, 2, 3, 4, 5}$

$set_n := {1, 2, 3, 4, 5}$

$for m in set_m do for n in set_n do SOL1[m, n] := fsolve({ODE11, ODE12}, {N, t__2}); t3Vals[m, n] := eval(t__3, [params[], SOL1[m, n][]]); t4Vals[m, n] := eval(t__4, [params[], SOL1[m, n][]]); if `and`(rhs(SOL1[m, n][1]) <= t4Vals[m, n], t3Vals[m, n] < rhs(params[25])) then d1[m, n] := eval(TCS__1, ['n' = n, 'm' = m, SOL1[m, n][]]); printf("Made it here at least once\n") end if end do end do$

SOL1(); t3Vals(); t4Vals()

(15)

Download badCalc2.mw

Looks like...

Answered: tomleslie 13876

July 22 2020

0 0

you haven't loaded the LinearAlgebra() package.

You need to enter either

with(LinearAlgebra):
Determinant( yourMatrix );

or

LinearAlgebra:-Determinant( yourMatrix );

When you submit a question you shoulld get into the habit of suplying an executable worksheet by using the big green up-arrow in the Mapleprimes toolba

Several typso...

Answered: tomleslie 13876

July 21 2020

0 3

(aside from the already-noted issue about labels) there are severaal typos - mainly because *sometimes* you have managed to enter an "inert" superscript, rather than mathematical exponentiation.

Fixing all of these, and generally tidying up, one ens up with a pair of equations in three variables, N, T, and t₂ The form of the equations is such that I think it very unlikely that solve() will return an analytic solution for N and t₂ in terms of T.

However as the attached shows, it is possible to use fsolve() to obtain numeric solutions for supplied numeric values of T.

>

i__m1 := proc (t) options operator, arrow; ((-c*t^2*`θ__m`^2+b*t*`θ__m`^2+2*c*t*`θ__m`-b*`θ__m`+`θ__m`^2-2*c)*exp(`θ__m`*t)*a*N^alpha*(lambda-1)/`θ__m`^3-(-b*`θ__m`+`θ__m`^2-2*c)*a*N^alpha*(lambda-1)/`θ__m`^3)*exp(-`θ__m`*t) end proc

>

i__m2 := proc (t) options operator, arrow; (-(-c*t^2*`θ__m`^2+b*t*`θ__m`^2+2*c*t*`θ__m`-b*`θ__m`+`θ__m`^2-2*c)*exp(`θ__m`*t)*a*N^alpha/`θ__m`^3+(-c*t__2^2*`θ__m`^2+b*t__2*`θ__m`^2+2*c*t__2*`θ__m`-b*`θ__m`+`θ__m`^2-2*c)*exp(`θ__m`*t__2)*a*N^alpha/`θ__m`^3)*exp(-`θ__m`*t) end proc

>

TC__m := A__m/(t__1+t__2)+(int(h__m*(i__m*t+1)*i__m1(t), t = 0 .. t__1))*(int(h__m*(i__m*t+1)*i__m2(t), t = 0 .. t__2))+P__m*I__om*m*(-(1/3)*a*c*N^alpha*M^3+(1/2)*a*b*N^alpha*M^2+a*N^alpha*M)/(t__1+t__2)+C__m*`θ__m`*(int(i__m1(t), t = 0 .. t__1)+int(i__m2(t), t = 0 .. t__2))/(t__1+t__2)

>

i__d := proc (t) options operator, arrow; (-(-c*t^2*`θ__d`^2+b*t*`θ__d`^2+2*c*t*`θ__d`-b*`θ__d`+`θ__d`^2-2*c)*a*N^alpha*exp(`θ__d`*t)/`θ__d`^3+(-c*t__3^2*`θ__d`^2+b*t__3*`θ__d`^2+2*c*t__3*`θ__d`-b*`θ__d`+`θ__d`^2-2*c)*a*N^alpha*exp(`θ__d`*t__3)/`θ__d`^3)*exp(-`θ__d`*t) end proc

>

TC__d1 := A__d*m/(t__1+t__2)+m*(int(h__d*(i__d*t+1)*i__d(t), t = 0 .. t__3))/(t__1+t__2)+P__d*I__OD*m*n*(-(1/3)*a*c*N^alpha*N^3+(1/2)*a*b*N^alpha*N^2+a*N^alpha*N)/(t__1+t__2)+P__m*`θ__m`*m*(int(i__d(t), t = 0 .. t__3))/(t__1+t__2)+P__m*I__c*m*(int(i__d(t), t = M .. t__3))/(t__1+t__2)-P__d*I__e*m*(-(1/3)*a*c*N^alpha*M^3+(1/2)*a*b*N^alpha*M^2+a*N^alpha*M)/(t__1+t__2)

>

i__r := proc (t) options operator, arrow; (-(-c*t^2*`θ__r`^2+b*t*`θ__r`^2+2*c*t*`θ__r`-b*`θ__r`+`θ__r`^2-2*c)*a*N^alpha*exp(`θ__r`*t)/`θ__r`^3+(-c*t__4^2*`θ__r`^2+b*t__4*`θ__r`^2+2*c*t__4*`θ__r`-b*`θ__r`+`θ__r`^2-2*c)*a*N^alpha*exp(`θ__r`*t__4)/`θ__r`^3)*exp(-`θ__r`*t) end proc

>

TC__r1 := A__r*m*n/(t__1+t__2)+m*n*(int(h__r*(i__r*t+1)*i__r(t), t = 0 .. t__4))/(t__1+t__2)+P__d*`θ__r`*m*n*(int(i__r(t), t = 0 .. t__4))/(t__1+t__2)+P__d*I__c*m*n*(int(i__r(t), t = N .. t__4))/(t__1+t__2)-P__r*I__e*m*n*(-(1/4)*a*c*N^alpha*N^4+(1/3)*a*b*N^alpha*N^3+(1/2)*a*N^alpha*N^2)/(t__1+t__2)

>

TC__r2 := A__r*m*n/(t__1+t__2)+m*n*(int(h__r*(i__r*t+1)*i__r(t), t = 0 .. t__4))/(t__1+t__2)+P__d*`θ__r`*m*n*(int(i__r(t), t = 0 .. t__4))/(t__1+t__2)-P__r*I__e*m*n*(-(1/4)*a*c*N^alpha*t__4^4+(1/3)*a*b*N^alpha*t__4^3+(1/2)*a*N^alpha*t__4^2+N-t__4-(1/3)*a*c*N^alpha*t__4^3+(1/2)*a*b*N^alpha*t__4^2+a*N^alpha*t__4)/(t__1+t__2)

>

t__1 := (1/6)*t__2*(2*b*t__2^2*`θ__m`-2*c*t__2^2+3*b*t__2+3*t__2*`θ__m`+6)/(lambda-1)

${b, c, lambda, t__2, theta__m}$

(1)

>

TCS1_param_1 := eval(TCS__1, [A__d = 150, A__m = 300, A__r = 50, C__m = 4, I__OD = .15, I__c = .3, I__e = .2, I__om = .1, M = 2, P__d = 10, P__m = 8, P__r = 12, a = 300, alpha = 0.2e-1, b = .15, c = .25, h__d = .3, h__m = .2, h__r = .5, i__d = .1, i__m = .1, i__r = .1, lambda = 3, m = 1, n = 4, `θ__d` = .12, `θ__m` = .15, `θ__r` = 0.5e-1])

${N, T, t__2}$

(2)

>

${N, T, t__2}$

(3)

>

ODE2 := diff(TCS1_param_1, t__2) = 0

${N, T, t__2}$

(4)

>

>	$seq([T = j, fsolve(eval({ODE1, ODE2}, T = j), {N, t__2})], j = 1 .. 10)$

$[T = 1, {N = 0.3169610207e-1, t__2 = -1.734267919}], [T = 2, {N = .2921764574, t__2 = 4.749386230}], [T = 3, {N = 1.118168228, t__2 = -2.511574550}], [T = 4, {N = 1.653457879, t__2 = -2.804382830}], [T = 5, {N = 1.982530174, t__2 = 4.748449437}], [T = 6, {N = 2.093187798, t__2 = 2.635830725}], [T = 7, {N = 2.446225114, t__2 = 4.743052309}], [T = 8, {N = 2.614318363, t__2 = 2.944845840}], [T = 9, {N = 2.867688641, t__2 = 3.155099720}], [T = 10, {N = 3.097722967, t__2 = -3.075355988}]$

(5)

>

Download getSol.mw

The best approach...

Answered: tomleslie 13876

July 20 2020

2 1

may depend on whether individual entries in the matrix depend on a single variable, or multiple variables. The attached, with a "toy" example handles both cases

Download matLim.mw

Is this...

Answered: tomleslie 13876

July 19 2020

0 1

what you were expecting?

As an aside, soething seems to have gone a bit wrng with rendering of long expressions on this site - the attached is perfectly legible within Maple - honest!

>

restart;

i__m1(t) = ((-c*t^2*`θ__m`^2+b*t*`θ__m`^2+2*c*t*`θ__m`-b*`θ__m`+`θ__m`^2-2*c)*exp(`θ__m`*t)*a*N^alpha*(lambda-1)/`θ__m`^3-(-b*`θ__m`+`θ__m`^2-2*c)*a*N^alpha*(lambda-1)/`θ__m`^3)*exp(-`θ__m`*t)

i__m1(t) = ((-c*t^2*theta__m^2+b*t*theta__m^2+2*c*t*theta__m-b*theta__m+theta__m^2-2*c)*exp(theta__m*t)*a*N^alpha*(lambda-1)/theta__m^3-(-b*theta__m+theta__m^2-2*c)*a*N^alpha*(lambda-1)/theta__m^3)*exp(-theta__m*t)

(1)

i__m2(t) = (-(-c*t^2*`θ__m`^2+b*t*`θ__m`^2+2*c*t*`θ__m`-b*`θ__m`+`θ__m`^2-2*c)*exp(`θ__m`*t)*a*N^alpha/`θ__m`^3+(-c*t__2^2*`θ__m`^2+b*t__2*`θ__m`^2+2*c*t__2*`θ__m`-b*`θ__m`+`θ__m`^2-2*c)*exp(`θ__m`*t__2)*a*N^alpha/`θ__m`^3)*exp(-`θ__m`*t)

i__m2(t) = (-(-c*t^2*theta__m^2+b*t*theta__m^2+2*c*t*theta__m-b*theta__m+theta__m^2-2*c)*exp(theta__m*t)*a*N^alpha/theta__m^3+(-c*t__2^2*theta__m^2+b*t__2*theta__m^2+2*c*t__2*theta__m-b*theta__m+theta__m^2-2*c)*exp(theta__m*t__2)*a*N^alpha/theta__m^3)*exp(-theta__m*t)

(2)

TC__m := A__m/(t__1+t__2)+(int(h__m*(i__m*t+1)*i__m1(t), t = 0 .. t__1))*(int(h__m*(i__m*t+1)*i__m2(t), t = 0 .. t__2))+P__m*I__om*m*(-(1/3)*a*c*N^alpha*M^3+(1/2)*a*b*N^alpha*M^2+a*N^alpha*M)/(t__1+t__2)+C__m*`θ__m`*(int(i__m1(t), t = 0 .. t__1)+int(i__m2(t), t = 0 .. t__2))/(t__1+t__2)

A__m/(t__1+t__2)+(int(h__m*(i__m*t+1)*i__m1(t), t = 0 .. t__1))*(int(h__m*(i__m*t+1)*i__m2(t), t = 0 .. t__2))+P__m*I__om*m*(-(1/3)*a*c*N^alpha*M^3+(1/2)*a*b*N^alpha*M^2+a*N^alpha*M)/(t__1+t__2)+C__m*theta__m*(int(i__m1(t), t = 0 .. t__1)+int(i__m2(t), t = 0 .. t__2))/(t__1+t__2)

(3)

i__d(t) = (-(-c*t^2*`θ__d`^2+b*t*`θ__d`^2+2*c*t*`θ__d`-b*`θ__d`+`θ__d`^2-2*c)*a*N^alpha*exp(`θ__d`*t)/`θ__d`^3+(-c*t__3^2*`θ__d`^2+b*t__3*`θ__d`^2+2*c*t__3*`θ__d`-b*`θ__d`+`θ__d`^2-2*c)*a*N^alpha*exp(`θ__d`*t__3)/`θ__d`^3)*exp(-`θ__d`*t)

i__d(t) = (-(-c*t^2*theta__d^2+b*t*theta__d^2+2*c*t*theta__d-b*theta__d+theta__d^2-2*c)*a*N^alpha*exp(theta__d*t)/theta__d^3+(-c*t__3^2*theta__d^2+b*t__3*theta__d^2+2*c*t__3*theta__d-b*theta__d+theta__d^2-2*c)*a*N^alpha*exp(theta__d*t__3)/theta__d^3)*exp(-theta__d*t)

(4)

TC__d1 := A__d*m/(t__1+t__2)+m*(int(h__d*(i__d*t+1)*i__d(t), t = 0 .. t__3))/(t__1+t__2)+P__d*I__OD*m*n*(-(1/3)*a*c*N^alpha*N^3+(1/2)*a*b*N^alpha*N^2+a*N^alpha*N)/(t__1+t__2)+P__m*`θ__m`*m*(int(i__d(t), t = 0 .. t__3))/(t__1+t__2)+P__m*I__c*m*(int(i__d(t), t = M .. t__3))/(t__1+t__2)-P__d*I__e*m*(-(1/3)*a*c*N^alpha*M^3+(1/2)*a*b*N^alpha*M^2+a*N^alpha*M)/(t__1+t__2)

A__d*m/(t__1+t__2)+m*(int(h__d*(i__d*t+1)*i__d(t), t = 0 .. t__3))/(t__1+t__2)+P__d*I__OD*m*n*(-(1/3)*a*c*N^alpha*N^3+(1/2)*a*b*N^alpha*N^2+a*N^alpha*N)/(t__1+t__2)+P__m*theta__m*m*(int(i__d(t), t = 0 .. t__3))/(t__1+t__2)+P__m*I__c*m*(int(i__d(t), t = M .. t__3))/(t__1+t__2)-P__d*I__e*m*(-(1/3)*a*c*N^alpha*M^3+(1/2)*a*b*N^alpha*M^2+a*N^alpha*M)/(t__1+t__2)

(5)

A__d*m/(t__1+t__2)+m*(int(h__d*(i__d*t+1)*i__d(t), t = 0 .. t__3))/(t__1+t__2)+P__d*I__OD*m*n*(-(1/3)*a*c*N^alpha*N^3+(1/2)*a*b*N^alpha*N^2+a*N^alpha*N)/(t__1+t__2)+P__m*theta__m*m*(int(i__d(t), t = 0 .. t__3))/(t__1+t__2)-P__d*I__e*m*(-(1/4)*a*c*N^alpha*t__3^4+(1/3)*a*b*N^alpha*t__3^3+(1/2)*a*N^alpha*t__3^2+M-t__3-(1/3)*a*c*N^alpha*t__3^3+(1/2)*a*b*N^alpha*t__3^2+a*N^alpha*t__3)/(t__1+t__2)

(6)

i__r(t) = (-(-c*t^2*`θ__r`^2+b*t*`θ__r`^2+2*c*t*`θ__r`-b*`θ__r`+`θ__r`^2-2*c)*a*N^alpha*exp(`θ__r`*t)/`θ__r`^3+(-c*t__4^2*`θ__r`^2+b*t__4*`θ__r`^2+2*c*t__4*`θ__r`-b*`θ__r`+`θ__r`^2-2*c)*a*N^alpha*exp(`θ__r`*t__4)/`θ__r`^3)*exp(-`θ__r`*t)

i__r(t) = (-(-c*t^2*theta__r^2+b*t*theta__r^2+2*c*t*theta__r-b*theta__r+theta__r^2-2*c)*a*N^alpha*exp(theta__r*t)/theta__r^3+(-c*t__4^2*theta__r^2+b*t__4*theta__r^2+2*c*t__4*theta__r-b*theta__r+theta__r^2-2*c)*a*N^alpha*exp(theta__r*t__4)/theta__r^3)*exp(-theta__r*t)

(7)

TC__r1 := A__r*m*n/(t__1+t__2)+m*n*(int(h__r*(i__r*t+1)*i__r(t), t = 0 .. t__4))/(t__1+t__2)+P__d*`θ__r`*m*n*(int(i__r(t), t = 0 .. t__4))/(t__1+t__2)+P__d*I__c*m*n*(int(i__r(t), t = N .. t__4))/(t__1+t__2)-P__r*I__e*m*n*(-(1/4)*a*c*N^alpha*N^4+(1/3)*a*b*N^alpha*N^3+(1/2)*a*N^alpha*N^2)/(t__1+t__2)

A__r*m*n/(t__1+t__2)+m*n*(int(h__r*(i__r*t+1)*i__r(t), t = 0 .. t__4))/(t__1+t__2)+P__d*theta__r*m*n*(int(i__r(t), t = 0 .. t__4))/(t__1+t__2)+P__d*I__c*m*n*(int(i__r(t), t = N .. t__4))/(t__1+t__2)-P__r*I__e*m*n*(-(1/4)*a*c*N^alpha*N^4+(1/3)*a*b*N^alpha*N^3+(1/2)*a*N^alpha*N^2)/(t__1+t__2)

(8)

TC__r2 := A__r*m*n/(t__1+t__2)+m*n*(int(h__r*(i__r*t+1)*i__r(t), t = 0 .. t__4))/(t__1+t__2)+P__d*`θ__r`*m*n*(int(i__r(t), t = 0 .. t__4))/(t__1+t__2)-P__r*I__e*m*n*(-(1/4)*a*c*N^alpha*t__4^4+(1/3)*a*b*N^alpha*t__4^3+(1/2)*a*N^alpha*t__4^2+N-t__4-(1/3)*a*c*N^alpha*t__4^3+(1/2)*a*b*N^alpha*t__4^2+a*N^alpha*t__4)/(t__1+t__2)

A__r*m*n/(t__1+t__2)+m*n*(int(h__r*(i__r*t+1)*i__r(t), t = 0 .. t__4))/(t__1+t__2)+P__d*theta__r*m*n*(int(i__r(t), t = 0 .. t__4))/(t__1+t__2)-P__r*I__e*m*n*(-(1/4)*a*c*N^alpha*t__4^4+(1/3)*a*b*N^alpha*t__4^3+(1/2)*a*N^alpha*t__4^2+N-t__4-(1/3)*a*c*N^alpha*t__4^3+(1/2)*a*b*N^alpha*t__4^2+a*N^alpha*t__4)/(t__1+t__2)

(9)

(10)

(11)

(12)

(13)

TCS := piecewise(M <= t__3 and N <= t__4, TCS__1, `and`(M <= t__3, N > t__4), TCS__2, `and`(M >= t__3, N <= t__4), TCS__3, `and`(M >= `t__3 `, N > t__4), TCS__4)

(14)

(15)

(16)

Download weirdDiff.mw

Another attempt...

Answered: tomleslie 13876

July 19 2020

1 1

It is not clear to me what output you desire - graphs or tables

Whichever of thed above you choose, how do you want the data organised?

for example

a graph of f(x), g(x), t(x) versus the independent variable x for a supplied value of the parameter ax?

a table of f(x), g(x), t(x) versus the parameter ax, for a supplied value of the independent variable 'x'?

The attached shows multiple possibilities - one of which ought(?) to meet your requirements.

>

m := .2; pa := 3.14*(1/3); aa := .1; bb := .3; ta := .2; kt := .4; h2 := 1+m*ax+bb*sin((2*3.14)*(ax-ta)); h1 := -1-m*ax-aa*sin((2*3.14)*(ax-ta)+pa); a2 := kt+aa*sin((2*3.14)*(ax-ta))+bb*sin((2*3.14)*(ax-ta)+pa)

(1)

>

1.5*(diff(diff(diff(diff(f(x), x), x), x), x))-.40*alpha*(6*(diff(diff(f(x), x), x))*(diff(diff(diff(f(x), x), x), x))^2+3*(diff(diff(f(x), x), x))^2*(diff(diff(diff(diff(f(x), x), x), x), x)))-.20*(diff(diff(f(x), x), x))+.20*(diff(t(x), x))+.20*(diff(g(x), x)) = 0

(2)

>

1.75*(diff(diff(t(x), x), x))+(diff(t(x), x))*(diff(g(x), x))+2*(diff(t(x), x))^2+.7*(diff(diff(f(x), x), x))^2-.7*(diff(diff(f(x), x), x))^4+.1*(diff(f(x), x))^2+.1*t(x) = 0

(3)

>

diff(diff(g(x), x), x)+2*(diff(diff(t(x), x), x))-.5*g(x) = 0

(4)

>

#
# Generate solution for values of 'ax' in the
# range 0..1 in step of 0.1, with the exceptions
# of ax=0.2 and ax=0.3 - which (for reasons I can't
# figure out) cause problems
#
  L:= {seq(j, j=0..1,0.1)} minus {0.2, 0.3}:
  NN := nops(L):
  for k from 1 to NN do
      R[L[k]]:= dsolve
                ( eval({bc, eq1,eq2,eq3}, ax= L[k]),
                  fcns,
                  type = numeric,
                  maxmesh=8192,
                  continuation=alpha,
                  output=listprocedure
                );
      lbound:=op(1, lhs(eval(bc, ax=L[k])[3]));
      ubound:=op(1, lhs(eval(bc, ax=L[k])[1]));
      plt[L[k]]:= odeplot
                  ( R[L[k]],
                    [ [ x, f(x) ],
                      [ x, g(x) ],
                      [ x, t(x) ]
                    ],
                    x = lbound..ubound
                  );
      vals[L[k]]:= Matrix
                   ( [ ['x', 'f(x)', 'g(x)', 't(x)'],
                       seq( eval( [x, f(x), g(x), t(x)],
                                  R[L[k]]
                                )(j),
                            j=lbound..ubound,(ubound-lbound)/20
                          )
                     ]
                   );
#
# Following command implies OP wants to know
#
# diff(f(x),x,x) evaluated at x=0
#
# No idea why this is significant - but if that
# is what is desired .......
#
      X1[L[k]]:=rhs(R[L[k]](0)[4]);
  end do:

>	############################################# # Some example output possibilities #############################################

>	# # Display all of the obtained solution curves # display( [entries(plt, 'nolist')]);

>

#
# Define a function which will plot f(x), g(x)
# t(x) as functions of the independent variable
# x, for the supplied value of the parameter ax
#
  gD0:= z-> `if`( member(z, L),
                  display
                  ( plt[z],
                    title=typeset( "f(x), g(x), t(x) versus x for ax = ", z),
                    titlefont=[times, bold, 20]
                  ),
                  cat( "Data unavailable for ax=",
                       convert(z,string)
                     )
                ):
#
# Example of use
#
  gD0(0.5);
  gD0(0.3);

(5)

>

#
# Define function which will return data for
# any supplied value of the parameter ax
#
  gD1:= z-> `if`( member(z, L),
                  vals[z],
                  cat( "Data unavailable for ax=",
                        convert(z,string)
                     )
                ):
#
# Examples of use
#
  interface(rtablesize=25):
  gD1(0.5);
  gD1(0.2);
  interface(rtablesize=10):

(6)

>

#
# Define a function to produce
#
#       f(x) vs ax,
#       g(x) vs ax,
#       t(x) vs ax,
#
# at a supplied value of the independent
# variable 'x'
#
  gD2:=z-> Matrix
            ( [ ['ax', 'f'(z), 'g'(z), 't'(z)],
                seq( eval( [L[k], f(x), g(x), t(x)],
                           R[L[k]]
                         )(z),
                     k=1..NN
                   )
              ]
            ):
#
# Examples of use
#
  gD2(-0.5);
  gD2(0);
  gD2(0.5);

Matrix(%id = 18446744074575034598)

(7)

>

#
# Define a function to plot
#
#       f(x) vs ax,
#       g(x) vs ax,
#       t(x) vs ax,
#
# at a supplied value of the independent
# variable 'x'
#
  moreOpts:= style=point,
             symbol=solidcircle,
             symbolsize=20,
             labels=[ typeset("ax"),
                      typeset("f(x),g(x),t(x)")
                    ],
             labelfont=[times, italic, 16]:
  cols:=[red, green, blue]:

  gD3:=z-> display
           ( [ seq
               ( plot
                 ( [seq( L[k], k=1..NN)],
                   [seq( eval( fcns[j], R[L[k]])(z), k=1..NN)],
                   color=cols[j],
                   moreOpts
                 ),
                 j=1..3
               )
            ],
            title=typeset( "f(x), g(x), t(x) versus ax at x = ", z),
            titlefont=[times, bold, 20]
           ):
#
# Examples of use
#
  gD3(-0.5);
  gD3(0.0);
  gD3(0.5);

>

Download odeCont3.mw

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

13876 Reputation

20 Badges

MaplePrimes Activity

These are answers submitted by tomleslie

Diififult to work out...

No analytic solution???...

That is ne "ugly" equation...

The answer is not 76 +2/25...

try the geom3d package...

Use the option 'listprocedure'...

Several typos...

Use the ImportMatrix() command...

Don't think...

Hmmmm...

Looks like...

Several typso...

The best approach...

Is this...

Another attempt...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again