tomleslie

5219 Reputation

15 Badges

9 years, 248 days

MaplePrimes Activity


These are replies submitted by tomleslie

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;

10

 

10

 

4

 

9

 

19

 

v__beta(r, theta, beta) = v__m(beta)*(1-r^2/R__max(theta, beta)^2)

 

(diff(v__r(r, theta, beta), r))*r*(19-r*cos(theta))+v__r(r, theta, beta)*(19-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;

R(beta) = 9+20*beta/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;

R__max(theta, beta) = (1/9)*R__max(theta, 0)*R(beta)

 

R__max(theta, 0) = 4*cos(theta)+(17+8*cos(2*theta))^(1/2)

(3)

eq38a:= subs(eq41, eq38);

R__max(theta, beta) = (1/9)*(4*cos(theta)+(17+8*cos(2*theta))^(1/2))*R(beta)

(4)

eq38b:=simplify(subs(eq39,eq38a));

R__max(theta, beta) = (1/9)*(9*Pi+20*beta)*(4*cos(theta)+(17+8*cos(2*theta))^(1/2))/Pi

(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));

R__max(theta, beta) = (1/36)*(18*Pi+40*beta)*(8*cos(theta)+2*(17+8*cos(2*theta))^(1/2))/Pi

(6)

eq44 := v__m(beta) = 1/(int(int((1-r^2/R__max(theta, beta)^2)*r, r = 0 .. R__max(theta, beta)), theta = 0 .. 2*Pi))

v__m(beta) = 4/(int(R__max(theta, beta)^2, theta = 0 .. 2*Pi))

(7)

eq44a:=simplify(expand(subs(eq38b,eq44)));

v__m(beta) = (162/25)*Pi/(9*Pi+20*beta)^2

(8)

eq47:=subs([eq44a,eq38b],eq43);

v__beta(r, theta, beta) = (162/25)*Pi*(1-81*r^2*Pi^2/((9*Pi+20*beta)^2*(4*cos(theta)+(17+8*cos(2*theta))^(1/2))^2))/(9*Pi+20*beta)^2

(9)

eq37a:=simplify(expand(subs(eq47, eq37)));

(1/5)*(-2361960*(cos(theta)*(9+16*cos(theta)^2)^(1/2)+4*cos(theta)^2+9/8)*(r*cos(theta)-19)*r*(Pi+(20/9)*beta)^5*(diff(v__r(r, theta, beta), r))-4723920*(cos(theta)*(9+16*cos(theta)^2)^(1/2)+4*cos(theta)^2+9/8)*(r*cos(theta)-19/2)*(Pi+(20/9)*beta)^5*v__r(r, theta, beta)-839808*(4*(Pi+(20/9)*beta)^2*cos(theta)^2+(9+16*cos(theta)^2)^(1/2)*(Pi+(20/9)*beta)^2*cos(theta)+(-(1/4)*r^2+9/8)*Pi^2+5*Pi*beta+(50/9)*beta^2)*Pi*r)/((9*Pi+20*beta)^5*(32*cos(theta)^2+8*cos(theta)*(9+16*cos(theta)^2)^(1/2)+9)) = 0

(10)

sol1:=simplify(dsolve(eq37a, v__r(r, theta,beta)), build);

v__r(r, theta, beta) = (Int(209952*(-2*(cos(2*theta)*r+r-38*cos(theta))*(Pi+(20/9)*beta)^2*(17+8*cos(2*theta))^(1/2)+152*(Pi+(20/9)*beta)^2*cos(2*theta)-4*(Pi+(20/9)*beta)^2*r*cos(3*theta)+(Pi^2*r^2-(33/2)*(Pi+(20/9)*beta)^2)*r*cos(theta)-19*Pi^2*r^2+(475/2)*(Pi+(20/9)*beta)^2)*Pi*r/(5*(9*Pi+20*beta)^5*(4*cos(2*theta)*r-152*cos(theta)+4*r)*(17+8*cos(2*theta))^(1/2)+5*(9*Pi+20*beta)^5*(33*r*cos(theta)+8*cos(3*theta)*r-304*cos(2*theta)-475)), r)+_F1(theta, beta))/(cos(theta)*r^2-19*r)

(11)

sol2:=eval(sol1, _F1(theta, beta)=0);

