Advanced Search

Gabriel samaila

35 Reputation

4 Badges

6 years, 124 days
Contact Gabriel samaila

MaplePrimes Activity


These are replies submitted by Gabriel samaila

@vv  very nice is running but iam ...

Gabriel samaila 35
October 09 2018
0 0

@vv 

very nice is running but iam not getting the values on table 1.

@tomleslie  thank you so much...

Gabriel samaila 35
October 09 2018
0 0

@tomleslie 

thank you so much seen and corrected

No=0...

Gabriel samaila 35
June 17 2018
0 0

@Carl Love 

ok, thank you so much for your efforts.

I’ll try it...

Gabriel samaila 35
June 17 2018
0 0

@Carl Love 

because these parameters were not in ref33 that is the results that the author compare his results with. I believe he has either make them zero or played with them when calculating the skin friction but Nb cannot be zero because is denominator..

Trial...

Gabriel samaila 35
June 17 2018
0 0

@Carl Love 

hi, have you tried using Nt, Br, and Sc to be zero while Nt a small value in the fsolve code you developed? In calculating the skin friction table?

 

Nt and Nb...

Gabriel samaila 35
May 29 2018
0 0

@Carl Love 

The author extend reference 33 and in that paper there was no Nb any Sc and Br. For the author to compare his skin friction with reference 33 he must have assume almost all these value to be zero except Nb which is the denominator be it be reduced to a small value. That is why the author didn’t put the value of those values in skin friction table 

From 3D plots...

Gabriel samaila 35
May 28 2018
0 0

@Carl Love 

i have sent them an email but yet the didn’t respond 

Range...

Gabriel samaila 35
May 28 2018
0 0

@Carl Love 

i think -1 to +1

Fsolve...

Gabriel samaila 35
May 27 2018
0 0

@Carl Love 

the default values are Nt and Nb=0.1, and Br and Sc =1 as in nussel number table 

table 3...

Gabriel samaila 35
May 25 2018
0 0

@Carl Love 

please dont get tied of me, i want to be guessing those four values untill i get them. i have tried to edit your code (sols1) but is not editing, how do i go about it?

 

 

 

nanofluid_BVP_(2).mw

Pass or fail...

Gabriel samaila 35
May 22 2018
0 0

@Carl Love 

i email the author on the Newton part and hasn’t responded I also solve the problem using partubation, every thing came up but the skin friction hasn’t and that is what my professor want to see.

I appreciate...

Gabriel samaila 35
May 22 2018
0 0

@Carl Love

you have solve my problem completely. The only thing I tried is not coming up is the skin friction and the nussel number.

Table 1, table2 and 3

no need to convert...

Gabriel samaila 35
May 21 2018
0 0

@Carl Love 

below is the code i develop for my problem, the code is fine but is does not solve skin friction correctly. 

                              NULL

> restart;
> eqn1 := {diff(B(y), y, y)+V[0]*P[m]*(diff(B(y), y))+H[a]*P[m]*(diff(U(y), y)) = 0, diff(phi(y), y, y)+N[t]*(diff(theta(y), y, y))/N[b]+V[0]*S[c]*(diff(phi(y), y)) = 0, diff(theta(y), y, y)+N[t]*(diff(theta(y), y))^2+V[0]*P[r]*(diff(theta(y), y))+N[b]*(diff(theta(y), y))*(diff(phi(y), y)) = 0, diff(U(y), y, y)+H[a]*(diff(B(y), y))+theta(y)+B[r]*phi(y)+V[0]*(diff(U(y), y)) = 0, B(0) = 0, U(0) = 0, U(1) = 0, phi(0) = 1, phi(1) = 0, theta(1) = 0, (D(B))(1) = 0, (D(theta))(0) = -1};
 /                                                                 
 |                                                                 
