mmcdara

3437 Reputation

16 Badges

5 years, 230 days

MaplePrimes Activity


These are answers submitted by mmcdara



I used Maple 2015.2

Main problem; the value of delta is not set, I arbitrarily ghoosed delta:=1:

I changed theta[p] into theta__p (not sure this has any influence but I fell safer to use theta__p which is a variable which has noth1ng to do with theta).
The main point is that your equations 1 and 2 do not contain the functtions theta(eta) nor theta__p(eta).
This mean the inititial system can be spilt into 2 sub-systems:

  • sys__12 made of eq1, eq2 and their bcs,
  • sys__34 made of eq3, eq4 and their bcs.

Thus:

  1. dsolve sys__12 (I use the option output=procedurelist)
     
  2. Define a procedure K wich takes as argument the name eta or a numerical value of eta and returns in this case the evaluations of [f(eta), F(eta), diff(f(eta), eta), diff(F(eta), eta)] for this value of theta.
     
  3. After having instanciate sys__34 with params, replace:
    1. f(eta) by K(eta)[1],
    2. F(eta) by K(eta)[2],
    3. diff(f(eta), eta) by K(eta)[3],
    4. diff(F(eta), eta) by K(eta)[4],
       
  4. Solve the resulting sys__34 system.
    You must use the option known=K to tell Maple that K(eta)[1]...K(eta)[4] are known elsewhere.know

    Without any precaution you should get ab error which advices you tou use the midpoint submethod either midrich or middefer.
    Whatever the submethod I used the system seems so badly conditionned that I wasn't able to reach a reasonable accuracy (I desperately set abserr=1 to obtain a solution). This what I titled my answer "almost done".

restart

with(plots):

eq1 := (2*eta*gamma+1)*(diff(f(eta), `$`(eta, 3)))+2*gamma*(diff(f(eta), `$`(eta, 2)))+f(eta)*(diff(f(eta), `$`(eta, 2)))-(diff(f(eta), eta))^2-(Q+S)*(diff(f(eta), eta))+beta*(diff(F(eta), eta)-(diff(f(eta), eta))) = 0:

eq2 := (diff(F(eta), `$`(eta, 2)))*F(eta)-(diff(F(eta), eta))^2+beta*(diff(f(eta), eta)-(diff(F(eta), eta))) = 0:

eq3 := (2*eta*gamma+1)*(1+Rd)*(diff(theta(eta), `$`(eta, 2))) = -Pr*((diff(theta(eta), eta))*f(eta)-2*(diff(f(eta), eta))*theta(eta))-gamma*(diff(theta(eta), eta))-N*Pr*betat*((theta__p(eta), eta)-theta(eta))-N*Pr*Ec*betat*(diff(F(eta), eta)-(diff(f(eta), eta)))-Pr*delta*theta(eta):

eq4 := 2*(diff(theta__p(eta), eta))*f(eta)-F(eta)*theta__p(eta)+betat*delta*(theta__p(eta)-theta(eta)) = 0:

bcs := f(0) = 0, (D(f))(0) = 1, (D(f))(5) = 0, (D(F))(5) = 0, F(5) = f(5), theta(0) = 1, theta(5) = 0, theta__p(5) = 0:

params := [Rd = .1, beta = .5, Q = .5, S = .5, gamma = .1, Pr = 6.2, N = .5, betat = .5, Ec = .1]:

delta := 1:  # delta not fixed !!!

# Main observation; eq1 and eq2 depend only upon f and F  

indets(eq1, function);
indets(eq2, function);

{F(eta), diff(F(eta), eta), diff(diff(diff(f(eta), eta), eta), eta), diff(diff(f(eta), eta), eta), diff(f(eta), eta), f(eta)}

 

{F(eta), diff(F(eta), eta), diff(diff(F(eta), eta), eta), diff(f(eta), eta), f(eta)}

(1)

