mmcdara

6279 Reputation

17 Badges

7 years, 325 days

MaplePrimes Activity


These are answers submitted by mmcdara

From the expresison of the Peng-Robinson EOS given here wiki here is a step-by-step procedure to find the result.
It mimics what you would do by hand, for what it worth...

Note that @dharr already gave you a faster way to get the desired result.

restart

Peng-Robinson_eos

eos := P = R*T/(v-b)-a*alpha/(v*(v+b)+b*(v-b))

P = R*T/(v-b)-a*alpha/(v*(v+b)+b*(v-b))

(1)

relation_6 := A = a*P/R^2/T^2;
relation_7 := B = b*P/R/T;
relation_8 := Z = P*v/(R*T);

A = a*P/(R^2*T^2)

 

B = b*P/(R*T)

 

Z = P*v/(R*T)

(2)

eos_1 := eval(eos, isolate(relation_6, a));

P = R*T/(v-b)-A*R^2*T^2*alpha/(P*(v*(v+b)+b*(v-b)))

(3)

eos_2 := eval(eos_1, isolate(relation_7, b));

P = R*T/(v-B*R*T/P)-A*R^2*T^2*alpha/(P*(v*(v+B*R*T/P)+B*R*T*(v-B*R*T/P)/P))

(4)

eos_3 := eval(eos_2, isolate(relation_8, P));

Z*R*T/v = R*T/(v-B*v/Z)-A*v*R*T*alpha/(Z*(v*(v+B*v/Z)+B*v*(v-B*v/Z)/Z))

(5)

eos_4 := numer(lhs(eos_3))*denom(rhs(eos_3)) - numer(rhs(eos_3))*denom(lhs(eos_3)) = 0

Z^2*R*T*v^2*(B-Z)*(B^2-2*B*Z-Z^2)-Z^2*R*T*v^2*(A*B*alpha-A*Z*alpha-B^2+2*B*Z+Z^2) = 0

(6)

eos_5 := collect(eos_4, Z)

R*T*v^2*Z^5+(B*R*T*v^2-R*T*v^2)*Z^4+(-3*R*T*v^2*B^2-R*T*v^2*(-A*alpha+2*B))*Z^3+(R*T*v^2*B^3-R*T*v^2*(A*B*alpha-B^2))*Z^2 = 0

(7)

eos_6 := factor(eos_5)

-Z^2*R*T*v^2*(A*B*alpha-A*Z*alpha-B^3+3*B^2*Z-B*Z^2-Z^3-B^2+2*B*Z+Z^2) = 0

(8)

# assuming Z^2*R*T*v^2 <> 0:

eos_7 := collect(op(-1, lhs(eos_6)), Z) = 0

-Z^3+(-B+1)*Z^2+(-A*alpha+3*B^2+2*B)*Z+A*alpha*B-B^3-B^2 = 0

(9)

 

Download EOS.mw

The two equations can be solved formally, independently the one from the other i,f we omit ic/bc conditions.
This leave four integration constants to fix.
Given your ic/bc conditions, only three of them can be found as a fonction of the fourth one which thus remains undefined.

To fix this last one you need one more ic/bc:
(I tried to keep the things clear wh-ile proceeding step by step, which can look lengthy and not very astute).

restart

``

a := Pi; 1; b := Pi; 1; lambda := 0.1e-1; 1; beta := 2.5; 1; x[1] := -1; 1; x[2] := 1; 1; y[1] := 1.5; 1; y[2] := 1.5; 1; alpha := 1; 1; Q[1] := 40; 1; Q[2] := 35; 1; n := 3

Pi

 

Pi

 

0.1e-1

 

2.5

 

-1

 

1

 

1.5

 

1.5

 

1

 

40

 

35

 

3

(1)

upsilon := (2*n-1)*Pi/(2*b);

5/2

(2)

EQ1 := diff(U[1, n](x), x, x)-upsilon^2*U[1, n](x) = -2*(int(Q[1]*Dirac(x-x[1])*Dirac(eta-y[1])*sin(upsilon*eta), eta = 0 .. b))/b

diff(diff(U[1, 3](x), x), x)-(25/4)*U[1, 3](x) = 14.55468946*Dirac(x+1.)

(3)

EQ2 := -(diff(U[2, n](x), x, x))-upsilon^2*U[2, n](x) = -2*(int(Q[2]*Dirac(x-x[2])*Dirac(eta-y[2])*sin(upsilon*eta), eta = 0 .. b))/b

-(diff(diff(U[2, 3](x), x), x))-(25/4)*U[2, 3](x) = 12.73535328*Dirac(x-1.)

(4)

 

bc := limit(U[2, n](x), x = infinity) < infinity, alpha*(diff(U[1, n](-a), x))-beta*U[1, n](-a) = 0, U[1, n](0) = U[2, n](0), (D(U[1, n]))(0) = lambda*(D(U[2, n]))(0)

