6279 Reputation

7 years, 325 days

From the expresison of the Peng-Rob...

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))
 (1)
 > relation_6 := A = a*P/R^2/T^2; relation_7 := B = b*P/R/T; relation_8 := Z = P*v/(R*T);
 (2)
 > eos_1 := eval(eos, isolate(relation_6, a));
 (3)
 > eos_2 := eval(eos_1, isolate(relation_7, b));
 (4)
 > eos_3 := eval(eos_2, isolate(relation_8, P));
 (5)
 > eos_4 := numer(lhs(eos_3))*denom(rhs(eos_3)) - numer(rhs(eos_3))*denom(lhs(eos_3)) = 0
 (6)
 > eos_5 := collect(eos_4, Z)
 (7)
 > eos_6 := factor(eos_5)
 (8)
 > # assuming Z^2*R*T*v^2 <> 0: eos_7 := collect(op(-1, lhs(eos_6)), Z) = 0
 (9)
 >

One condition is missing...

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).

 >
 >
 >
 (1)
 >
 (2)
 >
 (3)
 >
 (4)
 >
 >
 (5)
 > dsolve(EQ1); sol1 := unapply(rhs(%), x); sol2 := dsolve(EQ2); sol2 := unapply(eval(rhs(sol2), [_C1=_C3, _C2=_C4]), x);
 (6)
 > BC1 := select(has, {bc}, U[1, n]); BC2 := select(has, {bc}, U[2, n]);
 (7)
 > bc1 := remove(has, BC1, U[2, n]); bc2 := remove(has, BC2, U[1, n]);
 (8)
 > condition1 := eval(bc1, U[1, n](-Pi)=sol1(-Pi)); condition2 := eval(bc2, U[2, n](-Pi)=sol2(-Pi));
 (9)
 > BC := (BC1 minus bc1) union (BC2 minus bc2); condition34 :=  eval(BC, {U[1, n]=sol1, U[2, n]=sol2});
 (10)
 > Csol := solve({condition1[], condition34[]}, {_C1, _C2, _C3})
 (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
 (12)
 > SolC4free := unapply(eval([sol1(x), sol2(x)], Csol), x): print~(['U'[1, n](x), 'U'[2, n](x)] =~ %(x)):
 (13)
 > # You need another relation to determine _C4

Several errors...

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
 (1)
 > a1 := [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]:
 (3)
 > b1 :=eval(eq1, a1)
 (4)
 > c1 := dsolve({b1, F(-1) = -1, F(1) = 1, (D(F))(-1) = 0, (D(F))(1) = 0}, numeric); c1(0)
 (5)
 >

At least two ways...

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

One Way...

Here is a solution:

One_way.mw (UPDATED)

Formal and numerical resolutions are proposed, the latter approach is done in two different ways.
sol_1.mw

ImageTools ?...

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

The four images:

Statistics:-Trim seems to do what you wa...

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)) ):
 memory used=9.22MiB, alloc change=0 bytes, cpu time=64.00ms, real time=64.00ms, gc time=0ns
 >

varphi_c_c and var_phi_f not define......

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));
 (1)
 > limit(eval(ToInt, data), varphi=0, right); plot(eval(ToInt, data), eval(FromTo, data))
 > 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)

What a shame you changed the left bc !...

... 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 Newton...

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

 >
 >
 >
 (1)
 >
 >
 (2)
 >
 >
 (3)
 >

Use ListTools:-Search...

```ListTools:-Search(96, R);
```

Without even using maplemint one can see...

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

A few solutions...

It is likely that a lot of people will provide you different methods to tackle your problem.

• 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;
 (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);
 (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(%));
 (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;
 The vectors in set S are not independent
 > # And probably a lot of other methods can be used