sys__12 := eval([eq1, eq2, f(0) = 0, (D(f))(0) = 1, (D(f))(5) = 0, (D(F))(5) = 0, F(5) = f(5)], params):
sol__12 := dsolve(sys__12, numeric, output = procedurelist);

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(17, {(1) = .0, (2) = .2700295406680469, (3) = .5461814747766728, (4) = .8292588209972042, (5) = 1.1234150774760683, (6) = 1.431035339629634, (7) = 1.7496291223773444, (8) = 2.075914881269067, (9) = 2.408230078854766, (10) = 2.742205295476464, (11) = 3.0774349455365573, (12) = 3.4129672931189416, (13) = 3.74869641693143, (14) = 4.084468803510717, (15) = 4.417149968189033, (16) = 4.721816642175481, (17) = 5.0}, datatype = float[8], order = C_order); Y := Matrix(17, 5, {(1, 1) = .44776527313438064, (1, 2) = .33133727618764913, (1, 3) = .0, (1, 4) = 1.0, (1, 5) = -1.5719460272279495, (2, 1) = .5213624951662373, (2, 2) = .2217023909765899, (2, 3) = .22045952998207696, (2, 4) = .6583232859380348, (2, 5) = -1.0037671287227639, (3, 1) = .5718234259032432, (3, 2) = .14894670237473073, (3, 3) = .36912084196073325, (3, 4) = .43487740306649103, (3, 5) = -.6430155320436061, (4, 1) = .6066140717509945, (4, 2) = .10029249480024327, (4, 3) = .469902127251469, (4, 4) = .2879720893049211, (4, 5) = -.4130037744612695, (5, 1) = .6309127884296852, (5, 2) = 0.6725940922492528e-1, (5, 3) = .5391442607170165, (5, 4) = .19006598529443805, (5, 5) = -.26445645144099067, (6, 1) = .6478941036296227, (6, 2) = 0.4475730981457941e-1, (6, 3) = .586800586045253, (6, 4) = .12467762652820233, (6, 5) = -.1683906813056968, (7, 1) = .6595684318193582, (7, 2) = 0.29614870473976596e-1, (7, 3) = .6191402504248181, (7, 4) = 0.8157762164726536e-1, (7, 5) = -.10711429298647333, (8, 1) = .6674672167892126, (8, 2) = 0.1951756867756001e-1, (8, 3) = .6408263407473044, (8, 4) = 0.5344481967149114e-1, (8, 5) = -0.6842360971338068e-1, (9, 1) = .6727561100925025, (9, 2) = 0.12782783201441923e-1, (9, 3) = .6553122229469226, (9, 4) = 0.3508354643667022e-1, (9, 5) = -0.44010428510442734e-1, (10, 1) = .6762273833779221, (10, 2) = 0.8307953068430345e-2, (10, 3) = .6648932981595975, (10, 4) = 0.23142480667606512e-1, (10, 5) = -0.28681021784476287e-1, (11, 1) = .6784756350927372, (11, 2) = 0.5303277034954629e-2, (11, 3) = .6712427107983656, (11, 4) = 0.15279463138533988e-1, (11, 5) = -0.18963855581152253e-1, (12, 1) = .679892179218728, (12, 2) = 0.32714663294575444e-2, (12, 3) = .6754314367159898, (12, 4) = 0.10033737564104571e-1, (12, 5) = -0.12763512056209856e-1, (13, 1) = .6807440254039239, (13, 2) = 0.18925878798273792e-2, (13, 3) = .6781638938053991, (13, 4) = 0.6466085693279025e-2, (13, 5) = -0.878177393716817e-2, (14, 1) = .6812136245064051, (14, 2) = 0.9684873470550855e-3, (14, 3) = .6798934576626122, (14, 4) = 0.3978616235233124e-2, (14, 5) = -0.62238500751895555e-2, (15, 1) = .6814306856085619, (15, 2) = 0.3850239255045409e-3, (15, 3) = .6809058823403, (15, 4) = 0.21970566082543854e-2, (15, 5) = -0.4609984823372287e-2, (16, 1) = .6814973367424689, (16, 2) = 0.8800363841388453e-4, (16, 3) = .6813773865384698, (16, 4) = 0.9452622479871761e-3, (16, 5) = -0.36799735133912046e-2, (17, 1) = .6815055279510432, (17, 2) = .0, (17, 3) = .6815055279510432, (17, 4) = .0, (17, 5) = -0.3161424611013035e-2}, datatype = float[8], order = C_order); YP := Matrix(17, 5, {(1, 1) = .33133727618764913, (1, 2) = -.5014836674199502, (1, 3) = 1.0, (1, 4) = -1.5719460272279495, (1, 5) = 2.648720567351765, (2, 1) = .2217023909765899, (2, 2) = -.3244546719112376, (2, 3) = .6583232859380348, (2, 4) = -1.0037671287227639, (2, 5) = 1.643317863794819, (3, 1) = .14894670237473073, (3, 2) = -.21121945119123267, (3, 3) = .43487740306649103, (3, 4) = -.6430155320436061, (3, 5) = 1.0213465383823663, (4, 1) = .10029249480024327, (4, 2) = -.1381128738033511, (4, 3) = .2879720893049211, (4, 4) = -.4130037744612695, (4, 5) = .635940126024083, (5, 1) = 0.6725940922492528e-1, (5, 2) = -0.9015423517890793e-1, (5, 3) = .19006598529443805, (5, 4) = -.26445645144099067, (5, 5) = .3944415120283217, (6, 1) = 0.4475730981457941e-1, (6, 2) = -0.5858510111811698e-1, (6, 3) = .12467762652820233, (6, 4) = -.1683906813056968, (6, 5) = .24309630227876206, (7, 1) = 0.29614870473976596e-1, (7, 2) = -0.3806175951175528e-1, (7, 3) = 0.8157762164726536e-1, (7, 4) = -.10711429298647333, (7, 5) = .1496049121417643, (8, 1) = 0.1951756867756001e-1, (8, 2) = -0.2484420147202396e-1, (8, 3) = 0.5344481967149114e-1, (8, 4) = -0.6842360971338068e-1, (8, 5) = 0.9242420908739084e-1, (9, 1) = 0.12782783201441923e-1, (9, 2) = -0.1633130031287029e-1, (9, 3) = 0.3508354643667022e-1, (9, 4) = -0.44010428510442734e-1, (9, 5) = 0.57441143040287745e-1, (10, 1) = 0.8307953068430345e-2, (10, 2) = -0.10866524923457799e-1, (10, 3) = 0.23142480667606512e-1, (10, 4) = -0.28681021784476287e-1, (10, 5) = 0.3610169213484274e-1, (11, 1) = 0.5303277034954629e-2, (11, 2) = -0.7310458987672048e-2, (11, 3) = 0.15279463138533988e-1, (11, 4) = -0.18963855581152253e-1, (11, 5) = 0.229176337398394e-1, (12, 1) = 0.32714663294575444e-2, (12, 2) = -0.4957305332224504e-2, (12, 3) = 0.10033737564104571e-1, (12, 4) = -0.12763512056209856e-1, (12, 5) = 0.14673258499433704e-1, (13, 1) = 0.18925878798273792e-2, (13, 2) = -0.3353928837624718e-2, (13, 3) = 0.6466085693279025e-2, (13, 4) = -0.878177393716817e-2, (13, 5) = 0.9433680672377034e-2, (14, 1) = 0.9684873470550855e-3, (14, 2) = -0.2208009972550792e-2, (14, 3) = 0.3978616235233124e-2, (14, 4) = -0.62238500751895555e-2, (14, 5) = 0.6040988893516817e-2, (15, 1) = 0.3850239255045409e-3, (15, 2) = -0.13293620570413733e-2, (15, 3) = 0.21970566082543854e-2, (15, 4) = -0.4609984823372287e-2, (15, 5) = 0.38062804457790473e-2, (16, 1) = 0.8800363841388453e-4, (16, 2) = -0.6289409173556955e-3, (16, 3) = 0.9452622479871761e-3, (16, 4) = -0.36799735133912046e-2, (16, 5) = 0.2375189062515268e-2, (17, 1) = .0, (17, 2) = .0, (17, 3) = .0, (17, 4) = -0.3161424611013035e-2, (17, 5) = 0.13934066354042334e-2}, 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(17, {(1) = .0, (2) = .2700295406680469, (3) = .5461814747766728, (4) = .8292588209972042, (5) = 1.1234150774760683, (6) = 1.431035339629634, (7) = 1.7496291223773444, (8) = 2.075914881269067, (9) = 2.408230078854766, (10) = 2.742205295476464, (11) = 3.0774349455365573, (12) = 3.4129672931189416, (13) = 3.74869641693143, (14) = 4.084468803510717, (15) = 4.417149968189033, (16) = 4.721816642175481, (17) = 5.0}, datatype = float[8], order = C_order); Y := Matrix(17, 5, {(1, 1) = -0.5557170630311661e-9, (1, 2) = -0.8568444294619354e-8, (1, 3) = .0, (1, 4) = .0, (1, 5) = -0.996321158246133e-8, (2, 1) = -0.29775367430238932e-7, (2, 2) = 0.51448201034330047e-7, (2, 3) = -0.9338503564926378e-7, (2, 4) = 0.16119701358192726e-6, (2, 5) = -0.2661709248184363e-6, (3, 1) = -0.2891050343668424e-7, (3, 2) = 0.5344899476525872e-7, (3, 3) = -0.9176042528308786e-7, (3, 4) = 0.15799013994899224e-6, (3, 5) = -0.2499917438861638e-6, (4, 1) = -0.1927807403261658e-7, (4, 2) = 0.38893522217972736e-7, (4, 3) = -0.6314408631390969e-7, (4, 4) = 0.11040192879654927e-6, (4, 5) = -0.168450566115583e-6, (5, 1) = -0.9451108871656315e-8, (5, 2) = 0.23506561485624035e-7, (5, 3) = -0.3433544296514957e-7, (5, 4) = 0.6357224004387666e-7, (5, 5) = -0.9346927926922444e-7, (6, 1) = -0.17952890062909867e-8, (6, 2) = 0.11485578153704378e-7, (6, 3) = -0.12371911215185314e-7, (6, 4) = 0.2856188378954896e-7, (6, 5) = -0.40177124898242385e-7, (7, 1) = 0.33629876485582894e-8, (7, 2) = 0.3438144666793098e-8, (7, 3) = 0.19670495465459734e-8, (7, 4) = 0.6170587583276722e-8, (7, 5) = -0.7855938356527589e-8, (8, 1) = 0.6290794653972546e-8, (8, 2) = -0.10736074960720638e-8, (8, 3) = 0.969080021609951e-8, (8, 4) = -0.569460732839434e-8, (8, 5) = 0.8253813531999202e-8, (9, 1) = 0.7572046363127913e-8, (9, 2) = -0.3010207778743067e-8, (9, 3) = 0.12689878914657082e-7, (9, 4) = -0.1033127198991506e-7, (9, 5) = 0.13987204260214792e-7, (10, 1) = 0.7829327930419819e-8, (10, 2) = -0.3347426371048808e-8, (10, 3) = 0.12839070161904223e-7, (10, 4) = -0.1073328829776e-7, (10, 5) = 0.14076279806383137e-7, (11, 1) = 0.7582703863603602e-8, (11, 2) = -0.2891392769667636e-8, (11, 3) = 0.11656345534185066e-7, (11, 4) = -0.9198752266915943e-8, (11, 5) = 0.11835993527707034e-7, (12, 1) = 0.7166421536865248e-8, (12, 2) = -0.2148435088658628e-8, (12, 3) = 0.1007706561001513e-7, (12, 4) = -0.7061367623274055e-8, (12, 5) = 0.9048774297591804e-8, (13, 1) = 0.67611797018740185e-8, (13, 2) = -0.13976774396261053e-8, (13, 3) = 0.8595435301400648e-8, (13, 4) = -0.4983619985985624e-8, (13, 5) = 0.6526712448285493e-8, (14, 1) = 0.6440901841297378e-8, (14, 2) = -0.7706065262608924e-9, (14, 3) = 0.741967636293521e-8, (14, 4) = -0.3217163270911993e-8, (14, 5) = 0.454060837838052e-8, (15, 1) = 0.6216116456491702e-8, (15, 2) = -0.320565246765447e-9, (15, 3) = 0.6605099271231325e-8, (15, 4) = -0.18085873364690693e-8, (15, 5) = 0.3117357572183856e-8, (16, 1) = 0.6064304046619167e-8, (16, 2) = -0.6701463087834589e-10, (16, 3) = 0.6136443374376152e-8, (16, 4) = -0.7630971430389972e-9, (16, 5) = 0.2232644899187691e-8, (17, 1) = 0.5954749976109993e-8, (17, 2) = .0, (17, 3) = 0.5954749976109993e-8, (17, 4) = .0, (17, 5) = 0.17510675834850087e-8}, 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[17] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(2.661709248184363e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [5, 17, [F(eta), diff(F(eta), eta), f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[17] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[17] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(5, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(17, 5, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(5, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(17, 5, X, Y, outpoint, yout, L, V) end if; [eta = outpoint, seq('[F(eta), diff(F(eta), eta), f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta)]'[i] = yout[i], i = 1 .. 5)] 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[17] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(2.661709248184363e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [5, 17, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[17] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[17] 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(5, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(17, 5, 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(5, {(1) = 0., (2) = 0., (3) = 0., (4) = 0., (5) = 0.}); `dsolve/numeric/hermite`(17, 5, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 5)] end proc, (2) = Array(0..0, {}), (3) = [eta, F(eta), diff(F(eta), eta), f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta)], (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); [eta = res[1], seq('[F(eta), diff(F(eta), eta), f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta)]'[i] = res[i+1], i = 1 .. 5)] catch: error  end try end proc

 

[eta = 2., F(eta) = HFloat(0.6659115637409204), diff(F(eta), eta) = HFloat(0.021498918991851815), f(eta) = HFloat(0.6365649907511923), diff(f(eta), eta) = HFloat(0.058915598688308896), diff(diff(f(eta), eta), eta) = HFloat(-0.07584262138836992)]

(2)

K := proc(s)
  local Kf, KF, Kdf, KdF:
  if not type(evalf(s),'numeric') then
    'procname'(s);
  else
    return eval([f(eta), F(eta), diff(f(eta), eta), diff(F(eta), eta)], sol__12(s))
  end if:
end proc:
  
# Check
plot(['K(s)[1]', 'K(s)[2]'], s=1..5, color=[red, blue], legend=['f(eta)', 'F(eta)'], labels=[eta, ""]);
 

 

sys__34 := isolate(eq3, diff(theta(eta), eta, eta)),
           isolate(eq4, diff(theta__p(eta), eta)):
sys__34 := eval([sys__34, theta(0) = 1, theta(5) = 0, theta__p(5) = 0], params);

[diff(diff(theta(eta), eta), eta) = .9090909091*(-6.2*(diff(theta(eta), eta))*f(eta)+12.4*(diff(f(eta), eta))*theta(eta)-.1*(diff(theta(eta), eta))-1.550*theta__p(eta)-4.650*theta(eta)-.1550*(diff(F(eta), eta))+.1550*(diff(f(eta), eta)))/(.2*eta+1), diff(theta__p(eta), eta) = (1/2)*(F(eta)*theta__p(eta)-.5*theta__p(eta)+.5*theta(eta))/f(eta), theta(0) = 1, theta(5) = 0, theta__p(5) = 0]

(3)

sys__34 := eval(sys__34, [f(eta)='K(eta)[1]', F(eta)='K(eta)[2]', diff(f(eta), eta)='K(eta)[3]', diff(F(eta), eta)='K(eta)[4]']);

[diff(diff(theta(eta), eta), eta) = .9090909091*(-6.2*(diff(theta(eta), eta))*K(eta)[1]+12.4*K(eta)[3]*theta(eta)-.1*(diff(theta(eta), eta))-1.550*theta__p(eta)-4.650*theta(eta)-.1550*K(eta)[4]+.1550*K(eta)[3])/(.2*eta+1), diff(theta__p(eta), eta) = (1/2)*(K(eta)[2]*theta__p(eta)-.5*theta__p(eta)+.5*theta(eta))/K(eta)[1], theta(0) = 1, theta(5) = 0, theta__p(5) = 0]

(4)

sol__34 := dsolve(sys__34, numeric, output=procedurelist, known=K, method=bvp[middefer], maxmesh=512, abserr=1e-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) = .0, (2) = .3370404083306422, (3) = .7008767283468323, (4) = 1.1112972869729818, (5) = 1.609997064526372, (6) = 2.189609762048589, (7) = 2.812077057741318, (8) = 3.44210959171944, (9) = 4.073127629221379, (10) = 4.627447082063739, (11) = 5.0}, datatype = float[8], order = C_order); Y := Matrix(11, 3, {(1, 1) = 1.0, (1, 2) = -2.631180930496676, (1, 3) = -.5069651588422761, (2, 1) = .4023510464572303, (2, 2) = -.9152718910087276, (2, 3) = -.10719058799533322, (3, 1) = .1726623490489365, (3, 2) = -.34732166973349327, (3, 3) = -0.3642234782724219e-1, (4, 1) = 0.7245529317024713e-1, (4, 2) = -.14099234754733977, (4, 3) = -0.11955501922042339e-1, (5, 1) = 0.23698556314753913e-1, (5, 2) = -0.5454307897597913e-1, (5, 3) = -0.20139753068670863e-2, (6, 1) = 0.3049194442117786e-2, (6, 2) = -0.16709196762319765e-1, (6, 3) = 0.10235425366078611e-2, (7, 1) = -0.2719856786139026e-2, (7, 2) = -0.18268814161701507e-2, (7, 3) = 0.11518671538119254e-2, (8, 1) = -0.2614649585026932e-2, (8, 2) = 0.2160855283988688e-2, (8, 3) = 0.600072097033077e-3, (9, 1) = -0.12708908932571583e-2, (9, 2) = 0.20981630384967246e-2, (9, 3) = 0.18102423815965176e-3, (10, 1) = -0.3515015852953094e-3, (10, 2) = 0.12190191503928644e-2, (10, 3) = 0.23439573975186445e-4, (11, 1) = .0, (11, 2) = 0.6679696148136223e-3, (11, 3) = .0}, datatype = float[8], order = C_order); YP := Matrix(11, 3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (2, 1) = -.9152718910087276, (2, 2) = 2.4727536746209924, (2, 3) = .3760823990110259, (3, 1) = -.34732166973349327, (3, 2) = .7893743028047262, (3, 3) = 0.9662576107625089e-1, (4, 1) = -.14099234754733977, (4, 2) = .2662816956052928, (4, 3) = 0.3229410448051668e-1, (5, 1) = -0.5454307897597913e-1, (5, 2) = 0.9768863241923123e-1, (5, 3) = 0.9509064038998313e-2, (6, 1) = -0.16709196762319765e-1, (6, 2) = 0.37412194655996434e-1, (6, 3) = 0.13133552399215123e-2, (7, 1) = -0.18268814161701507e-2, (7, 2) = 0.1163249008646273e-1, (7, 3) = -0.867516162778065e-3, (8, 1) = 0.2160855283988688e-2, (8, 2) = 0.14317512514059364e-2, (8, 3) = -0.8874409868549002e-3, (9, 1) = 0.20981630384967246e-2, (9, 2) = -0.15103520383209646e-2, (9, 3) = -0.44321926460382834e-3, (10, 1) = 0.12190191503928644e-2, (10, 2) = -0.16528238604940286e-2, (10, 3) = -0.12586527735752815e-3, (11, 1) = 0.6679696148136223e-3, (11, 2) = -0.13132690311376787e-2, (11, 3) = .0}, 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) = .0, (2) = .3370404083306422, (3) = .7008767283468323, (4) = 1.1112972869729818, (5) = 1.609997064526372, (6) = 2.189609762048589, (7) = 2.812077057741318, (8) = 3.44210959171944, (9) = 4.073127629221379, (10) = 4.627447082063739, (11) = 5.0}, datatype = float[8], order = C_order); Y := Matrix(11, 3, {(1, 1) = .0, (1, 2) = -0.13944241604693985e-1, (1, 3) = 0.2199320251397907e-2, (2, 1) = 0.3122582429547549e-1, (2, 2) = -.10070764912778984, (2, 3) = -0.10906300241481343e-1, (3, 1) = 0.2396545586162674e-1, (3, 2) = -0.8321640362333171e-1, (3, 3) = 0.3933131209080537e-3, (4, 1) = 0.8640618859978494e-2, (4, 2) = -0.4791621713263419e-1, (4, 3) = 0.3229350888215585e-2, (5, 1) = -0.17331011525628486e-2, (5, 2) = -0.1902105958328692e-1, (5, 3) = 0.32869341286691613e-2, (6, 1) = -0.4847665589688804e-2, (6, 2) = -0.29053410124298754e-2, (6, 3) = 0.21332021422819417e-2, (7, 1) = -0.36149600965167117e-2, (7, 2) = 0.197695795668519e-2, (7, 3) = 0.9832255075634968e-3, (8, 1) = -0.17708414039592315e-2, (8, 2) = 0.2040124287083247e-2, (8, 3) = 0.32939453342591216e-3, (9, 1) = -0.60979066315007e-3, (9, 2) = 0.11477551977751616e-2, (9, 3) = 0.680826223469885e-4, (10, 1) = -0.13922033932913127e-3, (10, 2) = 0.5435456008739585e-3, (10, 3) = 0.5485745412833969e-5, (11, 1) = .0, (11, 2) = 0.29761609702088786e-3, (11, 3) = .0}, 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 2 elif outpoint = "error" then return .100707649127789844 elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [3, 11, [theta(eta), diff(theta(eta), eta), theta__p(eta)], 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(3, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(11, 3, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(3, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(11, 3, X, Y, outpoint, yout, L, V) end if; [eta = outpoint, seq('[theta(eta), diff(theta(eta), eta), theta__p(eta)]'[i] = yout[i], i = 1 .. 3)] 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 2 elif outpoint = "error" then return .100707649127789844 elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [3, 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 ) = (true), ( 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(3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(11, 3, 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 ) = (true), ( 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(3, {(1) = 0., (2) = 0., (3) = 0.}); `dsolve/numeric/hermite`(11, 3, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 3)] end proc, (2) = Array(0..0, {}), (3) = [eta, theta(eta), diff(theta(eta), eta), theta__p(eta)], (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); [eta = res[1], seq('[theta(eta), diff(theta(eta), eta), theta__p(eta)]'[i] = res[i+1], i = 1 .. 3)] catch: error  end try end proc

(5)

display(
  odeplot(sol__34, [eta, theta(eta)], eta=0..5, color=red, legend='theta(eta)'),
  odeplot(sol__34, [eta, theta__p(eta)], eta=0..5, color=blue, legend='theta__p(eta)')
)

 

 


 

Download MapleOde_mmcdara.mw

According to my previous reply (correct expression of bcs and ics).
You 'll see in the attached file that there is no more RootOf.

elctroWav_mmcdara.mw

I'm ,evertheless a little bit surprised by your boundary conditions.
What are R and rb?
It seems that rho varies between 0 and 1 and that R-rb is the inner radius and R the outer radius of a tube within which a wave propagates?
Why don't set dyrectly bcs at rho=R-rb and rho=R?
 

restart
interface(version)
Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 

   21 2015 Build ID 1097895
p := .5:
n := 10:
nseq := 2:
with(Statistics):
randomize():
UseHardwareFloats := false:
x1 := seq(Sample(RandomVariable(BernoulliDistribution(p)), n), i = 1 .. nseq):
x2 := seq(convert(x1[i], list), i = 1 .. nseq);
        [0., 1.0, 0., 0., 1.0, 0., 0., 0., 1.0, 1.0], 

          [1.0, 1.0, 1.0, 1.0, 0., 1.0, 0., 0., 0., 0.]
numboccur(x2[1], {0.})
                               6

An alternative to numboccur:
 

Tally(x2[1])
                       [0. = 6, 1.0 = 4]

 

How to find the answer to question c with MAPLE?
 

restart:
with(IntegrationTools):
J16 := Int(f(x), x=1..6);
J16 = Split(J16, 4);
isolate(%, Int(f(x), x=4..6));
eval(%, [Int(f(x), x=1..4) = 7, Int(f(x), x=1..6) = 4])

                         /6             
                        |               
                        |   f(x) dx = -3
                        |               
                       /4               


Extra: Answer to question f could be obtained this way

Qf := Int(3*f(x)-6*g(x), x=1..6):
Qf = Expand(Qf):
lhs(%) = eval(rhs(%), [Int(f(x), x=1..6) = 4, Int(g(x), x=1..6) = 5]);
                   /6                         
                  |                           
                  |   3 f(x) - 6 g(x) dx = -18
                  |                           
                 /1                           

 

@VK057 

Here are two methods, probably amid several others:

restart:

f := 2.0*cos(1.5*x + 4.0) + 2.2:

g := 2.0*cos(1.6*x + 3.85) + 2.0:

 

# Method 1

r1 := fsolve(f=g, x=-1..0);
r2 := fsolve(f=g, x= 0..4);

 

-.7967442905

 

3.834531622

(1)

# Method 2

eq := unapply(f-g, x);
r1 := RootFinding:-NextZero(x -> eq(x), -1);
r2 := RootFinding:-NextZero(x -> eq(x), r1);

proc (x) options operator, arrow; 2.0*cos(1.5*x+4.0)+.2-2.0*cos(1.6*x+3.85) end proc

 

-.7967442905

 

3.834531622

(2)

Download TwoMethods.mw

More generaly, if you want to select all the roots of f-g within some interval [x1, x2] (here [-8, 8]) you can do this

x1 := -8:
x2 :=  8:
r  := NULL:
while x1 < x2 do
  x1 := RootFinding:-NextZero(x -> eq(x), x1);
  r  := r, x1
end do:
r := select(`<`, [r], x2)
   [-6.502888690, -4.665889115, -2.362208506, -0.7967442905, 

     3.834531622, 5.406998312, 7.708001175]

Consider this loop

for i from 0 to N do:  
dsolve({equ[1][i],cond[1][i]},X[i](t));  
X[i](t):=rhs(%):  
end do;

and look to the ode you solve:

i:=1:
{equ[1][i],cond[1][i]}
      // d         \                                       
     { |--- X[1](t)| + 2 X[0](t) Y[0](t) = 0, X[1](0) = 0, 
      \\ dt        /                                       

                                                  \ 
       X[1](5) = 0, D(X[1])(0) = 0, D(X[1])(5) = 0 }
                                                  / 

This is a first order ode with 4 boundary conditions while only one should be given.
More of this the boundary condition cannot be a condition on the derivative of X[1](t).
As a result dsolve returns "nothing" and the command

rhs(%)

generates an error.

It is up to you to know what is the correct boundary condition 

X[1](0) = 0;
#or
X[1](5) = 0;

 

First point: 
A simple way to generate your system of equations is:
 

N   := 4:   # or 9 in your last reply
P   := 2*N:
sys := [ add(W__||n, n=1..N)=4, seq(add(W__||n*(zeta__||n)^(p-1), n=1..N)=2/p*modp(p, 2), p=2..P) ]:
print~(sys):  # check


Second point;
The coefficients 2/3, 2/5, ..., 2/(2*p+1) are (curiously ?) the values of the integral of the pth Legendre polynomial squared over -1 to 1:

seq([p, int(simplify(LegendreP(p, x),'LegendreP')^2, x=-1..1)], p=1..6)
        [   2]  [   2]  [   2]  [   2]  [   2 ]  [   2 ]
        [1, -], [2, -], [3, -], [4, -], [5, --], [6, --]
        [   3]  [   5]  [   7]  [   9]  [   11]  [   13]

Your set of equations has an an extremely strong resemblance (apart for the first relation W__1+...+W__N=4 instead of the "normal" W__1+...+W__N=2 relation) with the set of relations that enable computing the weights (W) and the locations (zeta) of the Gauss-Legendre points.
Is this a coincidental resemblance ?

For instance this code returns the weights and locations of the Gauss-Legendre points for a N points Gauss quadrature.

N   := any_strictly_positive_integer:
P   := 2*N:
sys := [ add(W__||n, n=1..N)=2, seq(add(W__||n*(zeta__||n)^(p-1), n=1..N)=2/p*modp(p, 2), p=2..P) ]:
print~(sys): #check
solve(sys);


These same points can be computed in a far more efficient way for large values of N
(see for instance Gauss–Legendre_quadrature or 9403469.pdf)

Thus, if your  "first" equation W__1+...+W__N=4 is indeed the correct one, one can imagine that the methods mentioned in this link could be adjusted to solve your specific system.


Third point;
If the "idea" suggested at the end of the previous point is a dead end, one can simplify your problem by assuming that there exist symmetries of "Gauss-Legendre points".
For instance, without loss of generality, and if N=4, on can assume that 

W__1 = W__4, W__2 = W__3, zeta__1 = -zeta__4, zeta__2 = -zeta__3

Next one can infer that W > 0 and that abs(zeta) <=1 (for the rhs of relations of odd ranks are decreasing as the powers of zeta are increasing).

For instance:

N   := 6:
P   := 2*N:
sys := [ add(W__||n, n=1..N)=N, seq(add(W__||n*(zeta__||n)^(p-1), n=1..N)=2/p*modp(p, 2), p=2..P) ]:

t0  := time():
sol := fsolve(sys): # no solution found
time()-t0;
                             16.204

#----------------------------------------------
hypothesis := [seq(W__||(N-i+1)=W__||i, i=1..N/2), seq(zeta__||(N-i+1)=-zeta__||i, i=1..N/2)]:
redsys     := eval(sys[[seq(1..P, 2)]], hypothesis):

t0  := time():
redsol := fsolve(redsys): # no solution found
time()-t0;

                             2.753
#----------------------------------------------
t0  := time():
redsol2 := fsolve(redsys, {seq(W__||i=0..4, i=1..N/2), seq(zeta__||i=-1..1, i=1..N/2)});
time()-t0;
{W__1 = 2.432885814, W__2 = 0.1793486330, W__3 = 0.3877655527, 

  zeta__1 = -0.08695232853, zeta__2 = -0.9294284522, 

  zeta__3 = -0.6423774108}
                             0.025


NO MIRACLE HERE: for N>6 no solution is found and some work remains to be done
For instance :

  1. Transform the problem into a minimization one.
  2. When N is odd, forcing the "central" zeta to 0  (zeta__i where i=ceil(N/2))
  3. Taking into account that for the values of zeta for N-1 and N are alternate in the sense that
    • -1 < zeta__1 < zeta[N-1, 1]
    • zeta[N-1, , k] < zeta__k < zeta[N-1, k+1] for k=2..N-1
    • zeta[N-1, N-1] < zeta__N < 1

This could provide better search ranges for the zeta

Here is procedure aimed to compute W and zeta for relatively large values of N which implements these ideas.
It returns:

  • the solution as a list of equalities of the form W__i =...  and  zeta__i = ...
  • the list of absolute errors when this solution is plugged into the system to solve.

The second argument in the Gale procedure is a switch to select the "natural" -1..1 range for zeta, or the "potentially optimal" ranges defined in point 3 above.
You will see in the attached file that the natural range option gives var better results.

GaLe.mw

restart

GaLe := proc(M, zoom)
  local N, P, sys, hypothesis, redsys, constr, sol, nm, xi, here, J, abserr:

  Digits := 20:
  N          := min(M, 6):
  P          := 2*N:
  sys        := [ add(W__||n, n=1..N)=4, seq(add(W__||n*(zeta__||n)^(p-1), n=1..N)=2/p*modp(p, 2), p=2..P) ]:
  hypothesis := [seq(W__||(N-i+1)=W__||i, i=1..N/2), seq(zeta__||(N-i+1)=-zeta__||i, i=1..N/2)];
  redsys     := eval(sys[[seq(1..P, 2)]], hypothesis):
  if N::odd then
    redsys := eval(redsys, zeta__||(ceil(N/2))=0):
  end if:
  constr     := {seq(W__||i=0..4, i=1..ceil(N/2)), seq(zeta__||i=-1..1, i=1..N/2)}:
  sol        := fsolve(evalf(redsys), constr);

  if M > 6 then
    for nm from N+1 to M do
      xi   := sort(eval([seq(zeta__||i, i=1..nm/2)], sol)):
      here := [-1..xi[1], seq(xi[k]..xi[k+1], k=1..nops(xi)-1)];
      if nm::even then here := [here[], op(2, here[-1])..0]: end if:

      P          := 2*nm:
      sys        := [ add(W__||n, n=1..nm)=nm, seq(add(W__||n*(zeta__||n)^(p-1), n=1..nm)=2/p*modp(p, 2), p=2..P) ]:
      hypothesis := [seq(W__||(nm-i+1)=W__||i, i=1..nm/2), seq(zeta__||(nm-i+1)=-zeta__||i, i=1..nm/2)];
      redsys     := eval(sys[[seq(1..P, 2)]], hypothesis):

      if nm::odd then
        redsys := eval(redsys, zeta__||(ceil(nm/2))=0):
        hypothesis := [hypothesis[], zeta__||(ceil(nm/2))=0];
      end if:

      J      := evalf(add(((lhs-rhs)~(redsys))^~2)):
      if zoom then
        constr := {
                    seq(W__||i >= 0, i=1..ceil(nm/2)),
                    seq(W__||i <= 4, i=1..ceil(nm/2)),
                    seq(zeta__||i >= op(1, here[i]), i=1..floor(nm/2)),
                    seq(zeta__||i <= op(2, here[i]), i=1..floor(nm/2))
                  }:
      else
        constr := {
                    seq(W__||i >= 0, i=1..ceil(nm/2)),
                    seq(W__||i <= 4, i=1..ceil(nm/2)),
                    seq(zeta__||i >= -1, i=1..floor(nm/2)),
                    seq(zeta__||i <= +1, i=1..floor(nm/2))
                  }:
      end if:
      sol := Optimization:-Minimize(J, constr, iterationlimit=1000)[2];
      if nm::odd then sol := [sol[], zeta__||(ceil(nm/2))=0] end if:
    end do:
    sol := [sol[], eval(hypothesis, sol)[]];
  end if:  
  if M <= 6 then
    sol := [sol[], eval(hypothesis, sol)[], `if`(M::odd, zeta__||(ceil(M/2))=0, NULL)];
  end if:
  abserr := abs~(evalf(eval((lhs-rhs)~(sys), sol))):
  return sol, abserr:
end proc:

# 2nd argument of procedure Gale:
#    false means zeta is searched between -1 and +1
#    true  means zeta is searched in the interval defined by the values obtained for N-1

N := 9:
sol_F := CodeTools:-Usage( GaLe(N, false) ):

print():

sol_T := CodeTools:-Usage( GaLe(N, true) ):

memory used=348.60MiB, alloc change=36.00MiB, cpu time=2.03s, real time=4.70s, gc time=150.48ms

 

 

memory used=140.04MiB, alloc change=0 bytes, cpu time=811.00ms, real time=814.00ms, gc time=72.94ms

 

dataplot(
  <[sol_F[2], sol_T[2]]>,
  color=[blue, red],
  legend=["False", "True"],
  labels=["eq number", "abserr"],
  labeldirections=[default, vertical]
);

 

ix_F    := sort( eval([seq(zeta__||i, i=1..N)], sol_F[1]), output=permutation):
zeta_FS := [ seq(eval(zeta__||(ix_F[i-N]), sol_F[1]), i=N+1..2*N) ]:
W_FS    := sort(  eval([seq(W__||i, i=1..N)], sol_F[1]) ):
W_FS    := [ W_FS[[seq(1..N, 2)]][], W_FS[[seq(N..1, -2)]][] ]:

#ix_T := sort( eval([seq(zeta__||i, i=1..N)], sol_F[1]), output=permutation):

sol_TS := sort(sol_T[1]):

Matrix(
  2*N, 3,
  (i,j) -> `if`(
             j=1,
             `if`(i <= N, W__||i, zeta__||(i-N)),
             `if`(
               j=2,
               `if`(i <= N, W_FS[i], zeta_FS[N-i]),
               `if`(i <= N, eval(W__||i, sol_TS), eval(zeta__||(i-N), sol_TS))
             )
           )
):
evalf[6](%):
< `<|>`(["", "False", "True"]), % >:

# DocumentTools:-Tabulate(%, width=30, alignment=left):

plots:-display(
  seq(
    plot([[zeta_FS[i], -0.1], [zeta_FS[i], 0.1]], thickness=7, color=blue, transparency=0.7, legend="False"),
    i=1..N
  ),
  seq(
    plot([[eval(zeta__||i, sol_T[1]), -0.1], [eval(zeta__||i, sol_T[1]), 0.1]], thickness=3, color=red, legend="True"),
    i=1..N
  ),
  view=[-1..1, -0.1..0.1],
  axis[2]=[tickmarks=[]],
  scaling=constrained
)

 

# Legendre polynomials compared to your "pesudo-Legendre" polynomials

Matrix(
  3, 3,
  (m,n) ->
    plots:-dualaxisplot(
      plot(simplify(LegendreP(3*(m-1)+n, x),'LegendreP'), x=-1..1, color=blue, legend="Legendre"),
      plot(
        expand(eval(mul((x-zeta__||i), i=1..3*(m-1)+n), GaLe(3*(m-1)+n, true)[1])),
        x=-1..1, color=red, legend="pseudo Legendre"
      ),
      title=cat("N = ", 3*(m-1)+n)
    )
):

# DocumentTools:-Tabulate(%, width=60):

 

 

 

 

Export your plot, either by clicking on it or programatically (see plotsetup), in jpeg, png or pdf, and insert thisfile into your docx document

You want to solve a system of 4 linear equations depending on 7 parameters (a, b, c, t, x, y, z) by condidering it as a system of 4 linear equations in 3 unknowns (a, b, c), parameterized by (t, x, y, z).
A priori it's not surprising that you got no solutions.

Try any of these instead:
 

solve({seq(G[i] = 0, i = 1 .. 4)});
{a = 2 y + 2 z, b = -2 x + 2 y - c, c = c, t = x - y - z, x = x, 

  y = y, z = z}
for u in [t, x, y, z] do
solve({seq(G[i] = 0, i = 1 .. 4)}, {a, b, c, u});
end do;
   {a = 2 y + 2 z, b = -2 x + 2 y - c, c = c, t = x - y - z}
   {a = 2 y + 2 z, b = -c - 2 z - 2 t, c = c, x = y + z + t}
   {a = -2 t + 2 x, b = -c - 2 z - 2 t, c = c, y = x - z - t}
   {a = -2 t + 2 x, b = -2 x + 2 y - c, c = c, z = x - y - t}


Other solution, it it makes sense to you, solve your system in the "least squares sense" :

vars := [a, b, c]:
M,V := GenerateMatrix(sys, vars):
<vars> = (M^+.M).(M^+.V);

Just do this
 

# write the general solutions
Un := (x, t, n) -> T(t, n)*X(x, n);

# Derivative of Un w.r.t t
DtUn := eval(diff(Un(x, t, n), t), t=0):

# Compute the Fourier coefficients
int(DtUn*sin(Pi*n*x/l), x=0..l):
int(psi*sin(Pi*n*x/l), x=0..l) assuming l::real:
isolate(%%=%, C1[n]);
    
                                  2 /                 /1     \\
          piecewise(l < 0, 0, 1) l  |sin(Pi n) - 2 sin|- Pi n||
                                    \                 \2     //
C1[n] = - -----------------------------------------------------
                2  2 /  1                         1       \    
              Pi  n  |- - a cos(Pi n) sin(Pi n) + - a Pi n|    
                     \  2                         2       /    

 

If you realy want to use the notation e^(...) instead of the "normal" one exp(...) you can do this

restart:
local e := exp(1):
uu := (x, y, t) -> e^(.2*t)*(e^(-x)+e^(-y)):
uu(.1, .1, .1);
                               2      
                          ------------
                                  0.08
                          (exp(1))    
evalf(%)
                          1.846232693



Or, if you are only interested in floating point values 

restart:
local e := exp(1.):
uu := (x, y, t) -> e^(.2*t)*(e^(-x)+e^(-y));
uu(.1, .1, .1)
                          1.846232693

But listening acer is a far better idea

EDITED  10/27 9:30 pm GMT

For a reason I don't understand it seems that you must transform your sine expression before using solve.

restart:
f := sin(2*t)-sin(3*t)+sin(4*t);     
g := expand(f, trig):
sol_g := solve([g, `and`(t >= 0, t <= 2*Pi)], allsolutions, explicit=true)

 /    1   \    /    5   \    /    2   \    /    4   \            
{ t = - Pi }, { t = - Pi }, { t = - Pi }, { t = - Pi }, {t = 0}, 
 \    3   /    \    3   /    \    3   /    \    3   /            

  {t = Pi}, {t = 2 Pi}

 

restart
ode := 1/r*diff(r*u(r), r)=2*a:
dsolve(ode, u(r));
                                 2      
                              a r  + _C1
                       u(r) = ----------
                                  r     

You can also use expand(%) after the dsolve command if you prefer

2 real solutions and 4 complex ones :
 

restart:

e1 := x^2+2*m*x-6*x-m^3;
e2 := (x-r)*(x-r^2);
sol := solve([coeffs(expand(e1-e2), x)], [m, r])

-m^3+2*m*x+x^2-6*x

 

(x-r)*(-r^2+x)

 

[[m = 2, r = -2], [m = -3, r = 3], [m = -(1/2)*RootOf(_Z^4+4*_Z^3-5*_Z^2-24*_Z+36)^2-(1/2)*RootOf(_Z^4+4*_Z^3-5*_Z^2-24*_Z+36)+3, r = RootOf(_Z^4+4*_Z^3-5*_Z^2-24*_Z+36)]]

(1)

# check
map(s -> simplify(eval(e1-e2, s)), sol)

[0, 0, 0]

(2)

simplify([allvalues(sol[3])])

[[m = 1/4+((3/4)*I)*3^(1/2)+(1/4)*(25-(4*I)*3^(1/2))^(1/2)-((1/4)*I)*3^(1/2)*(25-(4*I)*3^(1/2))^(1/2), r = -1+((1/2)*I)*3^(1/2)+(1/2)*(25-(4*I)*3^(1/2))^(1/2)], [m = 1/4+((3/4)*I)*3^(1/2)-(1/4)*(25-(4*I)*3^(1/2))^(1/2)+((1/4)*I)*3^(1/2)*(25-(4*I)*3^(1/2))^(1/2), r = -1+((1/2)*I)*3^(1/2)-(1/2)*(25-(4*I)*3^(1/2))^(1/2)], [m = 1/4-((3/4)*I)*3^(1/2)-(1/4)*(25+(4*I)*3^(1/2))^(1/2)-((1/4)*I)*3^(1/2)*(25+(4*I)*3^(1/2))^(1/2), r = -1-((1/2)*I)*3^(1/2)-(1/2)*(25+(4*I)*3^(1/2))^(1/2)], [m = 1/4-((3/4)*I)*3^(1/2)+(1/4)*(25+(4*I)*3^(1/2))^(1/2)+((1/4)*I)*3^(1/2)*(25+(4*I)*3^(1/2))^(1/2), r = -1-((1/2)*I)*3^(1/2)+(1/2)*(25+(4*I)*3^(1/2))^(1/2)]]

(3)

# check
map(s -> simplify(eval(e1-e2, s)), %)

[0, 0, 0, 0]

(4)

 


 

Download AllSolutions.mw

restart
phi := (1+sqrt(5))/2;

_phi := `#mi("&phi;")`;


plots[pointplot]([seq([n, sin(n*phi)], n = 1000 .. 2000)], symbol = point, axes = boxed, labels = [n, sin(_phi*n)], labeldirections = [horizontal, vertical])


(obtained with Maple 2015.2)

2 3 4 5 6 7 8 Last Page 4 of 28