v__r(r, theta, beta) = (Int(209952*(-2*(cos(2*theta)*r+r-38*cos(theta))*(Pi+(20/9)*beta)^2*(17+8*cos(2*theta))^(1/2)+152*(Pi+(20/9)*beta)^2*cos(2*theta)-4*(Pi+(20/9)*beta)^2*r*cos(3*theta)+(Pi^2*r^2-(33/2)*(Pi+(20/9)*beta)^2)*r*cos(theta)-19*Pi^2*r^2+(475/2)*(Pi+(20/9)*beta)^2)*Pi*r/(5*(9*Pi+20*beta)^5*(4*cos(2*theta)*r-152*cos(theta)+4*r)*(17+8*cos(2*theta))^(1/2)+5*(9*Pi+20*beta)^5*(33*r*cos(theta)+8*cos(3*theta)*r-304*cos(2*theta)-475)), r))/(cos(theta)*r^2-19*r)

(12)

eq0 := o__1 = (1-2*L__1/d)^2*(1+2*L__1/d+4*c*beta/(d*Pi))^2; eq00 := o__2 = 1-4*L__1*(cos(theta)*sqrt(1-4*L__1^2*sin(theta)^2/d^2)-L__1*cos(theta)/d)/d

o__1 = (1/25)*(9/5+4*beta/Pi)^2

 

o__2 = 1-(8/25)*cos(theta)*(25-16*sin(theta)^2)^(1/2)+(16/25)*cos(theta)

(13)

eq48 := v__r(r, theta, beta) = 8*(1+2*L__1/d)^2*c*r*(o__1-4*(r/d)^2*o__2)/(d^2*Pi*(1-2*L__1/d)^2*(1+2*L__1/d+4*c*beta/(d*Pi))^5*(R__f/d-r*cos(theta)/d))

v__r(r, theta, beta) = (324/5)*r*(o__1-(1/25)*r^2*o__2)/(Pi*(9/5+4*beta/Pi)^5*(19/10-(1/10)*r*cos(theta)))

(14)

eq48a := simplify(expand(subs([eq0, eq00], eq48)));

v__r(r, theta, beta) = -25920*(Pi^2*cos(theta)*(9+16*cos(theta)^2)^(1/2)*r^2-2*Pi^2*cos(theta)*r^2+(-(25/8)*r^2+81/8)*Pi^2+45*Pi*beta+50*beta^2)*Pi^2*r/((9*Pi+20*beta)^5*(r*cos(theta)-19))

(15)

#
# Check Maple's solution and OP's proposed solution
#
  simplify(expand(subs( sol2, eq37a)));
  simplify(expand(subs( eq48a, eq37a)));

13271040*((1/4)*(Pi^3*cos(theta)^2*r^2-(1/2)*Pi^3*r^2*cos(theta)+(-(1/2)*r^2+81/64)*Pi^3+((45/8)*beta-81/1600)*Pi^2+((25/4)*beta^2-(9/40)*beta)*Pi-(1/4)*beta^2)*cos(theta)*(9+16*cos(theta)^2)^(1/2)+Pi^3*cos(theta)^4*r^2-(1/2)*Pi^3*cos(theta)^3*r^2+((81/64-(7/32)*r^2)*Pi^3+((45/8)*beta-81/1600)*Pi^2+((25/4)*beta^2-(9/40)*beta)*Pi-(1/4)*beta^2)*cos(theta)^2-(9/64)*Pi^3*r^2*cos(theta)+(729/2048-(225/1024)*r^2)*Pi^3+((405/256)*beta+(81/25600)*r^2-729/51200)*Pi^2+((225/128)*beta^2-(81/1280)*beta)*Pi-(9/128)*beta^2)*Pi*r/((9*Pi+20*beta)^5*(32*cos(theta)^2+8*cos(theta)*(9+16*cos(theta)^2)^(1/2)+9)) = 0

(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) ));

v__r((1/9)*(9*Pi+20*beta)*(4*cos(theta)+(17+8*cos(2*theta))^(1/2))/Pi, theta, beta) = 0

 

v__r(r, theta, beta) = (52488/5)*Pi*(-8*cos(theta)*(Pi+(20/9)*beta)^2*(17+8*cos(2*theta))^(1/2)-16*(Pi+(20/9)*beta)^2*cos(2*theta)+(r^2-25)*Pi^2-(1000/9)*Pi*beta-(10000/81)*beta^2)*r/((9*Pi+20*beta)^5*(16*cos(2*theta)+25+8*(17+8*cos(2*theta))^(1/2)*cos(theta))*(r*cos(theta)-19))

 

0 = 0

 

v__r((1/9)*(9*Pi+20*beta)*(4*cos(theta)+(17+8*cos(2*theta))^(1/2))/Pi, theta, beta) = 0

(17)

 


Download odesAgain2.mw

 