< / d  / d      \\             / d      \             / d      \   
 ||--- |--- B(y)|| + V[0] P[m] |--- B(y)| + H[a] P[m] |--- U(y)| = 
 \\ dy \ dy     //             \ dy     /             \ dy     /   

                               / d  / d          \\
                          N[t] |--- |--- theta(y)||
     / d  / d        \\        \ dy \ dy         //
  0, |--- |--- phi(y)|| + -------------------------
     \ dy \ dy       //             N[b]           

               / d        \      / d  / d          \\
   + V[0] S[c] |--- phi(y)| = 0, |--- |--- theta(y)||
               \ dy       /      \ dy \ dy         //

                        2                           
          / d          \              / d          \
   + N[t] |--- theta(y)|  + V[0] P[r] |--- theta(y)|
          \ dy         /              \ dy         /

          / d          \ / d        \      / d  / d      \\
   + N[b] |--- theta(y)| |--- phi(y)| = 0, |--- |--- U(y)||
          \ dy         / \ dy       /      \ dy \ dy     //

          / d      \                                 / d      \   
   + H[a] |--- B(y)| + theta(y) + B[r] phi(y) + V[0] |--- U(y)| = 
          \ dy     /                                 \ dy     /   

  0, B(0) = 0, U(0) = 0, U(1) = 0, phi(0) = 1, phi(1) = 0, 

                                             \ 
                                             | 
                                              >
  theta(1) = 0, D(B)(1) = 0, D(theta)(0) = -1| 
                                             / 
> sys1 := eval(eqn1, {B[r] = 1, H[a] = 5, N[b] = .1, N[t] = .1, P[m] = .5, P[r] = .7, S[c] = 1, V[0] = .5});
 /                                                               
 |/ d  / d      \\        / d      \       / d      \      / d  /
< |--- |--- B(y)|| + 0.25 |--- B(y)| + 2.5 |--- U(y)| = 0, |--- |
 |\ dy \ dy     //        \ dy     /       \ dy     /      \ dy \
 \                                                               

   d        \\               / d  / d          \\
  --- phi(y)|| + 1.000000000 |--- |--- theta(y)||
   dy       //               \ dy \ dy         //

         / d        \      / d  / d          \\
   + 0.5 |--- phi(y)| = 0, |--- |--- theta(y)||
         \ dy       /      \ dy \ dy         //

                       2                      
         / d          \         / d          \
   + 0.1 |--- theta(y)|  + 0.35 |--- theta(y)|
         \ dy         /         \ dy         /

         / d          \ / d        \      / d  / d      \\
   + 0.1 |--- theta(y)| |--- phi(y)| = 0, |--- |--- U(y)||
         \ dy         / \ dy       /      \ dy \ dy     //

       / d      \                           / d      \      
   + 5 |--- B(y)| + theta(y) + phi(y) + 0.5 |--- U(y)| = 0, 
       \ dy     /                           \ dy     /      

  B(0) = 0, U(0) = 0, U(1) = 0, phi(0) = 1, phi(1) = 0, 

                                             \ 
                                             | 
  theta(1) = 0, D(B)(1) = 0, D(theta)(0) = -1 >
                                             | 
                                             / 
> 
> sol1 := dsolve(sys1, numeric, output = array([0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1.0]));

> 
> 
> 

> with(plots);
> sol1 := dsolve(sys1, type = numeric, range = 0 .. 1);
> 
> odeplot(sol1, [y, diff(U(y), y)], 0 .. 1, numpoints = 100);

> odeplot(sol1, [y, B(y)], 0 .. 1, numpoints = 25);

> odeplot(sol1, [y, theta(y)], 0 .. 1, numpoints = 25);


> 
> 
> odeplot(sol1, [y, phi(y)], 0 .. 1, numpoints = 25);

> 
> 

No need to convert...

Gabriel samaila 35
May 21 2018
0 0

@Carl Love 

i did that using the default numerical method but I could get the skin friction correct, I tried to used rk4 instead but error said you can solve second order equation using Runge Kutta 4 method . That is why I go by the first order.

thanks very helpful, i have fixed that....

Gabriel samaila 35
May 21 2018
0 0

@Carl Love 

I appreciate your effort. i did not know that maple is flexible like that.

But my answer is not comparable with the paper im reviewing.