limit(U[2, 3](x), x = infinity) < infinity, -2.5*U[1, 3](-Pi) = 0, U[1, 3](0) = U[2, 3](0), (D(U[1, 3]))(0) = 0.1e-1*(D(U[2, 3]))(0)

(5)

dsolve(EQ1);
sol1 := unapply(rhs(%), x);

sol2 := dsolve(EQ2);
sol2 := unapply(eval(rhs(sol2), [_C1=_C3, _C2=_C4]), x);

U[1, 3](x) = exp((5/2)*x)*_C2+exp(-(5/2)*x)*_C1+(727734473/250000000)*Heaviside(x+1)*(exp(5/2+(5/2)*x)-exp(-5/2-(5/2)*x))

 

proc (x) options operator, arrow; exp((5/2)*x)*_C2+exp(-(5/2)*x)*_C1+(727734473/250000000)*Heaviside(x+1)*(exp(5/2+(5/2)*x)-exp(-5/2-(5/2)*x)) end proc

 

U[2, 3](x) = sin((5/2)*x)*_C2+cos((5/2)*x)*_C1-(39797979/7812500)*Heaviside(x-1)*sin(-5/2+(5/2)*x)

 

proc (x) options operator, arrow; sin((5/2)*x)*_C4+cos((5/2)*x)*_C3-(39797979/7812500)*Heaviside(x-1)*sin(-5/2+(5/2)*x) end proc

(6)

BC1 := select(has, {bc}, U[1, n]);
BC2 := select(has, {bc}, U[2, n]);

{-2.5*U[1, 3](-Pi) = 0, U[1, 3](0) = U[2, 3](0), (D(U[1, 3]))(0) = 0.1e-1*(D(U[2, 3]))(0)}

 

{U[1, 3](0) = U[2, 3](0), (D(U[1, 3]))(0) = 0.1e-1*(D(U[2, 3]))(0), limit(U[2, 3](x), x = infinity) < infinity}

(7)

bc1 := remove(has, BC1, U[2, n]);
bc2 := remove(has, BC2, U[1, n]);
 

{-2.5*U[1, 3](-Pi) = 0}

 

{limit(U[2, 3](x), x = infinity) < infinity}

(8)

condition1 := eval(bc1, U[1, n](-Pi)=sol1(-Pi));
condition2 := eval(bc2, U[2, n](-Pi)=sol2(-Pi));

{-2.5*exp(-(5/2)*Pi)*_C2-2.5*exp((5/2)*Pi)*_C1 = 0}

 

{limit(U[2, 3](x), x = infinity) < infinity}

(9)

BC := (BC1 minus bc1) union (BC2 minus bc2);

condition34 :=  eval(BC, {U[1, n]=sol1, U[2, n]=sol2});

{U[1, 3](0) = U[2, 3](0), (D(U[1, 3]))(0) = 0.1e-1*(D(U[2, 3]))(0)}

 

{_C2+_C1+(727734473/250000000)*exp(5/2)-(727734473/250000000)*exp(-5/2) = _C3, (5/2)*_C2-(5/2)*_C1+(727734473/100000000)*exp(5/2)+(727734473/100000000)*exp(-5/2) = 0.2500000000e-1*_C4}

(10)

Csol := solve({condition1[], condition34[]}, {_C1, _C2, _C3})

{_C1 = 0.5380266007e-5-0.1507017048e-8*_C4, _C2 = -35.70142224+0.9999998493e-2*_C4, _C3 = -.4778779052+0.9999996986e-2*_C4}

(11)

collect(eval(sol2(x), Csol), _C4);

# All the functions having values in [-1, 1], their sum is clearli finite.
# Thus condition2 is always satisfied
 

(sin((5/2)*x)+0.9999996986e-2*cos((5/2)*x))*_C4-.4778779052*cos((5/2)*x)-5.094141312*Heaviside(x-1)*sin(-5/2+(5/2)*x)

(12)

SolC4free := unapply(eval([sol1(x), sol2(x)], Csol), x):
print~(['U'[1, n](x), 'U'[2, n](x)] =~ %(x)):

U[1, 3](x) = exp((5/2)*x)*(-35.70142224+0.9999998493e-2*_C4)+exp(-(5/2)*x)*(0.5380266007e-5-0.1507017048e-8*_C4)+(727734473/250000000)*Heaviside(x+1)*(exp(5/2+(5/2)*x)-exp(-5/2-(5/2)*x))

 

U[2, 3](x) = sin((5/2)*x)*_C4+cos((5/2)*x)*(-.4778779052+0.9999996986e-2*_C4)-(39797979/7812500)*Heaviside(x-1)*sin(-5/2+(5/2)*x)

(13)

# You need another relation to determine _C4

Download What_I_would_do.mw

I believe that when a command like

c1 := dsolve({b1, bcs}, numeric); 