@pooyan1990 

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;

10

 

10

 

4

 

9

 

19

 

v__beta(r, theta, beta) = v__m(beta)*(1-r^2/R__max(theta, beta)^2)

 

(diff(v__r(r, theta, beta), r))*r*(19-r*cos(theta))+v__r(r, theta, beta)*(19-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;

R(beta) = 9+20*beta/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;

R__max(theta, beta) = (1/9)*R__max(theta, 0)*R(beta)

 

R__max(theta, 0) = 4*cos(theta)+(17+8*cos(2*theta))^(1/2)

(3)

eq38a:= subs(eq41, eq38);

R__max(theta, beta) = (1/9)*(4*cos(theta)+(17+8*cos(2*theta))^(1/2))*R(beta)

(4)

eq38b:=simplify(subs(eq39,eq38a));

R__max(theta, beta) = (1/9)*(9*Pi+20*beta)*(4*cos(theta)+(17+8*cos(2*theta))^(1/2))/Pi

(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));

R__max(theta, beta) = (1/36)*(18*Pi+40*beta)*(8*cos(theta)+2*(17+8*cos(2*theta))^(1/2))/Pi

(6)

eq44 := v__m(beta) = 1/(int(int((1-r^2/R__max(theta, beta)^2)*r, r = 0 .. R__max(theta, beta)), theta = 0 .. 2*Pi))

v__m(beta) = 4/(int(R__max(theta, beta)^2, theta = 0 .. 2*Pi))

(7)

eq44a:=simplify(expand(subs(eq38b,eq44)));

v__m(beta) = (162/25)*Pi/(9*Pi+20*beta)^2

(8)

eq47:=subs([eq44a,eq38b],eq43);

v__beta(r, theta, beta) = (162/25)*Pi*(1-81*r^2*Pi^2/((9*Pi+20*beta)^2*(4*cos(theta)+(17+8*cos(2*theta))^(1/2))^2))/(9*Pi+20*beta)^2

(9)

eq37a:=simplify(expand(subs(eq47, eq37)));

(1/5)*(-2361960*r*(r*cos(theta)-19)*(Pi+(20/9)*beta)^5*(cos(theta)*(9+16*cos(theta)^2)^(1/2)+4*cos(theta)^2+9/8)*(diff(v__r(r, theta, beta), r))-4723920*(Pi+(20/9)*beta)^5*(r*cos(theta)-19/2)*(cos(theta)*(9+16*cos(theta)^2)^(1/2)+4*cos(theta)^2+9/8)*v__r(r, theta, beta)-839808*r*(4*(Pi+(20/9)*beta)^2*cos(theta)^2+(9+16*cos(theta)^2)^(1/2)*(Pi+(20/9)*beta)^2*cos(theta)+(-(1/4)*r^2+9/8)*Pi^2+5*Pi*beta+(50/9)*beta^2)*Pi)/((9*Pi+20*beta)^5*(32*cos(theta)^2+8*cos(theta)*(9+16*cos(theta)^2)^(1/2)+9)) = 0

(10)

sol1:=simplify(dsolve(eq37a, v__r(r, theta,beta)), build);

v__r(r, theta, beta) = (Int(209952*(-2*(r*cos(2*theta)+r-38*cos(theta))*(Pi+(20/9)*beta)^2*(17+8*cos(2*theta))^(1/2)+152*(Pi+(20/9)*beta)^2*cos(2*theta)-4*r*(Pi+(20/9)*beta)^2*cos(3*theta)+(Pi^2*r^2-(33/2)*(Pi+(20/9)*beta)^2)*r*cos(theta)-19*Pi^2*r^2+(475/2)*(Pi+(20/9)*beta)^2)*r*Pi/(5*(9*Pi+20*beta)^5*(4*r*cos(2*theta)-152*cos(theta)+4*r)*(17+8*cos(2*theta))^(1/2)+5*(9*Pi+20*beta)^5*(33*r*cos(theta)+8*cos(3*theta)*r-304*cos(2*theta)-475)), r)+_F1(theta, beta))/(cos(theta)*r^2-19*r)

(11)

sol2:=eval(sol1, _F1(theta, beta)=0);

