## 5473 Reputation

9 years, 338 days

## I should have noted...

Although I did the original worksheet in Maple 2016 (becuase that's what the OP has), the behaviour in Maple 2019 is identical :-(

## The concepts...

of namespace', function overloading and variable scope are basic features of any programming language. If you really do not understand these concepts then I respectfully suggest that you keep reading the attached toy example until enlightenment occurs

 > restart; # # Define som simple "top-level" functions # and check how they evaluate #   lambda:= x->x^2:   phi:= x->x^3:   pi:= x->x^4:   lambda(2), phi(2), pi(2); # # Now load the NumberTheory package which contains # its own definitions for the functions lambda, phi, # and pi. These will 'overload' the previous # definitions, so check how these now evaluate #   with(NumberTheory):   lambda(2), phi(2), pi(2); # # If you wish to retrieve/use the function definitions # performed at the top-level, rather than the overloaded # versions introduced by the NumberTheory package, then # prefix the function names with ':-', as in #   :-lambda(2), :-phi(2), :-pi(2);
 (1)
 >

## It's a version issue...

Attached is a minor revision of my original post which runs with no problems in Mpale 2018

 > restart:   kernelopts(version);   xmin:=1: xmax:=3: ymin:=2.5: ymax:=4.5:   f:= (nx, ny)-> plots:-pointplot                         ( [ seq                             ( seq                               ( [i, j],                                 i=xmin..xmax, (xmax-xmin)/nx                               ),                               j=ymin..ymax, (ymax-ymin)/ny                             )                           ],                           symbol=solidcircle,                           symbolsize=20,                           color=red,                           scaling=constrained                         ):
 (1)
 > # # A few test cases for different grid spacings #   f(4,4);   f(4,6);   f(10,3);
 >

## Can't reproduce...

provided commands are issued in the following sequence

fopen(yourFileName)
writeTo(yourFileName)
fclose(yourFileName)

The file will be opened (and "locked") when the fopen() command is executeded. It will remain open (and "locked") until the fclose() command is executed.

Only thing I can think of is that you are trying to access the file  before the executing fclose().

## It depends on your usage...

You could do it either way!

If you are going to run this on multiple data sets, then creating a procedure within the package to perform the plots will probably save you some typing

Within the package, I'd probably separate the "fitting" and the "plotting" into two separate procedures - because sometimes one might want the former and not the latter.

See the attached, which is a slight extension of the worksheet I posted previously

 > restart;   Mat:= module()                description "My Package";                option package;                export fitData, doPlots;                fitData:=proc( xd::list, yd::list, fittype::name,                               x::algebraic, c::algebraic                             )                             uses Statistics, CurveFitting:                             if   fittype = modfit                             then return Fit(c, xd, yd, x);                             elif fittype = lsqfit                             then return LeastSquares( xd, yd, x, curve=c)                             fi:                end proc;                doPlots:=proc( xd::list, yd::list, cur::algebraic )                             uses plots:                             display                             (  [ plot                                  ( cur,                                    x=min(xl)..max(xl),                                    color=blue,                                    title=typeset("Fit = ", cur)                                  ),                                  plot                                  ( xl, yl,                                    style=point,                                    symbol=solidcircle,                                    symbolsize=20,                                    color=red                                  )                                ],                                size=[600,600]                             );                         end proc;         end module:   with(Mat): # # test data #   xl:=[1,2,3,4]: yl:=[2,4,5,9]: # # Use Statistics:-Fit() with both linear # and quadratic models #   f1:= fitData(xl, yl, 'modfit', x, 'a*x+b');   doPlots(xl, yl, f1);   f2:= fitData(xl, yl, 'modfit', x, 'a*x^2+b*x+c');   doPlots(xl, yl, f2); # # Use CurveFitting:-LeastSquares() with # both linear and quadratic models, then # plot the fit #     f3:= fitData(xl, yl, 'lsqfit', x, 'a*x+b');   doPlots(xl, yl, f3);   f4:= fitData(xl, yl, 'lsqfit', x, 'a*x^2+b*x+c');   doPlots(xl, yl, f4);
 >

## Using Rouben's suggestion...

It is possible to produce a 'symbolic' solution for the ODE - see the attached - but maybe this is too complicated to be useful?

This particular problem seems to fall into the category of those which should be solved numerically!

 > restart; #B:=10; #H:=10; #L__1:=4; eq0:=R__i=(1/2)*B+L__1; eq01:=R__f=C+(1/2)*B+L__1; eq02:=theta__1=arctan((H/2)/(B/2 + L__1)); eq03:=theta__2=Pi - arctan((H/2)/(B/2 - L__1)); eq43:= v__beta(r,theta,beta)= v__m(beta)*(1-(r^2/((R__max(theta,beta))^2))); eq39:=R(beta)=R__i+2*beta*(R__f-R__i)/Pi; eq39a:=simplify(subs([eq0,eq01], eq39)); eq38:=R__max(theta,beta)=R__max(theta,0)*(R(beta)/R__i); # # Replace eq41 with the simple requirement that # R__max(theta,0)=phi(theta - Rouben's suggestion # eq41_orig:=R__max(theta,0)=piecewise(theta<=rhs(eq02),(B/2+L__1)/cos(theta),theta<=rhs(eq03),(H/2)/sin(theta),rhs(eq03)= 2 , H > 0 , B > 0 , L__1 > 0; eq44b:= v__m(beta)= rhs(eq44a) * B * H; eq47:=subs([eq44b,eq38a],eq43);
 (1)
 > eq37:=diff(v__r(r, theta, beta),r)*r*(R__f-r*cos(theta))+v__r(r,theta,beta)*(R__f-2*r*cos(theta))+r*diff(v__beta(r,theta,beta),beta)=0; eq37a:=simplify(expand(subs(eq47, eq37)));
 (2)
 > bc:=v__r(rhs(eq38a),theta,beta)=0;
 (3)
 > sol3:=simplify(dsolve([eq37a,bc], v__r(r, theta,beta)));
 (4)
 > # # Having obtained a solution in terms of phi(theta) # substitute the original piecewise definition from # eq41_orig # # The result is as ugly as sin! #   sol4:=eval( sol3, phi(theta)=rhs(eq41_orig));
 (5)
 > # # Evaluate this solution for specific B, H, L__1 # # This result is still ugly #   eval(sol4, [B=10, H=10, L__1=4]);
 (6)
 >

## Don't waste time...

OP has posted this question before, received an answer from me, and then deleted the whole thread. Why? I do not know!

As well as the multiple recursive assignments, the OP has no understanding of of the distinction between  inert and "ert" subscripts, so the whole worksheet is littered with confusions between varName__n and varName[n]

I fixed all of this once, but since the OP has deleted my original fix, I don't propose to bother doing it again

## The obvious solution...

is to supply the required boundary condition to the dsolve() command, which I have done in the final execution group of the attached

 > restart; d := 10; c := 10; L__1 := 4; R__i := (1/2)*d+L__1; R__f := c+(1/2)*d+L__1; eq43:= v__beta(r,theta,beta)= v__m(beta)*(1-(r^2/(R__max(theta,beta)^2))); eq37:=diff(v__r(r, theta, beta),r)*r*(R__f-r*cos(theta))+v__r(r,theta,beta)*(R__f-2*r*cos(theta))+r*diff(v__beta(r,theta,beta),beta)=0;
 (1)
 > eq39:=R(beta)=R__i+2*beta*(R__f-R__i)/Pi;
 (2)
 > eq38:=R__max(theta,beta)=R__max(theta,0)*(R(beta)/R__i); eq41:=R__max(theta,0)=((2*L__1*cos(theta)) + sqrt(d^2 - 2*(L__1)^2 + 2*(L__1)^2*cos(2*theta)))/2;
 (3)
 > eq38a:= subs(eq41, eq38);
 (4)
 > eq38b:=simplify(subs(eq39,eq38a));
 (5)
 > eq42:=R__max(theta,beta)=( (((d+2*L__1)*Pi) - 2*beta*(d+2*L__1-2*R__f))*(((2*L__1*cos(theta)) + sqrt(d^2 - 2*(L__1)^2 + 2*(L__1)^2*cos(2*theta)))))/(2*Pi*(d+2*L__1));
 (6)
 >
 (7)
 > eq44a:=simplify(expand(subs(eq38b,eq44)));
 (8)
 > eq47:=subs([eq44a,eq38b],eq43);
 (9)
 > eq37a:=simplify(expand(subs(eq47, eq37)));
 (10)
 > sol1:=simplify(dsolve(eq37a, v__r(r, theta,beta)), build);
 (11)
 > sol2:=eval(sol1, _F1(theta, beta)=0);
 (12)
 >
 (13)
 >
 (14)
 > eq48a := simplify(expand(subs([eq0, eq00], eq48)));
 (15)
 > # # Check Maple's solution and OP's proposed solution #   simplify(expand(subs( sol2, eq37a)));   simplify(expand(subs( eq48a, eq37a)));
 (16)
 > # # Construct a boundary condition involving R__max #   bc:=v__r(simplify(rhs(eq42)), theta,beta)=0; # # Solve the ODE with this bc #   sol3:=simplify(dsolve([eq37a,bc], v__r(r, theta,beta))); # # verify that this is indeed a solution #   simplify(expand(subs(sol3, eq37a))); # # Verify that the boundary condition is satisfied #   simplify(eval(sol3, r=op([1,1],bc) ));
 (17)
 >

## A quick check...

You state

eq48a is from the paper and sol2 should be equal to eq48a.

Now one (or both??) of these should be the solution of the ODE you have obtained (ie eq37a) - so try substituting both eq48a and sol2 in the ODE eq37a and see what happens - as in the attached.

Substituting 'sol2' in eq37a returns 0=0, and eq48a does not, so the former is a solution and the latter is not. Now it may be possible that some restrictions/assumptions on the variables r, theta, beta (which you are not applying) will convert eq48a to a "solution" - but only you know what these conditions/restrictions might be

 > restart; d := 10; c := 10; L__1 := 4; R__i := (1/2)*d+L__1; R__f := c+(1/2)*d+L__1; eq43:= v__beta(r,theta,beta)= v__m(beta)*(1-(r^2/(R__max(theta,beta)^2))); eq37:=diff(v__r(r, theta, beta),r)*r*(R__f-r*cos(theta))+v__r(r,theta,beta)*(R__f-2*r*cos(theta))+r*diff(v__beta(r,theta,beta),beta)=0;
 (1)
 > eq39:=R(beta)=R__i+2*beta*(R__f-R__i)/Pi;
 (2)
 > eq38:=R__max(theta,beta)=R__max(theta,0)*(R(beta)/R__i); eq41:=R__max(theta,0)=((2*L__1*cos(theta)) + sqrt(d^2 - 2*(L__1)^2 + 2*(L__1)^2*cos(2*theta)))/2;
 (3)
 > eq38a:= subs(eq41, eq38);
 (4)
 > eq38b:=simplify(subs(eq39,eq38a));
 (5)
 > eq42:=R__max(theta,beta)=( (((d+2*L__1)*Pi) - 2*beta*(d+2*L__1-2*R__f))*(((2*L__1*cos(theta)) + sqrt(d^2 - 2*(L__1)^2 + 2*(L__1)^2*cos(2*theta)))))/(2*Pi*(d+2*L__1));
 (6)
 >
 (7)
 > eq44a:=simplify(expand(subs(eq38b,eq44)));
 (8)
 > eq47:=subs([eq44a,eq38b],eq43);
 (9)
 > eq37a:=simplify(expand(subs(eq47, eq37)));
 (10)
 > sol1:=simplify(dsolve(eq37a, v__r(r, theta,beta)), build);
 (11)
 > sol2:=eval(sol1, _F1(theta, beta)=0);
 (12)
 >
 (13)
 >
 (14)
 >
 (15)
 > simplify(expand(subs( eq48a, eq37a)));
 (16)
 > simplify(expand(subs( sol2, eq37a)));
 (17)
 >

## So...

I have shown you how to solve the equation system

You now want to use different definitions for some of the equations. Why don't you try following the logic of my original post post and typing in the relevant equations yourself?

## Post the worksheet...

use the big green up-arrow in the Mapleprimes toolbar to upload the problem worksheet - otherwise all anyone can do is guess!

## In Maple 18...

The attached seems to run correctly in Maple 18.

I had to make even more edits which you will have to check. The biggest problem was the use of 2-D math input, whose idiosyncrasies divide opinion even now. In Maple 18 (ie five year-old software), the implementation of 2-D math input appears to be even more temperamental! So much so, that I gave up fighting with it and converted the whole worksheet to 1-D input; simple-to-use and pretty much bulletproof.

One other "feature" which slightly bothers me is that if I believe the "labels" which you have supplied for the two calculations, then it would seem that these calculations are actually different. The final two "labels" in the first calculation are "Number of Available Vacancies" and  "Employed Individuals E(t)". The last two labels in the second calculation are "Unemployment U(t)", and  "Poverty P[R](t)" - so definitely different.

Since all of the dependent variables are coupled, this discrepancy (if real, rather than just a typo invoving label names) will "propagate" to all other dependent variables, and (maybe?) good agreement betwee the two calcuations is not to be expected???

Anyhow, you may want to consider the attached. (and since someone has fixed the in-line display, at least you can see that everything seems to be working)

 > restart:
 > with(plots):   interface(version);   local gamma:
 (1)
 > # # Define parameters #   M__h:= .50: beta[o]:= 0.34e-1: beta[1]:= 0.25e-1:   mu[r]:= 0.4e-3: sigma:= .7902: alpha:= .11:   psi:= 0.136e-3: xi:= 0.5e-1: gamma:= .7:   M__c:= .636: mu[b]:= 0.5e-2: omega:= .134:
 > # # Some initial values #   B[0]:= .50: C[0]:= .30: DD[0]:= .21: E[0]:= .14:   F[0]:= .70: G[0]:= .45: H[0]:= .14:
 > # # Define ODEs, initial conditions, and solve #   ODEs:= { diff(J[1](T), T) = M__h-beta[1]*psi*(J[1](T)+J[2](T))*J[7](T)-sigma*psi*beta[1]*J[4](T)*J[7](T)-mu[r]*J[1](T),            diff(J[2](T), T) = beta[1]*psi*(J[1](T)+J[2](T))*J[7](T)-(alpha+xi+mu[r])*J[2](T),            diff(J[3](T), T) = alpha*J[2](T)-(omega+mu[r])*J[3](T),            diff(J[4](T), T) = omega*J[3](T)-(gamma+mu[r])*J[4](T),            diff(J[5](T), T) = gamma*J[4](T)+sigma*psi*beta[1]*J[5](T)*J[7](T)-mu[r]*J[5](T),            diff(J[6](T), T) = M__c-psi*beta[o]*J[6](T)*J[3](T)-mu[b]*J[6](T),            diff(J[7](T), T) = psi*beta[o]*J[6](T)*J[3](T)-mu[b]*J[7](T)          }:   ic1:= { J[1](0) = B[0], J[2](0) = C[0], J[3](0) = DD[0],           J[4](0) = E[0], J[5](0) = F[0], J[6](0) = G[0],           J[7](0) = H[0]          }:   sol1:= dsolve( `union`(ODEs, ic1),                   numeric                 ):
 > # # Generate a simple display of numerical results # from the above calculation. # # Need to juggle the default display limits # (then reset them to defaults) #   interface(displayprecision=4):   interface(rtablesize=20):   M:= Matrix( [ lhs~(sol1(0)),                 seq( rhs~(sol1(j)), j=0..10)               ]             );   interface(displayprecision=-1):   interface(rtablesize=10):
 (2)
 > # # Generate plots of the above solution #   plotopts1:= T = 0 .. 10, color = red, labeldirections = [horizontal, vertical]:   ylabels1:= [ "Gross National Income G[I](t)",                "Gross National Savings G[S](t)",                "Gross Domestic Products G[D](t)",                "Public Expenditure on Education P[E](t)",                "Energy Consumption E[C](t)",                "Number of Available Vacancies",                "Employed Individuals E(t)"              ]:   v:= [ seq          ( odeplot            ( sol1,              [T, J[k](T)],              plotopts1,              labels = ["Time in years", ylabels1[k]]            ),            k = 1 .. 7          )        ]:   display(v[1]);   display(v[2]);   display(v[3]);   display(v[4]);   display(v[5]);   display(v[6]);   display(v[7]);
 > for k from 0 to 3 do       B[k+1]:= (-B[k]*mu[r]+M__h-beta[1]*psi*add(B[k]*H[k-j], j=0..k)-beta[1]*psi*add(C[k]*H[k-j], j=0..k)-beta[1]*sigma*psi*add(F[k]*G[k-j], j=0..k))/(k+1);       C[k+1]:= (-(alpha+mu[r])*C[k]+beta[1]*psi*add(B[k]*H[k-j], j=0..k)+beta[1]*psi*add(C[k]*H[k-j], j=0..k))/(k+1);       DD[k+1]:= (alpha*C[k]-(omega+sigma+mu[r])*DD[k])/(k+1);       E[k+1]:= (omega*DD[k]-(gamma+mu[r])*E[k])/(k+1);       F[k+1]:= (gamma*E[k]-mu[r]*F[k]-beta[1]*sigma*psi*add(F[k]*G[k-j], j = 0 .. k))/(k+1);       G[k+1]:= (M__c-mu[b]*G[k]-beta[o]*psi*add(G[k]*C[k-j], j = 0 .. k))/(k+1);       H[k+1]:= (mu[b]*H[k]-beta[o]*psi*add(G[k]*C[k-j], j = 0 .. k))/(k+1):   end do:   g:= X-> unapply(sum(X[i]*z^i, i=0..1), z):   f:= X-> unapply(sum(X[i]*z^i, i=0..4), z):   fns:=['B','C', 'DD', 'E', 'F', 'G', 'H']:
 > # # Generate numerical results from this second calculation #   interface(displayprecision=4):   interface(rtablesize=20):   Matrix( [ ['T', fns[]],             seq( [k, g(fns[1])(k),             seq( f(fns[j])(k), j=2..7)], k=0..10)           ]         );   interface(displayprecision=-1):   interface(rtablesize=10):
 (3)
 > # # Generate plots from this second calculation #   plotopts2:=0..10, color=blue, thickness=4, linestyle=2,labeldirections = [horizontal, vertical]:   ylabels2:=[ "National Gross Income GI(t)",               "Gross National Savings G[S](t)",               "Gross Domestic Product G[D](t)",               "Public Expenditure on Education P[E](t)",               "Energy Consumption E[C](t)",               "Unemployment U(t)",               "Poverty P[R](t)"             ]:   a:= [ plot         ( g(fns[1]),           plotopts2,           labels=["Time in years", ylabels2[1]]         ),         seq          ( plot            ( f(fns[k]),              plotopts2,              labels = ["Time in years", ylabels2[k]]            ),            k = 2 .. 7          )        ]:   display(a[1]);   display(a[2]);   display(a[3]);   display(a[4]);   display(a[5]);   display(a[6]);   display(a[7]);
 > # # Combine displays from the two calculations. Note that # (if I believe the labels) only the first five entries # in each represent the *same* quantity. The last two # entries in the first calculation are # # "Number of Available Vacancies", # "Employed Individuals E(t)" # # whereas the last two entries in the second calculation # are # # "Unemployment U(t)", # "Poverty P[R](t)" # # So definitely not the same thing!! I assume therefore # it only makes sense to compare the first five dependent # variables in each calculation # # And because all of the variables are coupled this may # have some impact on the calculations themselves. No way # I can tell!! # # Anyhow, for what it is worth #   display( [v[1],a[1]] );   display( [v[2],a[2]] );   display( [v[3],a[3]] );   display( [v[4],a[4]] );   display( [v[5],a[5]] );
 >

## Modification...

In my previous post the final set of graphs rendered in a rather odd way when the x-coordinate is close to zero. I'm still not sure why this is happening. In the attached I have slightly modified the way the final calculations are done. The two methods provide the "same" answers, but in the attached the plot rendering is much better. (I'm still not sure why!)

Anyhow for what it is worth, check out the attached, which I still can't display on this site because someone has broken the "inline display" capability -again!

someODEs2.mw