based on my calculation, using newton rapson method i got u1=0,u2=0, u3=1, and u4=0.64. though i may be wrong.

can any one look through my code and paper(dimensionless)? check the dimensionless equation , boundary conditions and the transformation of second order to first order.

sunday.mwsunday.mwsunday.mw
 

NULL

> 

restart
> 

dsys := {diff(x[1](y), y) = 1, diff(x[2](y), y) = x[3](y), diff(x[3](y), y) = -x[6](y)-H[a]*x[5](y)-V[0]*x[3](y)-B[r]*x[8](y), diff(x[4](y), y) = x[5](y), diff(x[5](y), y) = -H[a]*P[m]*x[3](y)-V[0]*P[m]*x[5](y), diff(x[6](y), y) = x[7](y), diff(x[7](y), y) = -N[t]*x[7](y)^2-V[0]*P[r]*x[7](y)-N[b]*x[7](y)*x[9](y), diff(x[8](y), y) = x[9](y), diff(x[9](y), y) = -V[0]*S[c]*x[9](y)+N[t]^2*x[7](y)^2/N[b]+N[t]*P[r]*V[0]*x[7](y)/N[b]+N[t]*x[7](y)*x[9](y), x[1](0) = 0, x[2](0) = 0, x[3](0) = 1, x[4](0) = 0, x[5](0) = 1, x[6](0) = .1, x[7](0) = -1, x[8](0) = 1, x[9](0) = 0};

{diff(x[1](y), y) = 1, diff(x[2](y), y) = x[3](y), diff(x[3](y), y) = -x[6](y)-H[a]*x[5](y)-V[0]*x[3](y)-B[r]*x[8](y), diff(x[4](y), y) = x[5](y), diff(x[5](y), y) = -H[a]*P[m]*x[3](y)-V[0]*P[m]*x[5](y), diff(x[6](y), y) = x[7](y), diff(x[7](y), y) = -N[t]*x[7](y)^2-V[0]*P[r]*x[7](y)-N[b]*x[7](y)*x[9](y), diff(x[8](y), y) = x[9](y), diff(x[9](y), y) = -V[0]*S[c]*x[9](y)+N[t]^2*x[7](y)^2/N[b]+N[t]*P[r]*V[0]*x[7](y)/N[b]+N[t]*x[7](y)*x[9](y), x[1](0) = 0, x[2](0) = 0, x[3](0) = 1, x[4](0) = 0, x[5](0) = 1, x[6](0) = .1, x[7](0) = -1, x[8](0) = 1, x[9](0) = 0}

 (1)
> 

``
> 

sys1 := eval(dsys, {B[r] = 1, H[a] = 5, N[b] = .1, N[t] = .1, P[m] = .8, P[r] = 10, S[c] = 1, V[0] = .2});

{diff(x[1](y), y) = 1, diff(x[2](y), y) = x[3](y), diff(x[3](y), y) = -x[6](y)-5*x[5](y)-.2*x[3](y)-x[8](y), diff(x[4](y), y) = x[5](y), diff(x[5](y), y) = -4.0*x[3](y)-.16*x[5](y), diff(x[6](y), y) = x[7](y), diff(x[7](y), y) = -.1*x[7](y)^2-2.0*x[7](y)-.1*x[7](y)*x[9](y), diff(x[8](y), y) = x[9](y), diff(x[9](y), y) = -.2*x[9](y)+.1000000000*x[7](y)^2+2.000000000*x[7](y)+.1*x[7](y)*x[9](y), x[1](0) = 0, x[2](0) = 0, x[3](0) = 1, x[4](0) = 0, x[5](0) = 1, x[6](0) = .1, x[7](0) = -1, x[8](0) = 1, x[9](0) = 0}

 (2)
> 

sol := dsolve(sys1, numeric, method = classical[rk4], stepsize = 0.1e-1);