v__r(r, theta, beta) = (Int(209952*(-2*(r*cos(2*theta)+r-38*cos(theta))*(Pi+(20/9)*beta)^2*(17+8*cos(2*theta))^(1/2)+152*(Pi+(20/9)*beta)^2*cos(2*theta)-4*r*(Pi+(20/9)*beta)^2*cos(3*theta)+(Pi^2*r^2-(33/2)*(Pi+(20/9)*beta)^2)*r*cos(theta)-19*Pi^2*r^2+(475/2)*(Pi+(20/9)*beta)^2)*r*Pi/(5*(9*Pi+20*beta)^5*(4*r*cos(2*theta)-152*cos(theta)+4*r)*(17+8*cos(2*theta))^(1/2)+5*(9*Pi+20*beta)^5*(33*r*cos(theta)+8*cos(3*theta)*r-304*cos(2*theta)-475)), r))/(cos(theta)*r^2-19*r)

(12)

eq0 := o__1 = (1-2*L__1/d)^2*(1+2*L__1/d+4*c*beta/(d*Pi))^2; eq00 := o__2 = 1-4*L__1*(cos(theta)*sqrt(1-4*L__1^2*sin(theta)^2/d^2)-L__1*cos(theta)/d)/d

o__1 = (1/25)*(9/5+4*beta/Pi)^2

 

o__2 = 1-(8/25)*cos(theta)*(25-16*sin(theta)^2)^(1/2)+(16/25)*cos(theta)

(13)

eq48 := v__r(r, theta, beta) = 8*(1+2*L__1/d)^2*c*r*(o__1-4*(r/d)^2*o__2)/(d^2*Pi*(1-2*L__1/d)^2*(1+2*L__1/d+4*c*beta/(d*Pi))^5*(R__f/d-r*cos(theta)/d))

v__r(r, theta, beta) = (324/5)*r*(o__1-(1/25)*r^2*o__2)/(Pi*(9/5+4*beta/Pi)^5*(19/10-(1/10)*r*cos(theta)))

(14)

eq48a := simplify(expand(subs([eq0, eq00], eq48)))

v__r(r, theta, beta) = -25920*r*(Pi^2*cos(theta)*(9+16*cos(theta)^2)^(1/2)*r^2-2*Pi^2*cos(theta)*r^2+(-(25/8)*r^2+81/8)*Pi^2+45*Pi*beta+50*beta^2)*Pi^2/((9*Pi+20*beta)^5*(r*cos(theta)-19))

(15)

simplify(expand(subs( eq48a, eq37a)));

13271040*r*((1/4)*cos(theta)*(Pi^3*cos(theta)^2*r^2-(1/2)*Pi^3*r^2*cos(theta)+(-(1/2)*r^2+81/64)*Pi^3+((45/8)*beta-81/1600)*Pi^2+((25/4)*beta^2-(9/40)*beta)*Pi-(1/4)*beta^2)*(9+16*cos(theta)^2)^(1/2)+Pi^3*cos(theta)^4*r^2-(1/2)*Pi^3*cos(theta)^3*r^2+((-(7/32)*r^2+81/64)*Pi^3+((45/8)*beta-81/1600)*Pi^2+((25/4)*beta^2-(9/40)*beta)*Pi-(1/4)*beta^2)*cos(theta)^2-(9/64)*Pi^3*r^2*cos(theta)+(729/2048-(225/1024)*r^2)*Pi^3+((405/256)*beta+(81/25600)*r^2-729/51200)*Pi^2+((225/128)*beta^2-(81/1280)*beta)*Pi-(9/128)*beta^2)*Pi/((9*Pi+20*beta)^5*(32*cos(theta)^2+8*cos(theta)*(9+16*cos(theta)^2)^(1/2)+9)) = 0

(16)

simplify(expand(subs( sol2, eq37a)));

0 = 0

(17)

 


 

Download odesAgain.mw

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?

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

 

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:

`Standard Worksheet Interface, Maple 18.02, Windows 7, October 20 2014 Build ID 991181`

(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):