returns an error, a good behavior is to print bcs and b1:

  • you will then see that bcs is unassigned (you wrote bcs instead of Bcs, or Bcs instead of bcs, it's uo to you to correct this),
  • and b1 (b2...) are not what you expect:
    Look to the mu and rho terms indexed bu f(x).
    This comes for you writting mu[f] and rho[f] (for instance instead of mu__f and rho__f)  and the consion (f=f(x)) that thiw writting induces.

There are several ways to fix this.
The simpler, IMO, is to do this (here for a single equation... I think you transpose yourself to the others)

restart

with(plots):

eq1 := mu[hnf]*rho[f]*(diff(F(x), x, x, x, x))/(rho[hnf]*mu[f])+3*alpha*(diff(F(x), x, x))+alpha*eta*(diff(F(x), x, x, x))-2*R*F(x)*(diff(F(x), x, x, x))-rho[f]*M*(diff(F(x), x, x))/rho[hnf] = 0

997.1*(diff(diff(diff(diff(F(x), x), x), x), x))/((1-phi[1])^2.5*(1-phi[2])^2.5*(1-phi[2](997.1+5322.9*phi[1])+3970*phi[2]))+3*alpha*(diff(diff(F(x), x), x))+alpha*eta*(diff(diff(diff(F(x), x), x), x))-2*R*F(x)*(diff(diff(diff(F(x), x), x), x))-997.1*M*(diff(diff(F(x), x), x))/(1-phi[2](997.1+5322.9*phi[1])+3970*phi[2]) = 0

(1)

a1 := [phi[1] = 0.1e-1, phi[2] = 0.1e-1, R = 1.0, M = 1, alpha = -1, eta = 1]

[phi[1] = 0.1e-1, phi[2] = 0.1e-1, R = 1.0, M = 1, alpha = -1, eta = 1]

(2)

# You forgot assigning mu[f]
# I choosed:
mu[f] := 1;

rho[f]  := 997.1:
rho[s1] := 6320:
rho[s2] := 3970:
c[p][f] := 4180:
c[p][s1] := 531.5:
c[p][s2] := 765:
k[s1] := 76.5:
k[s2] := 40:
k[f] := .613:
k[hnf] := k[nf]*(k[s2]+2*k[nf]-2*phi[2]*(k[nf]-k[s2]))/(k[s2]+2*k[nf]+phi[2]*(k[nf]-k[s2])):
k[nf] := k[f]*(k[s1]+2*k[f]-2*phi[1]*(k[f]-k[s1]))/(k[s1]+2*k[f]+phi[1]*(k[f]-k[s1])):
mu[hnf] := mu[f]/((1-phi[1])^2.5*(1-phi[2])^2.5):
rho[hnf] := (1-phi[2])((1-phi[1])*rho[f]+phi[1]*rho[s1])+phi[2]*rho[s2]:

1

(3)

b1 :=eval(eq1, a1)

25.76766428*(diff(diff(diff(diff(F(x), x), x), x), x))-27.50479233*(diff(diff(F(x), x), x))-(diff(diff(diff(F(x), x), x), x))-2.0*F(x)*(diff(diff(diff(F(x), x), x), x)) = 0

(4)

c1 := dsolve({b1, F(-1) = -1, F(1) = 1, (D(F))(-1) = 0, (D(F))(1) = 0}, numeric);
c1(0)

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(11, {(1) = -1.0, (2) = -.8300102949676604, (3) = -.64658629648752, (4) = -.44482312709262306, (5) = -.22232215574426015, (6) = 0.15250483747107274e-1, (7) = .24946665186051725, (8) = .4639899121532438, (9) = .6564323348024983, (10) = .8322625498127126, (11) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(11, 4, {(1, 1) = -1.0, (1, 2) = .0, (1, 3) = 3.192175801161809, (1, 4) = -4.282041992577323, (2, 1) = -.9572654340777135, (2, 2) = .4835342091017062, (2, 3) = 2.5121494492784278, (2, 4) = -3.7402606312251065, (3, 1) = -.8300352541330445, (3, 2) = .8841027104724593, (3, 3) = 1.8692043139533685, (3, 4) = -3.292670479766869, (4, 1) = -.6179839126514046, (4, 2) = 1.1968100545740237, (4, 3) = 1.242136735575425, (4, 4) = -2.94762809896349, (5, 1) = -.3262342001087303, (5, 2) = 1.4023621621988174, (5, 3) = .613603778343822, (5, 4) = -2.729301994490605, (6, 1) = 0.1820734167649253e-1, (6, 2) = 1.4720751087907518, (6, 3) = -0.24456516682766806e-1, (6, 4) = -2.6719253988214358, (7, 1) = .35656320693026283, (7, 2) = 1.392374561630252, (7, 3) = -.6607611929418454, (7, 4) = -2.791085177884434, (8, 1) = .6353718392299886, (8, 2) = 1.184636970385788, (8, 3) = -1.28563165989483, (8, 4) = -3.061416614334482, (9, 1) = .8358016089843219, (9, 2) = .8783720658960743, (9, 3) = -1.9096996176864884, (9, 4) = -3.448512042797582, (10, 1) = .9574998312478541, (10, 2) = .4869409158526661, (10, 3) = -2.5569793740766595, (10, 4) = -3.936854358364328, (11, 1) = 1.0, (11, 2) = .0, (11, 3) = -3.265818370018978, (11, 4) = -4.538381846985902}, datatype = float[8], order = C_order); YP := Matrix(11, 4, {(1, 1) = .0, (1, 2) = 3.192175801161809, (1, 3) = -4.282041992577323, (1, 4) = 3.573555347655438, (2, 1) = .4835342091017062, (2, 2) = 2.5121494492784278, (2, 3) = -3.7402606312251065, (2, 4) = 2.81425324074925, (3, 1) = .8841027104724593, (3, 2) = 1.8692043139533685, (3, 3) = -3.292670479766869, (3, 4) = 2.079562608876083, (4, 1) = 1.1968100545740237, (4, 2) = 1.242136735575425, (4, 3) = -2.94762809896349, (4, 4) = 1.3528683846147953, (5, 1) = 1.4023621621988174, (5, 2) = .613603778343822, (5, 3) = -2.729301994490605, (5, 4) = .6181594743990174, (6, 1) = 1.4720751087907518, (6, 2) = -0.24456516682766806e-1, (6, 3) = -2.6719253988214358, (6, 4) = -.1335741606703747, (7, 1) = 1.392374561630252, (7, 2) = -.6607611929418454, (7, 3) = -2.791085177884434, (7, 4) = -.8908677512914607, (8, 1) = 1.184636970385788, (8, 2) = -1.28563165989483, (8, 3) = -3.061416614334482, (8, 4) = -1.6420861348793063, (9, 1) = .8783720658960743, (9, 2) = -1.9096996176864884, (9, 3) = -3.448512042797582, (9, 4) = -2.3959853946020426, (10, 1) = .4869409158526661, (10, 2) = -2.5569793740766595, (10, 3) = -3.936854358364328, (10, 4) = -3.1747198703412507, (11, 1) = .0, (11, 2) = -3.265818370018978, (11, 3) = -4.538381846985902, (11, 4) = -4.014364688696904}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(11, {(1) = -1.0, (2) = -.8300102949676604, (3) = -.64658629648752, (4) = -.44482312709262306, (5) = -.22232215574426015, (6) = 0.15250483747107274e-1, (7) = .24946665186051725, (8) = .4639899121532438, (9) = .6564323348024983, (10) = .8322625498127126, (11) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(11, 4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = -0.783523289563683e-7, (1, 4) = 0.7366704178590536e-7, (2, 1) = 0.2911481441576516e-8, (2, 2) = -0.15958984565276724e-7, (2, 3) = -0.6384361270402364e-7, (2, 4) = 0.58046226638656524e-7, (3, 1) = 0.36654860061536426e-8, (3, 2) = -0.3100302176659618e-7, (3, 3) = -0.4977384154160686e-7, (3, 4) = 0.4211149953956398e-7, (4, 1) = 0.26894747375389394e-8, (4, 2) = -0.4507455567166456e-7, (4, 3) = -0.35491956774726486e-7, (4, 4) = 0.26450286379618158e-7, (5, 1) = 0.10170405390149387e-8, (5, 2) = -0.5655691698445183e-7, (5, 3) = -0.19647113910452217e-7, (5, 4) = 0.13392577583337413e-7, (6, 1) = 0.8651628438285204e-10, (6, 2) = -0.6073535104451812e-7, (6, 3) = -0.1031332037431753e-8, (6, 4) = 0.8970963692462515e-8, (7, 1) = -0.9899664033687283e-9, (7, 2) = -0.54747329531663056e-7, (7, 3) = 0.1749707327574523e-7, (7, 4) = 0.16780034100927893e-7, (8, 1) = -0.27942941245322632e-8, (8, 2) = -0.4330823755475724e-7, (8, 3) = 0.3338267394836564e-7, (8, 4) = 0.30968772526141964e-7, (9, 1) = -0.3882953260480258e-8, (9, 2) = -0.3019873197276184e-7, (9, 3) = 0.4787482229906417e-7, (9, 4) = 0.46718003309339046e-7, (10, 1) = -0.3217412416134993e-8, (10, 2) = -0.16015748539050472e-7, (10, 3) = 0.6247946054645264e-7, (10, 4) = 0.6295704686039219e-7, (11, 1) = .0, (11, 2) = .0, (11, 3) = 0.7832429590649342e-7, (11, 4) = 0.7978731193395113e-7}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[11] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(7.978731193395113e-8) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [4, 11, [F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[11] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[11] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(4, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(11, 4, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(4, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(11, 4, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)]'[i] = yout[i], i = 1 .. 4)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[11] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(7.978731193395113e-8) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [4, 11, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[11] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[11] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(11, 4, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(4, {(1) = 0., (2) = 0., (3) = 0., (4) = 0.}); `dsolve/numeric/hermite`(11, 4, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 4)] end proc, (2) = Array(0..0, {}), (3) = [x, F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[F(x), diff(F(x), x), diff(diff(F(x), x), x), diff(diff(diff(F(x), x), x), x)]'[i] = res[i+1], i = 1 .. 4)] catch: error  end try end proc

 

[x = 0., F(x) = HFloat(-0.0042437804179444585), diff(F(x), x) = HFloat(1.4721374397210758), diff(diff(F(x), x), x) = HFloat(0.016277974284338716), diff(diff(diff(F(x), x), x), x) = HFloat(-2.6702559502576326)]

(5)

``

Download A_correct_syntax.mw

restart

T := (1+(0.1388888889e-3*(1-y))*(1+y)*(15.57859618+2.919247663*beta-6.091152226*y^2+0.2836741238e-2*y+0.8630974840e-2*y^2*(0.6026785714e-1*y^2+.5*beta)-(.74322432*(0.4041321429e-2*y^2-.1676400000*beta-1.546875000))*y+0.5201703585e-3*y^6+0.4315487420e-2*y^4*beta+0.9797142858e-1*y^2*(-2.48236004+2.28352528*beta)-3*y*(-0.3114853125e-1*y^2*(0.8035714286e-1*y^2+.5*beta)+.3135477600*y^2+.3229905171-.2691297675*beta)+(0.4683861600e-1*(1.012500000*y^2+4.5*beta))*y^2-(.3025000000*(-.4526280000*y^2-2.514600000*beta+33.75000000))*y))*(-(0.1388888889e-1*(1-y))*(-9.223735354+.1058510135*y^2+0.28316841e-2*beta+.1420729992*y-0.2831684659e-1*y^2*(-0.7500000000e-1*y^2+.5*beta)+(.6096*(0.3771900000e-1*y^2-.1676400000*beta-4.5))*y)+(0.1388888889e-1*(1+y))*(-9.223735354+.1058510135*y^2+0.28316841e-2*beta+.1420729992*y-0.2831684659e-1*y^2*(-0.7500000000e-1*y^2+.5*beta)+(.6096*(0.3771900000e-1*y^2-.1676400000*beta-4.5))*y)-(0.1388888889e-1*(1+y))*(1-y)*(.2117020270*y-2.601127001-0.5663369318e-1*y*(-0.7500000000e-1*y^2+.5*beta)+0.4247526988e-2*y^3+0.6898050720e-1*y^2-.1021933440*beta)-.3048000000*y-.5500000000)-(0.3472222222e-2*(1-y))*(-9.223735354+.1058510135*y^2+0.28316841e-2*beta+.1420729992*y-0.2831684659e-1*y^2*(-0.7500000000e-1*y^2+.5*beta)+(.6096*(0.3771900000e-1*y^2-.1676400000*beta-4.5))*y)+(0.3472222222e-2*(1+y))*(-9.223735354+.1058510135*y^2+0.28316841e-2*beta+.1420729992*y-0.2831684659e-1*y^2*(-0.7500000000e-1*y^2+.5*beta)+(.6096*(0.3771900000e-1*y^2-.1676400000*beta-4.5))*y)-(0.3472222222e-2*(1+y))*(1-y)*(.2117020270*y-2.601127001-0.5663369318e-1*y*(-0.7500000000e-1*y^2+.5*beta)+0.4247526988e-2*y^3+0.6898050720e-1*y^2-.1021933440*beta)-0.7620000000e-1*y-.1375000000-0.8333333333e-1*beta*(-(0.1388888889e-1*(1-y))*(-9.223735354+.1058510135*y^2+0.28316841e-2*beta+.1420729992*y-0.2831684659e-1*y^2*(-0.7500000000e-1*y^2+.5*beta)+(.6096*(0.3771900000e-1*y^2-.1676400000*beta-4.5))*y)+(0.1388888889e-1*(1+y))*(-9.223735354+.1058510135*y^2+0.28316841e-2*beta+.1420729992*y-0.2831684659e-1*y^2*(-0.7500000000e-1*y^2+.5*beta)+(.6096*(0.3771900000e-1*y^2-.1676400000*beta-4.5))*y)-(0.1388888889e-1*(1+y))*(1-y)*(.2117020270*y-2.601127001-0.5663369318e-1*y*(-0.7500000000e-1*y^2+.5*beta)+0.4247526988e-2*y^3+0.6898050720e-1*y^2-.1021933440*beta)-.3048000000*y-.5500000000)^3:

with(plots):

kernelopts(version):

minimum := eval(T, [y=0.5, beta=-3]):

c := contourplot(T, beta = -3 .. 3, y = -.5 .. .5, filledregions=true, coloring=["Blue", "Yellow"]):
p := plot3d(T, beta = -3 .. 3, y = -.5 .. .5, style = surfacecontour, colorscheme=["Blue", "Yellow"]):
phi := plottools:-transform((x, y) -> [x, y, minimum*1.2]):
display(p, phi(c), orientation=[130, 73, -34])

 

maximum := eval(T, [y=-0.5, beta=3]):
cont    := [seq(minimum..maximum, (maximum-minimum)/10)]:

c := contourplot(T, beta = -3 .. 3, y = -.5 .. .5, contours=cont, color=gray):
p := plot3d(T, beta = -3 .. 3, y = -.5 .. .5, style = surfacecontour, contours=cont, colorscheme=["Blue", "Green", "Yellow", "Red"]):
phi := plottools:-transform((x, y) -> [x, y, minimum*1.2]):
display(p, phi(p), phi(c), orientation=[130, 73, -34])

 

# this corresponds the attached picture you provide in your question

p := plot3d(T, beta = -3 .. 3, y = -.5 .. .5, style = patch, colorscheme=["Blue", "Green", "Yellow", "Red"]):
phi := plottools:-transform((x, y) -> [x, y, minimum*1.2]):
display(p, phi(p), orientation=[130, 73, -34])

 

 

Download At_least_two_ways.mw

Here is a solution:

One_way.mw (UPDATED)

... here is something that could help you.
Formal and numerical resolutions are proposed, the latter approach is done in two different ways.
sol_1.mw

Did you try ImageTools?

restart:

with(ImageTools):

# Source of the original image: Maple_leaves
# Downloaded in directory f with name Maple.jpg

f := cat(currentdir(), "/Desktop/Maple.jpg"):
img := Read(f):
Embed(img):

Embed(Rotate(img, 90))

Embed(Scale(img, 1/3, 2))

Embed(Flip(img, vertical))

 

Download Images.mw

The four images:

Statistics:-Trim seems to do what you want:

restart

with(Statistics):

N  := 10^5:
X  := Sample(Uniform(0, 1), N):

t0 := time():

k  := 1000:          # 1000 elements to select
# the 1000 smaller
CodeTools:-Usage( Trim(X, 0, 100*k/N) ):

memory used=2.87MiB, alloc change=0 bytes, cpu time=10.00ms, real time=126.00ms, gc time=0ns

 

randomize():

X  := Sample(Uniform(0, 1), N):
# the sorted 1000 smaller
CodeTools:-Usage( sort(Trim(X, 0, 100*k/N)) ):
 

memory used=1.55MiB, alloc change=1.53MiB, cpu time=2.00ms, real time=2.00ms, gc time=0ns

 

# Using sort:
randomize():
X  := Sample(Uniform(0, 1), N):
CodeTools:-Usage( sort(X)[1..k] ):

memory used=0.85MiB, alloc change=7.33MiB, cpu time=16.00ms, real time=5.00ms, gc time=0ns

 

alias(ST = StringTools):

N := 10^5:
K := 10^3:
text := (ST:-Explode@ST:-Random)(N, 'alnum'):numelems(text);

CodeTools:-Usage( ST:-Char~(Trim(ST:-Ord~(text), 0, 100*K/N)) ):

100000

 

memory used=9.22MiB, alloc change=0 bytes, cpu time=64.00ms, real time=64.00ms, gc time=0ns

 

 

Download PartialSorting.mw

@Sky-Bj 

varphi_c_c and var_phi_f not defined (see @Axel Vogt )
See here

restart

V := m^4*(varphi/M)^p:

V1 := diff(V, varphi):

V2 := diff(V1, varphi):

f := Zeta * (varphi^2):

f1 := diff(f, varphi):

f2 := diff(f1, varphi):

R:= simplify((V/3-f1*V1/3*V)/((1-kappa^2*f)/12*kappa^2+f1/V)):

ToInt := 3*V1*kappa^2*((2*V*V1)/3 - f1^2*V1*R/(3*V) - f1*V1^2/(3*V))/(V*(-f*kappa^2 + 1)*(-R*f1 - 2*V1)):
FromTo := varphi=varphi__c..varphi__f:   # Use your own range

data := [p=2, Zeta=1/6,M = 1, m = 1, kappa=1, varphi__c=0, varphi__f=1]: # Use your own range
meth := method=_DEFAULT: # Change it

evalf(eval(Int(ToInt, FromTo, meth), data));

Float(infinity)

(1)

limit(eval(ToInt, data), varphi=0, right);
plot(eval(ToInt, data), eval(FromTo, data))

infinity

 

 

data := [p=2, Zeta=1/6,M = 1, m = 1, kappa=1, varphi__c=1e-2, varphi__f=1]:
evalf(eval(Int(ToInt, FromTo, meth), data));

2.063819612

(2)

NULL

Download How_to_integrate.mw

@Vortex 

... for one can find the exact solution of your initial problem.
This is not that immediate, but it is possible, at least after a few trials and failures and with the help of an intermediate plot:
dsolve_exact.mw

solution := tanh((1/4)*sqrt(2)*x+(15/4)*sqrt(2))

The exact solution of your "new" probem (z(0)=0, z(15)=1) is here
dsolve_exact_new_problem.mw

solution := tanh((1/4)*sqrt(2)*x)

NewtonsMethod's help page:
The NewtonsMethod(f(x), x=a) command returns the result of applying 5 iterations of Newton's method for approximating a root.
So do not expect this command gives a converged result, it has to be used mainly to illutrate graphically how the Newton's method works (more: the Newton's method not always converges).

Bisection's method:
As f:=sin(x)+1 is never negative this method cannot be applied (it requires that the target function has opposite signs at the bounds of the initial interval.

A correct example

restart

f := sin(x)+.95; 1; plot(f, x = -3 .. 2)

sin(x)+.95

 

 

sols := fsolve(f, x = -2 .. 2)

-1.253235898

(1)

with(Student[Calculus1]):

NewtonsMethod(f, x = .5)

-1.253235897

(2)

with(Student[NumericalAnalysis]):

omega := [-1.5, 0]:

-1.253235897

(3)

 


Download help_roots_corrected.mw

ListTools:-Search(96, R);

Without even using maplemint one can see a lot of errors in your procedures.

In the attached (yellow highlighted red text) file some of these errors are fixed but it remains a lot of ambiguities.
Having said this, I have the feeling that what you are trying to achieve could be done more simply:

  • Procedure CC seems to be the discrtization of a first order ode with a forward Euler scheme?
  • Procedure normalisatie could be written in a single line.
  • Procedure integraal is a numeric integral with trapezoid rule (am I right?)
  • Finally optimalisatie seems to be a procedure to find the 0 of some function with a dichotomy strategy (am I right?)

Perhaps you should simply explain what your problem is so that you can receive more appropriate assistance?

Anyway here is my analysis of youy file

restart

C__0 := 1000*10^(-9):
Q := .8*((1/6)*10^(-7)):
R := 289.1*10^(-6):
V__r := 8*10^(-6):
m__b := 1.5*10^(-3):
rho := 290.:

Gamma__i := 3.01*10^(-6)*((1/6)*10^(-8)):

T__ex := [0, 30, 60, 120, 240, 360]*~60:
Cl__ex := [1, .56, .45, .32, .18, .11]:

CC := proc (C__0, Q, V__r, m__b, rho, R, Gamma__i)
  local t, dC, C__i, i;
  t := Vector(360,i -> 60*i);
  C__i := Vector(360, 0);
  C__i[1] := C__0;
  for i to numelems(t)-1 do
    dC := -Q*(1-exp(-3*m__b*sqrt(Gamma__i/(rho*t[i+1]))/(sqrt(Pi)*Q*R)))*C__i[i]/V__r;
    C__i[i+1] := C__i[i]+60*dC
  end do;
  return t, C__i
end proc:

t, C := CC(C__0, Q, V__r, m__b, rho, R, Gamma__i)

t, C := Vector(4, {(1) = ` 1 .. 360 `*Vector[column], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order}), Vector(4, {(1) = ` 1 .. 360 `*Vector[column], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

(1)

normalisatie := proc (C__i, C__0)
  local Z, i;
  Z := Vector(360, 0);
  for i to numelems(Z) do
    Z[i] := C__i[i]/C__0
  end do;
  return Z
end proc:


integraal := proc (t, r)
  local Integraal, i; # Change integraal into Integraal to avoid conflicts with the name of the procedure
  Integraal := 0;
  for i to numelems(t)-1 do
    Integraal := Integraal+((1/2)*r[i+1]+(1/2)*r[i])*(t[i+1]-t[i])
  end do:
  return Integraal
end proc:


optimalisatie := proc (t__ex, c__ex, C__0, Q, V__r, m__b, rho, R, Gamma__i_min, Gamma__i_max)
  local Gamma__i, t__sim, C__i, int__sim, int__ex, C__i_min, C__i_max, C__i0_min, C__i0_max, C__i0, int__sim_min, int__sim_max, i;
  local gamma__i_min, gamma__i_max;   # As gamma__i_min and Gamma__i_max are arguments of optimalisatiea
                                      # they cannot be assigned within this procedure.
  gamma__i_min := copy(Gamma__i_min): # what I understand
  gamma__i_max := copy(Gamma__i_max): # what I understand

  Gamma__i     := (1/2)*gamma__i_min+(1/2)*gamma__i_max;
  t__sim, C__i := CC(C__0, Q, V__r, m__b, rho, R, Gamma__i);
  int__sim     := integraal(t__sim, C__i);
  int__ex      := integraal(t__ex, c__ex);



if true then
  i := 0:  # You probably missed to initialize i
  while i < 10 and 1/10 < abs(int__sim-int__ex) do
    Gamma__i := (1/2)*gamma__i_min+(1/2)*gamma__i_max;
    t__sim, C__i_min := CC(C__0, Q, V__r, m__b, rho, R, gamma__i_min);
#print(`-----`, 'C__i_min'=C__i_min[1..5]^+);
    t__sim, C__i_max := CC(C__0, Q, V__r, m__b, rho, R, gamma__i_max);
#print(`-----`, 'C__i_min'=C__i_max[1..5]^+);
    t__sim, C__i     := CC(C__0, Q, V__r, m__b, rho, R, Gamma__i);
    C__i0_min        := normalisatie(C__i_min, C__0);
    C__i0_max        := normalisatie(C__i_max, C__0);
    C__i0            := normalisatie(C__i, C__0);
    int__sim_min     := integraal(t__sim, C__i0_min);
    int__sim_max     := integraal(t__sim, C__i0_max);
    int__sim         := integraal(t__sim, C__i0);


print('(int__sim-int__ex)*(int__sim_max-int__ex)' = (int__sim-int__ex)*(int__sim_max-int__ex));
    if 0 < (int__sim-int__ex)*(int__sim_max-int__ex) then
      gamma__i_max := Gamma__i;
    else
      gamma__i_min := Gamma__i;
    end if:
    i := i+1:
  end do;
end if:

  return Gamma__i
end proc:

Gamma__i_min := 0:
Gamma__i_max := 10: # You missed to initialize this variable



# Note that (int__sim-int__ex)*(int__sim_max-int__ex) keeps a constant value
# all along the computations

Gamma__i := optimalisatie(T__ex, Cl__ex, C__0, Q, V__r, m__b, rho, R, Gamma__i_min, Gamma__i_max)

(int__sim-int__ex)*(int__sim_max-int__ex) = 35676729.00

 

(int__sim-int__ex)*(int__sim_max-int__ex) = 35676729.00

 

(int__sim-int__ex)*(int__sim_max-int__ex) = 35676729.00

 

(int__sim-int__ex)*(int__sim_max-int__ex) = 35676729.00

 

(int__sim-int__ex)*(int__sim_max-int__ex) = 35676729.00

 

(int__sim-int__ex)*(int__sim_max-int__ex) = 35676729.00

 

(int__sim-int__ex)*(int__sim_max-int__ex) = 35676729.00

 

(int__sim-int__ex)*(int__sim_max-int__ex) = 35676729.00

 

(int__sim-int__ex)*(int__sim_max-int__ex) = 35676729.00

 

(int__sim-int__ex)*(int__sim_max-int__ex) = 35676729.00

 

5/512

(2)

 

Download A_few_corrections.mw

It is likely that a lot of people will provide you different methods to tackle your problem.
In the meantime, here are some answers to your question:

  • Method 1: to use in the case S contains a two vectors, each of the form (numeric)*x + (numeric)), returns the coefficients of f in the basis S (in case S is not abasis look at the end of the attached file).
  • Method 2: same restriction as Method 1, checks if S is a basis.
  • Method 3: same restriction and same purpose as Method 2, using LinearAlgebra:-Basis a can be easily extended to a set of N N-dimensional vectors.

restart:

S := {x+4, 3*x-7}:
f := -5*x + 10:

# Method 1: (constructive method which gives the coefficients

try
  solve(identity(a[1]*S[1]+a[2]*S[2]=f, x), [a[1], a[2]])[];
  check = eval(a[1]*S[1]+a[2]*S[2] - f, %);
catch:
  printf("The vectors in set S are not independent")
end try;

[a[1] = -5/19, a[2] = -30/19]

 

check = 0

(1)

# Method 2:
# Let's say U is the Unity "1"

SU := map(s -> add([selectremove(has, s, x)] *~ [1, U]), S);

# If SU = 0 iif x=0 and U=0 the two vectors of S are independent, so
# f can be rewritten as a linear combination of both of them.

solve(SU);

{x+4*U, 3*x-7*U}

 

{U = 0, x = 0}

(2)

# Method 3:
# Use of LinearAlgebra:-Basis: if the output of this function contains as many
# vectors has S, then the vectors of S form are linearly independent and
# f can be rewritten as a linear combination of both of them.

SU := map(s -> <[selectremove(has, s, x)]>, S);
LinearAlgebra:-Basis([SU[]]);

is(numelems(S)=numelems(%));

SU := {Vector(2, {(1) = x, (2) = 4}), Vector(2, {(1) = 3*x, (2) = -7})}

 

[Vector(2, {(1) = x, (2) = 4}), Vector(2, {(1) = 3*x, (2) = -7})]

 

true

(3)

# Method 3: a counter example

SC := {x+4, 2*x+8}:
SU := map(s -> <[selectremove(has, s, x)]>, SC);
LinearAlgebra:-Basis([SU[]]);

# Note that Method 1 doesn't finding a solution means the two vectors of S are
# not independent.

solve(identity(a[1]*SC[1]+a[2]*SC[2]=f, x), [a[1], a[2]])[];
try
  solve(identity(a[1]*SC[1]+a[2]*SC[2]=f, x), [a[1], a[2]])[];
  check = eval(a[1]*SC[1]+a[2]*SC[2] - f, %);
catch:
  printf("The vectors in set S are not independent")
end try;

SU := {Vector(2, {(1) = x, (2) = 4}), Vector(2, {(1) = 2*x, (2) = 8})}

 

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

 

The vectors in set S are not independent

 

# And probably a lot of other methods can be used
 

Download A_few_methods.mw

4 5 6 7 8 9 10 Last Page 6 of 53