proc (x_classical) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_classical) else _xout := evalf(x_classical) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _n, _y0, _yn, _yl, _ctl, _octl, _reinit, _errcd, _fcn, _i, _yini, _pars, _ini, _par; option `Copyright (c) 2002 by the University of Waterloo. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _ctl := Array(1..9, {(1) = 9, (2) = 0., (3) = 0., (4) = 0.1e-1, (5) = 50000, (6) = 0, (7) = 4, (8) = 1, (9) = 0}); _octl := Array(1..9, {(1) = 9, (2) = 0., (3) = 0., (4) = 0.1e-1, (5) = 50000, (6) = 0, (7) = 4, (8) = 1, (9) = 0}); _n := trunc(_ctl[1]); _yini := Array(0..9, {(1) = 0., (2) = 0., (3) = 0., (4) = 1., (5) = 0., (6) = 1., (7) = .1, (8) = -1., (9) = 1.}); _y0 := Array(0..9, {(1) = 0., (2) = 0., (3) = 0., (4) = 1., (5) = 0., (6) = 1., (7) = .1, (8) = -1., (9) = 1.}); _yl := Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}); _yn := Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}); _fcn := proc (N, X, Y, YP) option `[Y[1] = x[1](y), Y[2] = x[2](y), Y[3] = x[3](y), Y[4] = x[4](y), Y[5] = x[5](y), Y[6] = x[6](y), Y[7] = x[7](y), Y[8] = x[8](y), Y[9] = x[9](y)]`; YP[1] := 1; YP[3] := -Y[6]-5*Y[5]-.2*Y[3]-Y[8]; YP[5] := -4.0*Y[3]-.16*Y[5]; YP[7] := -.1*Y[7]^2-2.0*Y[7]-.1*Y[7]*Y[9]; YP[9] := -.2*Y[9]+.1000000000*Y[7]^2+2.000000000*Y[7]+.1*Y[7]*Y[9]; YP[2] := Y[3]; YP[4] := Y[5]; YP[6] := Y[7]; YP[8] := Y[9]; 0 end proc; _pars := []; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then return _y0[0] elif _xout = "method" then return "classical" elif _xout = "numfun" then return round(_ctl[6]) elif _xout = "initial" then return [seq(_yini[_i], _i = 0 .. _n)] elif _xout = "parameters" then return [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_yini[_i], _i = 0 .. _n)], [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _ctl[3]-_y0[0] <> 0. then _xout := _ctl[3] elif _ctl[2]-_y0[0] <> 0. then _xout := _ctl[2] else error "no information is available on last computed point" end if elif _xout = "enginedata" then return eval(_octl, 1) elif _xout = "function" then return eval(_fcn, 1) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _yini) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n, _ini, _yini, _pars) end if; if _pars <> [] then _par := {seq(rhs(_pars[_i]) = _yini[_n+_i], _i = 1 .. nops(_pars))}; for _i from 0 to _n do _y0[_i] := subs(_par, _yini[_i]) end do; for _i from _n+1 to _n+nops(_pars) do _y0[_i] := _yini[_i] end do; for _i to nops(_pars) do _yn[_n+_i] := _y0[_n+_i] end do else for _i from 0 to _n do _y0[_i] := _yini[_i] end do end if; _octl[2] := _y0[0]; _octl[3] := _y0[0]; for _i to 9 do _ctl[_i] := _octl[_i] end do; for _i to _n+nops(_pars) do _yl[_i] := _y0[_i] end do; if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [seq(_yini[_i], _i = 0 .. _n)] elif _xout = "parameters" then return [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] else return [seq(_yini[_i], _i = 0 .. _n), op(_pars)], [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] end if else return "procname" end if end if; if _xout-_y0[0] = 0. then return [seq(_y0[_i], _i = 0 .. _n)] end if; _reinit := false; if _xin <> "last" then if 0 < 0 and `dsolve/numeric/checkglobals`(0, table( [ ] ), _pars, _n, _yini) then _reinit := true; if _pars <> [] then _par := {seq(rhs(_pars[_i]) = _yini[_n+_i], _i = 1 .. nops(_pars))}; for _i from 0 to _n do _y0[_i] := subs(_par, _yini[_i]) end do; for _i from _n+1 to _n+nops(_pars) do _y0[_i] := _yini[_i] end do; for _i to nops(_pars) do _yn[_n+_i] := _y0[_n+_i] end do else for _i from 0 to _n do _y0[_i] := _yini[_i] end do end if end if; if _pars <> [] and select(type, {seq(_yini[_n+_i], _i = 1 .. nops(_pars))}, 'undefined') <> {} then error "parameters must be initialized before solution can be computed" end if end if; if not _reinit and _xout-_ctl[3] = 0. then [_ctl[3], seq(_yn[_i], _i = 1 .. _n)] elif not _reinit and _xout-_ctl[2] = 0. then [_ctl[2], seq(_yl[_i], _i = 1 .. _n)] else if _ctl[2]-_y0[0] = 0. or sign(_ctl[2]-_y0[0]) <> sign(_xout-_ctl[2]) or abs(_xout-_y0[0]) < abs(_xout-_ctl[2]) or _reinit then for _i to 9 do _ctl[_i] := _octl[_i] end do; for _i to _n+nops(_pars) do _yl[_i] := _y0[_i] end do; for _i to nops(_pars) do _yn[_n+_i] := _y0[_n+_i] end do end if; _ctl[3] := _xout; if Digits <= evalhf(Digits) then try _errcd := evalhf(`dsolve/numeric/classical`(_fcn, var(_ctl), _y0[0], var(_yl), var(_yn), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), var(Array(1..3, 1..9, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 0})))) catch: userinfo(2, `dsolve/debug`, print(`Exception in classical:`, [lastexception])); if searchtext('evalhf', lastexception[2]) <> 0 or searchtext('real', lastexception[2]) <> 0 or searchtext('hardware', lastexception[2]) <> 0 then _errcd := `dsolve/numeric/classical`(_fcn, _ctl, _y0[0], _yl, _yn, Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..3, 1..9, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 0})) else _ctl[3] := _y0[0]; _errcd := [lastexception][2 .. -1] end if end try else try _errcd := `dsolve/numeric/classical`(_fcn, _ctl, _y0[0], _yl, _yn, Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..3, 1..9, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 0})) catch: _ctl[3] := _y0[0]; _errcd := [lastexception][2 .. -1] end try end if; if type(_errcd, 'list') or _errcd < 0 then userinfo(2, {dsolve, `dsolve/classical`}, `Last values returned:`); userinfo(2, {dsolve, `dsolve/classical`}, ` t =`, _ctl[2]); userinfo(2, {dsolve, `dsolve/classical`}, ` y =`, _yl[1]); for _i from 2 to _n do userinfo(2, {dsolve, `dsolve/classical`}, ` `, _yl[_i]) end do; if type(_errcd, 'list') then error op(_errcd) elif _errcd+1. = 0. then Rounding := `if`(_y0[0] < _xout, -infinity, infinity); error "cannot evaluate the solution past %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_ctl[3]) else Rounding := `if`(_y0[0] < _xout, -infinity, infinity); error "cannot evaluate the solution past %1, unknown error code returned from classical %2", evalf[8](_ctl[3]), trunc(_errcd) end if end if; if _Env_smart_dsolve_numeric = true then if _y0[0] < _xout and procname("right") < _xout then procname("right") := _xout elif _xout < _y0[0] and _xout < procname("left") then procname("left") := _xout end if end if; [_xout, seq(_yn[_i], _i = 1 .. _n)] end if end proc, (2) = Array(0..0, {}), (3) = [y, x[1](y), x[2](y), x[3](y), x[4](y), x[5](y), x[6](y), x[7](y), x[8](y), x[9](y)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_classical, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_classical, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_classical, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_classical, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_classical), 'string') = rhs(x_classical); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_classical), 'string') = rhs(x_classical)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_classical) else _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_classical) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc

 (3)
> 

plots:-odeplot(sol, [[y, x[6](y)]], y = 0 .. 1, color = ["Red"])

 
> 

``
> 

``

``


 

Download sunday.mw
 

NULL

> 

restart
> 

dsys := {diff(x[1](y), y) = 1, diff(x[2](y), y) = x[3](y), diff(x[3](y), y) = -x[6](y)-H[a]*x[5](y)-V[0]*x[3](y)-B[r]*x[8](y), diff(x[4](y), y) = x[5](y), diff(x[5](y), y) = -H[a]*P[m]*x[3](y)-V[0]*P[m]*x[5](y), diff(x[6](y), y) = x[7](y), diff(x[7](y), y) = -N[t]*x[7](y)^2-V[0]*P[r]*x[7](y)-N[b]*x[7](y)*x[9](y), diff(x[8](y), y) = x[9](y), diff(x[9](y), y) = -V[0]*S[c]*x[9](y)+N[t]^2*x[7](y)^2/N[b]+N[t]*P[r]*V[0]*x[7](y)/N[b]+N[t]*x[7](y)*x[9](y), x[1](0) = 0, x[2](0) = 0, x[3](0) = 1, x[4](0) = 0, x[5](0) = 1, x[6](0) = .1, x[7](0) = -1, x[8](0) = 1, x[9](0) = 0};

{diff(x[1](y), y) = 1, diff(x[2](y), y) = x[3](y), diff(x[3](y), y) = -x[6](y)-H[a]*x[5](y)-V[0]*x[3](y)-B[r]*x[8](y), diff(x[4](y), y) = x[5](y), diff(x[5](y), y) = -H[a]*P[m]*x[3](y)-V[0]*P[m]*x[5](y), diff(x[6](y), y) = x[7](y), diff(x[7](y), y) = -N[t]*x[7](y)^2-V[0]*P[r]*x[7](y)-N[b]*x[7](y)*x[9](y), diff(x[8](y), y) = x[9](y), diff(x[9](y), y) = -V[0]*S[c]*x[9](y)+N[t]^2*x[7](y)^2/N[b]+N[t]*P[r]*V[0]*x[7](y)/N[b]+N[t]*x[7](y)*x[9](y), x[1](0) = 0, x[2](0) = 0, x[3](0) = 1, x[4](0) = 0, x[5](0) = 1, x[6](0) = .1, x[7](0) = -1, x[8](0) = 1, x[9](0) = 0}

 (1)
> 

``
> 

sys1 := eval(dsys, {B[r] = 1, H[a] = 5, N[b] = .1, N[t] = .1, P[m] = .8, P[r] = 10, S[c] = 1, V[0] = .2});

{diff(x[1](y), y) = 1, diff(x[2](y), y) = x[3](y), diff(x[3](y), y) = -x[6](y)-5*x[5](y)-.2*x[3](y)-x[8](y), diff(x[4](y), y) = x[5](y), diff(x[5](y), y) = -4.0*x[3](y)-.16*x[5](y), diff(x[6](y), y) = x[7](y), diff(x[7](y), y) = -.1*x[7](y)^2-2.0*x[7](y)-.1*x[7](y)*x[9](y), diff(x[8](y), y) = x[9](y), diff(x[9](y), y) = -.2*x[9](y)+.1000000000*x[7](y)^2+2.000000000*x[7](y)+.1*x[7](y)*x[9](y), x[1](0) = 0, x[2](0) = 0, x[3](0) = 1, x[4](0) = 0, x[5](0) = 1, x[6](0) = .1, x[7](0) = -1, x[8](0) = 1, x[9](0) = 0}

 (2)
> 

sol := dsolve(sys1, numeric, method = classical[rk4], stepsize = 0.1e-1);

proc (x_classical) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_classical) else _xout := evalf(x_classical) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _n, _y0, _yn, _yl, _ctl, _octl, _reinit, _errcd, _fcn, _i, _yini, _pars, _ini, _par; option `Copyright (c) 2002 by the University of Waterloo. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _ctl := Array(1..9, {(1) = 9, (2) = 0., (3) = 0., (4) = 0.1e-1, (5) = 50000, (6) = 0, (7) = 4, (8) = 1, (9) = 0}); _octl := Array(1..9, {(1) = 9, (2) = 0., (3) = 0., (4) = 0.1e-1, (5) = 50000, (6) = 0, (7) = 4, (8) = 1, (9) = 0}); _n := trunc(_ctl[1]); _yini := Array(0..9, {(1) = 0., (2) = 0., (3) = 0., (4) = 1., (5) = 0., (6) = 1., (7) = .1, (8) = -1., (9) = 1.}); _y0 := Array(0..9, {(1) = 0., (2) = 0., (3) = 0., (4) = 1., (5) = 0., (6) = 1., (7) = .1, (8) = -1., (9) = 1.}); _yl := Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}); _yn := Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}); _fcn := proc (N, X, Y, YP) option `[Y[1] = x[1](y), Y[2] = x[2](y), Y[3] = x[3](y), Y[4] = x[4](y), Y[5] = x[5](y), Y[6] = x[6](y), Y[7] = x[7](y), Y[8] = x[8](y), Y[9] = x[9](y)]`; YP[1] := 1; YP[3] := -Y[6]-5*Y[5]-.2*Y[3]-Y[8]; YP[5] := -4.0*Y[3]-.16*Y[5]; YP[7] := -.1*Y[7]^2-2.0*Y[7]-.1*Y[7]*Y[9]; YP[9] := -.2*Y[9]+.1000000000*Y[7]^2+2.000000000*Y[7]+.1*Y[7]*Y[9]; YP[2] := Y[3]; YP[4] := Y[5]; YP[6] := Y[7]; YP[8] := Y[9]; 0 end proc; _pars := []; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then return _y0[0] elif _xout = "method" then return "classical" elif _xout = "numfun" then return round(_ctl[6]) elif _xout = "initial" then return [seq(_yini[_i], _i = 0 .. _n)] elif _xout = "parameters" then return [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_yini[_i], _i = 0 .. _n)], [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _ctl[3]-_y0[0] <> 0. then _xout := _ctl[3] elif _ctl[2]-_y0[0] <> 0. then _xout := _ctl[2] else error "no information is available on last computed point" end if elif _xout = "enginedata" then return eval(_octl, 1) elif _xout = "function" then return eval(_fcn, 1) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _yini) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n, _ini, _yini, _pars) end if; if _pars <> [] then _par := {seq(rhs(_pars[_i]) = _yini[_n+_i], _i = 1 .. nops(_pars))}; for _i from 0 to _n do _y0[_i] := subs(_par, _yini[_i]) end do; for _i from _n+1 to _n+nops(_pars) do _y0[_i] := _yini[_i] end do; for _i to nops(_pars) do _yn[_n+_i] := _y0[_n+_i] end do else for _i from 0 to _n do _y0[_i] := _yini[_i] end do end if; _octl[2] := _y0[0]; _octl[3] := _y0[0]; for _i to 9 do _ctl[_i] := _octl[_i] end do; for _i to _n+nops(_pars) do _yl[_i] := _y0[_i] end do; if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [seq(_yini[_i], _i = 0 .. _n)] elif _xout = "parameters" then return [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] else return [seq(_yini[_i], _i = 0 .. _n), op(_pars)], [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] end if else return "procname" end if end if; if _xout-_y0[0] = 0. then return [seq(_y0[_i], _i = 0 .. _n)] end if; _reinit := false; if _xin <> "last" then if 0 < 0 and `dsolve/numeric/checkglobals`(0, table( [ ] ), _pars, _n, _yini) then _reinit := true; if _pars <> [] then _par := {seq(rhs(_pars[_i]) = _yini[_n+_i], _i = 1 .. nops(_pars))}; for _i from 0 to _n do _y0[_i] := subs(_par, _yini[_i]) end do; for _i from _n+1 to _n+nops(_pars) do _y0[_i] := _yini[_i] end do; for _i to nops(_pars) do _yn[_n+_i] := _y0[_n+_i] end do else for _i from 0 to _n do _y0[_i] := _yini[_i] end do end if end if; if _pars <> [] and select(type, {seq(_yini[_n+_i], _i = 1 .. nops(_pars))}, 'undefined') <> {} then error "parameters must be initialized before solution can be computed" end if end if; if not _reinit and _xout-_ctl[3] = 0. then [_ctl[3], seq(_yn[_i], _i = 1 .. _n)] elif not _reinit and _xout-_ctl[2] = 0. then [_ctl[2], seq(_yl[_i], _i = 1 .. _n)] else if _ctl[2]-_y0[0] = 0. or sign(_ctl[2]-_y0[0]) <> sign(_xout-_ctl[2]) or abs(_xout-_y0[0]) < abs(_xout-_ctl[2]) or _reinit then for _i to 9 do _ctl[_i] := _octl[_i] end do; for _i to _n+nops(_pars) do _yl[_i] := _y0[_i] end do; for _i to nops(_pars) do _yn[_n+_i] := _y0[_n+_i] end do end if; _ctl[3] := _xout; if Digits <= evalhf(Digits) then try _errcd := evalhf(`dsolve/numeric/classical`(_fcn, var(_ctl), _y0[0], var(_yl), var(_yn), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), var(Array(1..3, 1..9, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 0})))) catch: userinfo(2, `dsolve/debug`, print(`Exception in classical:`, [lastexception])); if searchtext('evalhf', lastexception[2]) <> 0 or searchtext('real', lastexception[2]) <> 0 or searchtext('hardware', lastexception[2]) <> 0 then _errcd := `dsolve/numeric/classical`(_fcn, _ctl, _y0[0], _yl, _yn, Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..3, 1..9, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 0})) else _ctl[3] := _y0[0]; _errcd := [lastexception][2 .. -1] end if end try else try _errcd := `dsolve/numeric/classical`(_fcn, _ctl, _y0[0], _yl, _yn, Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..3, 1..9, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 0})) catch: _ctl[3] := _y0[0]; _errcd := [lastexception][2 .. -1] end try end if; if type(_errcd, 'list') or _errcd < 0 then userinfo(2, {dsolve, `dsolve/classical`}, `Last values returned:`); userinfo(2, {dsolve, `dsolve/classical`}, ` t =`, _ctl[2]); userinfo(2, {dsolve, `dsolve/classical`}, ` y =`, _yl[1]); for _i from 2 to _n do userinfo(2, {dsolve, `dsolve/classical`}, ` `, _yl[_i]) end do; if type(_errcd, 'list') then error op(_errcd) elif _errcd+1. = 0. then Rounding := `if`(_y0[0] < _xout, -infinity, infinity); error "cannot evaluate the solution past %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_ctl[3]) else Rounding := `if`(_y0[0] < _xout, -infinity, infinity); error "cannot evaluate the solution past %1, unknown error code returned from classical %2", evalf[8](_ctl[3]), trunc(_errcd) end if end if; if _Env_smart_dsolve_numeric = true then if _y0[0] < _xout and procname("right") < _xout then procname("right") := _xout elif _xout < _y0[0] and _xout < procname("left") then procname("left") := _xout end if end if; [_xout, seq(_yn[_i], _i = 1 .. _n)] end if end proc, (2) = Array(0..0, {}), (3) = [y, x[1](y), x[2](y), x[3](y), x[4](y), x[5](y), x[6](y), x[7](y), x[8](y), x[9](y)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_classical, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_classical, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_classical, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_classical, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_classical), 'string') = rhs(x_classical); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_classical), 'string') = rhs(x_classical)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_classical) else _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_classical) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc

 (3)
> 

plots:-odeplot(sol, [[y, x[6](y)]], y = 0 .. 1, color = ["Red"])

 
> 

``
> 

``

``


 

Download sunday.mw

 

 

DIMENSIONLESS.pdfDIMENSIONLESS.pdfDIMENSIONLESS.pdfDIMENSIONLESS.pdfDIMENSIONLESS.pdfDIMENSIONLESS.pdf

1 2 3 4 Page 3 of 4