M := Matrix(12, 8, {(1, 1) = T, (1, 2) = J[1](T), (1, 3) = J[2](T), (1, 4) = J[3](T), (1, 5) = J[4](T), (1, 6) = J[5](T), (1, 7) = J[6](T), (1, 8) = J[7](T), (2, 1) = 0., (2, 2) = .5, (2, 3) = .3, (2, 4) = .21, (2, 5) = .14, (2, 6) = .7, (2, 7) = .45, (2, 8) = .14, (3, 1) = 1., (3, 2) = .999699523826566, (3, 3) = .255541352976897, (3, 4) = .21206859553329, (3, 5) = 0.898799904012718e-1, (3, 6) = .778114947832419, (3, 7) = 1.08216751409771, (3, 8) = .139302495316784, (4, 1) = 2., (4, 2) = 1.49919900731425, (4, 3) = .217671541999843, (4, 4) = .209656728291875, (4, 5) = 0.649490252974368e-1, (4, 6) = .830991859590075, (4, 7) = 1.71118146741969, (4, 8) = .138609081598146, (5, 1) = 3., (5, 2) = 1.998498523054, (5, 3) = .185414121860713, (5, 4) = .203953318724857, (5, 5) = 0.521558916668485e-1, (5, 6) = .871164327211628, (5, 7) = 2.33705762890861, (5, 8) = .137919697947623, (6, 1) = 4., (6, 2) = 2.49759814759995, (6, 3) = .157937353308399, (6, 4) = .195905077686603, (6, 5) = 0.451233125882032e-1, (6, 6) = .904634733867582, (6, 7) = 2.95981170927178, (6, 8) = .137234263359935, (7, 1) = 5., (7, 2) = 2.99649795956702, (7, 3) = .13453275778874, (7, 4) = .186260703845659, (7, 5) = 0.407659190040016e-1, (7, 6) = .934226826671212, (7, 7) = 3.57945935204539, (7, 8) = .136552686047415, (8, 1) = 6., (8, 2) = 3.49519803867417, (8, 3) = .114596851235069, (8, 4) = .175607485224099, (8, 5) = 0.376153958656721e-1, (8, 6) = .961236526196306, (8, 7) = 4.19601612726572, (8, 8) = .135874870157663, (9, 1) = 7., (9, 2) = 3.99369846525699, (9, 3) = 0.976155831384228e-1, (9, 4) = .164401567091675, (9, 5) = 0.349921862832695e-1, (9, 6) = .986240652200646, (9, 7) = 4.8094975271433, (9, 8) = .135200720486221, (10, 1) = 8., (10, 2) = 4.49199932001408, (10, 3) = 0.831510845320758e-1, (10, 4) = .152992930806212, (10, 5) = 0.325970591587624e-1, (10, 6) = 1.00948948775473, (10, 7) = 5.4199189632638, (10, 8) = .134530145660569, (11, 1) = 9., (11, 2) = 4.99010068387068, (11, 3) = 0.708303753168741e-1, (11, 4) = .141645981916411, (11, 5) = 0.30310421412715e-1, (11, 6) = 1.03109464256827, (11, 7) = 6.02729576493829, (11, 8) = .13386306017299, (12, 1) = 10., (12, 2) = 5.48800263790178, (12, 3) = 0.603357524172864e-1, (12, 4) = .130556463866947, (12, 5) = 0.280945993212309e-1, (12, 6) = 1.05111646607856, (12, 7) = 6.63164317840872, (12, 8) = .133199385555971})

(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):

Matrix(12, 8, {(1, 1) = T, (1, 2) = B, (1, 3) = C, (1, 4) = DD, (1, 5) = E, (1, 6) = F, (1, 7) = G, (1, 8) = H, (2, 1) = 0, (2, 2) = .50, (2, 3) = .30, (2, 4) = .21, (2, 5) = .14, (2, 6) = .70, (2, 7) = .45, (2, 8) = .14, (3, 1) = 1, (3, 2) = .9997987729, (3, 3) = .2686433294, (3, 4) = .1043456400, (3, 5) = 0.8306455626e-1, (3, 6) = .7764343230, (3, 7) = 1.770372716, (3, 8) = .1407005386, (4, 1) = 2, (4, 2) = 1.499597546, (4, 3) = .2405660980, (4, 4) = 0.8230565704e-1, (4, 5) = 0.4320908694e-1, (4, 6) = .8231900200, (4, 7) = 7.214709275, (4, 8) = .1414020975, (5, 1) = 3, (5, 2) = 1.999396319, (5, 3) = .2154323836, (5, 4) = .1956610723, (5, 5) = -0.651753148e-2, (5, 6) = .8597593744, (5, 7) = 23.76632611, (5, 8) = .1421015341, (6, 1) = 4, (6, 2) = 2.499195092, (6, 3) = .1929512830, (6, 4) = .620052072, (6, 5) = -.1112389777, (6, 6) = .9058423115, (6, 7) = 62.21819391, (6, 8) = .1427941256, (7, 1) = 5, (7, 2) = 2.998993864, (7, 3) = .1728769115, (7, 4) = 1.654978010, (7, 5) = -.3342514853, (7, 6) = .9813463960, (7, 7) = 137.1729376, (7, 8) = .1434735694, (8, 1) = 6, (8, 2) = 3.498792637, (8, 3) = .1550084044, (8, 4) = 3.723797405, (8, 5) = -.7570238435, (8, 6) = 1.106386836, (8, 7) = 267.0428363, (8, 8) = .1441319826, (9, 1) = 7, (9, 2) = 3.998591410, (9, 3) = .1391899156, (9, 4) = 7.373727938, (9, 5) = -1.479197396, (9, 6) = 1.301286479, (9, 7) = 474.0498237, (9, 8) = .1447599022, (10, 1) = 8, (10, 2) = 4.498390183, (10, 3) = .1253106180, (10, 4) = 13.27584647, (10, 5) = -2.618586042, (10, 6) = 1.586575815, (10, 7) = 784.2254871, (10, 8) = .1453462858, (11, 1) = 9, (11, 2) = 4.998188956, (11, 3) = .1133047038, (11, 4) = 22.22508901, (11, 5) = -4.311176236, (11, 6) = 1.982992970, (11, 7) = 1227.411068, (11, 8) = .1458785105, (12, 1) = 10, (12, 2) = 5.497987729, (12, 3) = .1031513842, (12, 4) = 35.14025076, (12, 5) = -6.711126987, (12, 6) = 2.511483720, (12, 7) = 1837.257464, (12, 8) = .1463423735})

(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]] );

 

 

 

 

 

 

Download someODEs4.mw

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

@David Sycamore 

Just use the download link in my previous post and run it!

@michalkvasnicka 

I think it depends on the contents of the pathdef.m file. I'm not sure if you can have more than one of these!

@David Sycamore 

Pretty trivial to do by stripping down my previous worksheet, as in the atached

  restart;
  interface(version);
#
# Define the upper limit for the odd numbers
# to be tested and the upper limit for the
# even addends
#
  Nodd:=1000:
  Neven:=10000:
#
# Procedure which does the work
#
  doWork:= proc( p::posint, N::posint )
                 local #
                       # Edited this to run in Maple 18
                       # whihc means it will probably run
                       # in every subsequent version
                       #
                       # f:=NumberTheory:-PrimeFactors(p),
                        f:=numtheory:-factorset(p),
                       j:
                 for j from 2 by 2 to N do
                     if   isprime~(f+~j)={true}
                     then #
                          # Return the lowest even number
                          #
                            return j;
                     fi;
                 od;
                 return #
                        # No even addend found: return a zero
                        #
                          0;
           end proc:

`Standard Worksheet Interface, Maple 2019.1, Windows 7, May 21 2019 Build ID 1399874`

(1)

#
# Just get a list of the even addends for
# the odd numbers 3,5,7.....Nodd. Only even
# addends up to Neven are checked
#
 [ seq( doWork(j, Neven), j=3..Nodd,2)]

[2, 2, 4, 2, 2, 4, 2, 2, 4, 4, 6, 2, 2, 2, 6, 2, 6, 4, 4, 2, 4, 2, 6, 4, 2, 6, 2, 4, 2, 6, 4, 6, 4, 8, 2, 6, 2, 6, 4, 2, 6, 2, 2, 8, 4, 10, 12, 4, 2, 2, 4, 0, 2, 4, 4, 14, 6, 4, 6, 2, 2, 2, 4, 4, 6, 4, 2, 2, 10, 14, 6, 2, 4, 2, 6, 2, 6, 6, 8, 6, 4, 2, 6, 4, 4, 6, 6, 2, 2, 10, 10, 6, 2, 4, 2, 4, 0, 2, 12, 4, 12, 2, 8, 12, 12, 2, 18, 6, 10, 6, 4, 2, 2, 4, 0, 6, 6, 4, 2, 10, 2, 6, 4, 14, 6, 6, 2, 6, 4, 2, 6, 6, 8, 2, 6, 4, 2, 4, 10, 2, 10, 0, 6, 2, 4, 14, 2, 2, 6, 4, 2, 6, 4, 4, 2, 4, 0, 14, 2, 2, 12, 6, 4, 6, 6, 4, 6, 10, 14, 6, 4, 8, 2, 4, 4, 6, 2, 0, 8, 4, 2, 6, 6, 2, 6, 6, 2, 18, 4, 4, 6, 6, 4, 8, 6, 8, 18, 4, 4, 8, 6, 2, 6, 10, 2, 12, 6, 10, 2, 10, 14, 2, 6, 0, 2, 6, 2, 18, 4, 4, 6, 8, 2, 8, 2, 16, 6, 4, 2, 2, 4, 0, 12, 4, 10, 18, 12, 8, 8, 4, 0, 6, 4, 4, 8, 2, 2, 12, 4, 14, 6, 2, 4, 12, 6, 4, 6, 6, 8, 2, 18, 0, 6, 6, 2, 6, 2, 2, 6, 6, 10, 18, 10, 10, 12, 4, 0, 6, 4, 2, 6, 14, 4, 2, 6, 2, 6, 10, 4, 6, 6, 0, 6, 10, 2, 6, 6, 28, 2, 6, 4, 2, 6, 0, 6, 4, 2, 2, 12, 8, 12, 2, 0, 6, 10, 16, 12, 4, 2, 2, 4, 0, 6, 2, 10, 6, 6, 10, 2, 12, 0, 12, 8, 4, 6, 4, 2, 6, 4, 2, 8, 2, 4, 6, 10, 0, 12, 2, 8, 8, 4, 14, 6, 10, 4, 6, 6, 2, 8, 4, 10, 2, 6, 2, 24, 6, 0, 6, 4, 4, 8, 2, 14, 6, 6, 20, 6, 4, 8, 8, 4, 2, 24, 4, 14, 14, 6, 4, 12, 2, 2, 6, 10, 8, 24, 6, 8, 12, 6, 8, 6, 6, 2, 2, 10, 10, 18, 4, 4, 2, 4, 2, 2, 10, 4, 6, 6, 10, 14, 2, 2, 6, 6, 10, 6, 4, 0, 2, 4, 0, 14, 6, 2, 18, 4, 4, 6, 4, 14, 2, 4, 2, 20, 4, 2, 12, 2, 0, 12, 6, 4, 12, 4, 2, 8, 6, 0, 6, 10, 4, 18, 6, 4, 8, 4, 2, 2, 4, 4, 6, 6, 0, 6, 6, 14, 14, 2, 2, 12, 6, 2, 6, 4, 0, 6, 10, 0, 6, 8, 4, 8, 2, 0, 18, 6, 16, 12, 12, 4]

(2)

 

Download prmeProb4.mw

@Kitonum 

one can achieve the same thing using LPSolve, by adding the option 'integervariables', as in

Optimization:-LPSolve(P.V, {seq(C.V <=~ R), seq(V >=~ 0)}, maximize, integervariables=[F,K]);

whiihc produces

[8500, [F = 50, K = 150]]

 

@student_md 

pdetest() outputs a list whose entries are the "result/discrepancy" are the outcomes of substituting the "obtained_solution" into the PDE system and BCs, ICs

You have a single PDE, two BCs and a single IC, so there are a total of four checks to be performed. The output will be a list of four values. If everything is perfect, the output will be a list of four zeros. In a PDE system with (say) four PDEs, six BCs and three ICs, the output list woukd have thirteen entries.

Any non-zero entry indicates that the provided solution does not fulfil the PDE/BCs/IC. The order in the output list is the same order as the functions are given to pdetest(). So, for your case

pdetest(solution, [PDE, BCs, IC])

  1. first list entry is essentially obtained by substituting "solution" into the "PDE", and working out the difference between the lhs() and rhs() of the result. This ought to be '0', if "solution" satisfies the PDE
  2. subsequent three entries are obtained by evaluating the "solution" under the same conditions as the BCs/ICs and again evaluating rhs()-lhs() for each of the generated equations

You can experiment with this behaviour by testing the "solution" against "slightly" different BCs/ICs than those which were used to produce it. In general this will produce a nonzero entry (showing the resulting discrepancy) at the appropriate place in the output list. You can even change the "solution" (I'd advise "slightly"). This will produce a non-zero first entry - the PDE is not satisfied, and generally one or more non-zeros entries for the BCs/ICs

Experimenting with solutions PDEs/BCs/ICs in this way can often assisy you if you are "guessing" a solution (whihc we have all done!). You can immediately see which PDEs/BCs/ICs are unsatisfied,  by how much, and adjust your "guess" appropriately

In the attached I have added an execution group which shows the effect of testing the obtained solution against "somewhat modified" BCs/ICs, as well as changing the obtained "solution" itself

restart:
with(PDEtools):
PDE :=  diff(y(x,t), t)-diff(y(x,t), x,x,t)-diff(y(x,t), x$2)+ diff(y(x,t), x)+y(x,t)*diff(y(x,t),x)=exp(-t)*(cos(x)-sin(x)+1/2*exp(-t)*sin(2*x));

diff(y(x, t), t)-(diff(diff(diff(y(x, t), t), x), x))-(diff(diff(y(x, t), x), x))+diff(y(x, t), x)+y(x, t)*(diff(y(x, t), x)) = exp(-t)*(cos(x)-sin(x)+(1/2)*exp(-t)*sin(2*x))

(1)

# Initial/boundary conditions
  BCs:=y(0,t) = 0, y(Pi,t)=0;
  ICs:=y(x,0) =sin(x) ;

y(0, t) = 0, y(Pi, t) = 0

 

y(x, 0) = sin(x)

(2)

pdsolve({PDE, BCs,ICs}, y(x,t), HINT=`*`);  # NULL. Maple does not find an exact solution

num_sol := pdsolve(PDE, {BCs,ICs}, numeric);
num_sol:-plot3d(x=0..1, t=0..1);

_m900620480

 

 

  exact_solution:=exp(-t)*sin(x);
  Test1:=pdetest(y(x,t)=exact_solution,[PDE, BCs,ICs]);

 
 

exp(-t)*sin(x)

 

[0, 0, 0, 0]

(3)

#
# "Adjust"  one of the BCs
#
  BCs2:=y(0,t) = 1, y(Pi,t)=0:
  Test2:=pdetest(y(x,t)=exact_solution,[PDE, BCs2,ICs]);
#
# "Adjust"  another one of the BCs
#
  BCs3:=y(0,t) = 0, y(Pi,t)=1:
  Test3:=pdetest(y(x,t)=exact_solution,[PDE, BCs3,ICs]);
#
# "Adjust"  the IC
#
  ICs2:=y(x,0) =cos(x):
  Test4:=pdetest(y(x,t)=exact_solution,[PDE, BCs,ICs2]);
#
# "Adjust" the exact solution a couple of different ways
#
  Test5:=pdetest(y(x,t)=2*exact_solution,[PDE, BCs,ICs]);
  Test6:=pdetest(y(x,t)=2+exact_solution,[PDE, BCs,ICs]);

[0, -1, 0, 0]

 

[0, 0, -1, 0]

 

[0, 0, 0, (-1/2+(1/2)*I)*exp(-I*x)+(-1/2-(1/2)*I)*exp(I*x)]

 

[((1/2)*I)*exp(-t+I*x)+(1/2)*exp(-t+I*x)-((1/2)*I)*exp(-t-I*x)+(1/2)*exp(-t-I*x)+(3/2)*exp(-2*t)*sin(2*x), 0, 0, sin(x)]

 

[2*exp(-t)*cos(x), 2, 2, 2]

(4)

 

Download testPDE.mw

 

@Carl Love 

You might need parse~(ImportMatrix() ), but using a matrix with pretty "random" symbolic entries, the following works for me.

(You will need to change the file path for your installation);

A:=Matrix([ [x^2, int(x*cos(x),x)],
            [piecewise( x<5, 2*x, 3*x/2), diff(y(x),x$2)*exp(y(x))=0]
          ]
         );
ExportMatrix( "C:/Users/TomLeslie/Desktop/testMat", A);

Matrix(2, 2, {(1, 1) = x^2, (1, 2) = cos(x)+x*sin(x), (2, 1) = piecewise(x < 5, 2*x, (3/2)*x), (2, 2) = (diff(diff(y(x), x), x))*exp(y(x)) = 0})

 

82

(1)

restart;
B:=parse~(ImportMatrix("C:/Users/TomLeslie/Desktop/testMat"));

Matrix(2, 2, {(1, 1) = x^2, (1, 2) = cos(x)+x*sin(x), (2, 1) = piecewise(x < 5, 2*x, (3/2)*x), (2, 2) = (diff(diff(y(x), x), x))*exp(y(x)) = 0})

(2)

 


Download mattest.mw

 

to debug this problem if you posted the offending worksheet using the big green up-arrow in the Mapleprimes toolbar

It is really easy to "work around" because (in the above worksheet)

member( 10, [entries(t1, nolist)]);

will work. But this really shouldn't be necessary

@brian bovril 

The Omega ratio is explained in

https://en.wikipedia.org/wiki/Omega_ratio

and one of the references in that page

Kapsos, Michalis; Zymler, Steve; Christofides, Nicos; Rustem, Berç (Summer 2014). "Optimizing the Omega Ratio using Linear Programming" (PDF). Journal of Computational Finance. 17 (4): 49–57.

gives a more comprensive description of the linear programming approach.

As for the data used in Samir Khan's worksheet, I have no more idea than you about exactly what is meant by the numbers in each column. If I had to guess, based on "real-word" plausibility, then I'd probably assume that the numbers were percentage return per month - purely on the basis that moving up and down a few % per month seems "realistic"

As for trying your own portfolio - I don't give financial advice!! :-)

1 2 3 4 5 6 7 Last Page 1 of 147