Advanced Search

Rouben Rostamian

7711 Reputation

22 Badges

19 years, 3 days
Contact Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

Ctrl+U...

Rouben Rostamian 7711
July 05 2019
0 0

You may view (but not edit) the HTML code of any web page (not only MaplePrimes) by pressing Ctrl and U keys simultaneously.  This works on Firefox and Chromebrowser.  I suspect that it will work similarly on other browsers.

Beware that the code shown is that of the complete web page, including the navigation bar and the side boxes.  The text of your message would be a very small part of the whole thing.  You may search for an identifiable string within that page to locate your message within it.  In viewing the source code of your current message, I searched for "observe" and located the message part immediately.  It reads:

<span id="ctl00_MainContent_body" class="mainBody document"><p>When writing a post on MaplePrimes it is possible to access the source code, which allows to introduce/edit many attributes. Is it possible to observe the source of another post?</p>

 

Solution and comments...

Rouben Rostamian 7711
July 04 2019
1 3

While I was writing my solution, Carl Love posted his.  I am posing my writeup anyway since its comments may be helpful to you.

You have this initial-boundary-value problem in the unknown u(t, x):
"(&PartialD; u)/(&PartialD; t)=-(&PartialD; u)/(&PartialD; x)+g(t,x)",
"u(0,x)=f(x), u(t,0) = alpha(t),"
which you wish to solve it through Picard iterations.  To do that, we
integrate the PDE with respect to t

"u(t,x)=f(x)-(&int;)[0]^(t)(&PartialD; u(s,x))/(&PartialD; x)- g(s,x) &DifferentialD;s, "
and then set up the following two-step Picard iteration scheme
step 1:   u__n(t, x) = f(x)-(int(`&PartialD;`(`u__n-1`(s, x))/`&PartialD;`(x)-g(s, x), s = 0 .. t))
step 2:   u__n(t, x) = u__n(t, x)-u__n(t, 0)+alpha(t)
with "n=1,2,..., "and an arbitrary initial guess u__0(x).

Important!  This iteration scheme implies that "`u__n`(t,0)=alpha(t),"
which is good, but u__n(0, x) = f(x)-f(0)+alpha(0), which is not good.
Getting the desired initial condition u__n(0, x) = f(x) requires
the compatibility condition f(0) = alpha(0) to hold!

The proc doit()  implements the iteration scheme.  It receives a u
and returns the next u.

> 

restart;
> 

doit := proc(u::procedure)
  # next line needs to be fixed.  See note at the bottom.
  f(x) - int(D[2](u[n-1])(s,x) - g(s,x), s=0..t);
  eval(%, x=0) - alpha(t);
  return unapply(%%-%, [t,x]);
end proc:

Let's apply this to your specific problem with an arbitrarily picked

initial guess of u__0(t, x) = t*x^4.

> 

g := (t,x) -> 2 + t + x;

proc (t, x) options operator, arrow; 2+t+x end proc

> 

f := x -> 1 + x;

proc (x) options operator, arrow; x+1 end proc

> 

alpha := t -> 1 + t;

proc (t) options operator, arrow; t+1 end proc

Verify compatibility:

> 

is(f(0)=alpha(0));

true

Now iterate:

> 

u[0] := (t,x) -> t*x^4;
for n from 1 to 5 do
  u[n] := doit(u[n-1]);
end do;

proc (t, x) options operator, arrow; t*x^4 end proc

proc (t, x) options operator, arrow; x+1-(1/2)*(4*x^3-1)*t^2+x*t+t-(1/2)*t^2 end proc

proc (t, x) options operator, arrow; 2*t^3*x^2+x*t+t+x+1 end proc

proc (t, x) options operator, arrow; -t^4*x+x*t+t+x+1 end proc

proc (t, x) options operator, arrow; x*t+t+x+1 end proc

proc (t, x) options operator, arrow; x*t+t+x+1 end proc

In this particular case the iteration arrives at the exact answer
at step u__4 since u__5 = u__4.

You should be able to verify by hand that u__4 satisfies the PDE as well
as the initial and boundary conditions.

 

Note added later:
The line printed in green is wrong -- it produces the correct result purely by accident.  It should have been:
f(x) - int(D[2](u)(s,x) - g(s,x), s=0..t);
You will need to fix that line in the attached worksheet.

 

Download picard-for-pde.mw

Solving a PDE with a moving boundary...

Rouben Rostamian 7711
July 02 2019
6 8
> 

restart;
> 

with(PDEtools):

The original PDE

> 

pde := diff(c(z,t),t) + diff(w*c(z,t),z) = diff(c(z,t),z,z);

diff(c(z, t), t)+w*(diff(c(z, t), z)) = diff(diff(c(z, t), z), z)

The space domain is the shrinking interval w*t, h.

> 

ibc := { c(z,0)=0, c(w*t,t)=cs, D[1](c)(h,t)=0 };

{c(z, 0) = 0, c(w*t, t) = cs, (D[1](c))(h, t) = 0}

Change the space variable from z to Z, where Z lives in the fixed interval 0, 1.
The transformation is "z=Z*h + (1-Z)*w*t."
Then Z = 0 corresponds to z = w*t,  and Z = 1 corresponds to z = h.

> 

change_of_vars := {z=Z*h+(1-Z)*w*t, c(z,t)=C(Z,t)};

{z = Z*h+(1-Z)*w*t, c(z, t) = C(Z, t)}

Here is the transformed PDE

> 

PDE := dchange(change_of_vars, pde, [C(Z,t), Z]);

diff(C(Z, t), t)+w*(diff(C(Z, t), Z))/(-t*w+h) = (diff(diff(C(Z, t), Z), Z))/(-t*w+h)^2

and here are the transformed initial and boundary conditions:

> 

IBC := dchange(change_of_vars, ibc, [C(Z,t), Z]);

{(D[1](C))(1, t)/(-t*w+h) = 0, C(0, t) = cs, C(Z, 0) = 0}

I pick some arbitrary values for the parameters.  Change as needed.

> 

w := 0.5: h := 1:  cs := 1:
pdsol := pdsolve(PDE, IBC, numeric, spacestep=0.01):

Here is the plot of the concentration in the Z, t variables:

> 

p := pdsol:-plot3d(C(Z,t), Z=0..1, t=0..1, orientation=[-45,60,0]);

Apply the transformation z = Z*h+(1-Z)*w*t to obtain the plot of the
concentration in the originalz, t variables.

> 

plottools:-transform((Z,t,C)->[Z*h + (1-Z)*w*t, t, C])(p):
plots:-display(%, labels=[z,t,c],
    shading=zhue, style=surfacecontour, orientation=[-45,60,0]);


 

Download mw.mw

Solve it with pdsolve...

Rouben Rostamian 7711
July 01 2019
2 1

Consider an ODE of the form diff(x(t), t) = f(x(t), mu) for the unknown x(t), where mu is a parameter.
We wish to make a 3D plot of the solution x as a function of t and mu.

To do that, view the ODE as a PDE in the unknown x(t, mu), that is, "(&PartialD; x)/(&PartialD; t)=f(x,mu)".  Solve that
PDE, and then 3D-plot the solution "x(t,mu)."  Here is how it works.

> 

restart;

Your ODE is

> 

diff(B[1](t),t) = piecewise(t < 1000, kaC*(R-B[1](t)) - kd1*B[1](t), -kd1*B[1](t));

diff(B[1](t), t) = piecewise(t < 1000, kaC*(R-B[1](t))-kd1*B[1](t), -kd1*B[1](t))

Arrange the right-hand side of the ODE into the form A__1(t)+A__2(t)*B__1(t) because
that's what pdsolve() prefers.

> 

de := diff(B[1](t,kd1),t) = piecewise(t<1000, kaC*R,0)
                           -piecewise(t<1000, kaC+kd1, kd1)*B[1](t,kd1);

de := diff(B[1](t, kd1), t) = piecewise(t < 1000, kaC*R, 0)-piecewise(t < 1000, kaC+kd1, kd1)*B[1](t, kd1)

Pick your initial condition.  Change as needed..

> 

ic := B[1](0,kd1) = 0;

B[1](0, kd1) = 0

Fix the values of R and kaC:

> 

myde := eval(de, {R=1, kaC=6e-1});

myde := diff(B[1](t, kd1), t) = piecewise(t < 1000, .6, 0)-piecewise(t < 1000, .6+kd1, kd1)*B[1](t, kd1)

Solve the PDE:

> 

pdsol := pdsolve(myde, {ic}, numeric, time=0..2000, range=0..1);

_m139769764691840

> 

pdsol:-plot3d(t=0..2000, kd1=0..1);


 

Download mw.mw

Solution...

Rouben Rostamian 7711
June 28 2019
1 3

Looking at your maplee.docx file, it appears that you are
interested in solving

diff(f(eta), eta, eta, eta)+f(eta)*(diff(f(eta), eta, eta)) = 0, where f = f(eta),

with the boundary conditions
"(lim)f '(eta)=0",  " (lim)f '(eta)=1",     limit(-eta*f+eta, eta = infinity) = 0.

If so, then there is no reason to appeal to shooting.
Let Maple's dsolve() do the work.

> 

restart;
> 

de := diff(f(eta), eta$3) + diff(f(eta), eta$2)*f(eta) = 0;

diff(diff(diff(f(eta), eta), eta), eta)+(diff(diff(f(eta), eta), eta))*f(eta) = 0

> 

a := 10;  # surrogate for infinity

10

> 

dsol := dsolve({de, D(f)(-a)=0, D(f)(a)=1, f(a)=a}, numeric);

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(27, {(1) = -10.0, (2) = -9.232690925174811, (3) = -8.400948293356127, (4) = -7.568968590305362, (5) = -6.737727888889229, (6) = -5.907313057964531, (7) = -5.078204592534601, (8) = -4.251594150448721, (9) = -3.429391562535225, (10) = -2.6141287117372447, (11) = -1.8085471664745767, (12) = -1.0149443418684818, (13) = -.23460448462138273, (14) = .5322867133237555, (15) = 1.286245362574061, (16) = 2.0299285351172367, (17) = 2.7657746379629105, (18) = 3.4960588686445333, (19) = 4.222709593140783, (20) = 4.9471575201173215, (21) = 5.670371337996262, (22) = 6.392927283827333, (23) = 7.1151455505991725, (24) = 7.8371254600843265, (25) = 8.558285432902858, (26) = 9.279151313463982, (27) = 10.0}, datatype = float[8], order = C_order); Y := Matrix(27, 3, {(1, 1) = -.8738699745131914, (1, 2) = .0, (1, 3) = 0.19453618355416916e-3, (2, 1) = -.8737974379458486, (2, 2) = 0.2126559515304102e-3, (2, 3) = 0.38036413282775076e-3, (3, 1) = -.8734503794278492, (3, 2) = 0.6777103278353436e-3, (3, 3) = 0.7866519226334487e-3, (4, 1) = -.8725343958452793, (4, 2) = 0.1639670425821428e-2, (4, 3) = 0.16264464956168914e-2, (5, 1) = -.8704450947492267, (5, 2) = 0.3624927883294863e-2, (5, 3) = 0.335660042998554e-2, (6, 1) = -.8659405026702968, (6, 2) = 0.7710994731576229e-2, (6, 3) = 0.6904143252000777e-2, (7, 1) = -.8564897103613516, (7, 2) = 0.16073894873647625e-1, (7, 3) = 0.14106526134330544e-1, (8, 1) = -.8370205542219876, (8, 2) = 0.3299189209715141e-1, (8, 3) = 0.2843241617985271e-1, (9, 1) = -.7976730782373457, (9, 2) = 0.6644809848653324e-1, (9, 3) = 0.557806032742662e-1, (10, 1) = -.7203134591786324, (10, 2) = .12985832536945638, (10, 3) = .10392802789486236, (11, 1) = -.5747724484150437, (11, 2) = .2412044745700165, (11, 3) = .17615298103946722, (12, 1) = -.31941706728599945, (12, 2) = .41268706362224855, (12, 3) = .253453131738646, (13, 1) = 0.8460907160454192e-1, (13, 2) = .6265186144300298, (13, 3) = .2808159637909694, (14, 1) = .6433980071002114, (14, 2) = .8221742083467642, (14, 3) = .21447524251902275, (15, 1) = 1.3137106766127817, (15, 2) = .9417404576451969, (15, 3) = .10313792024850699, (16, 1) = 2.034509148074173, (16, 2) = .9876095153044471, (16, 3) = 0.2975926876438234e-1, (17, 1) = 2.7662734631043637, (17, 2) = .9983401352762191, (17, 3) = 0.5090030834695261e-2, (18, 1) = 3.4960935964279356, (18, 2) = .9998619837867148, (18, 3) = 0.5172370314613091e-3, (19, 1) = 4.222711112789611, (19, 2) = .9999929502456798, (19, 3) = 0.31293924121315814e-4, (20, 1) = 4.947157561890147, (20, 2) = .999999778465008, (20, 3) = 0.1136651162622266e-5, (21, 1) = 5.670371338671851, (21, 2) = .9999999959678437, (21, 3) = 0.23487052888232733e-7, (22, 1) = 6.392927283836079, (22, 2) = .9999999999485296, (22, 3) = 0.3247918498716828e-9, (23, 1) = 7.115145550598878, (23, 2) = 1.0000000000040963, (23, 3) = -0.2718938301063687e-10, (24, 1) = 7.837125460084901, (24, 2) = .9999999999982201, (24, 3) = 0.1158736274671473e-10, (25, 1) = 8.558285432902982, (25, 2) = 1.0000000000005311, (25, 3) = -0.5177298091061657e-11, (26, 1) = 9.27915131346414, (26, 2) = .9999999999995605, (26, 3) = 0.2484296054720659e-11, (27, 1) = 10.0, (27, 2) = 1.0, (27, 3) = -0.12649772507404847e-11}, datatype = float[8], order = C_order); YP := Matrix(27, 3, {(1, 1) = .0, (1, 2) = 0.19453618355416916e-3, (1, 3) = 0.16999932976437534e-3, (2, 1) = 0.2126559515304102e-3, (2, 2) = 0.38036413282775076e-3, (2, 3) = 0.33236120475138304e-3, (3, 1) = 0.6777103278353436e-3, (3, 2) = 0.7866519226334487e-3, (3, 3) = 0.6871014203018328e-3, (4, 1) = 0.1639670425821428e-2, (4, 2) = 0.16264464956168914e-2, (4, 3) = 0.1419130510427756e-2, (5, 1) = 0.3624927883294863e-2, (5, 2) = 0.335660042998554e-2, (5, 3) = 0.29217363793140588e-2, (6, 1) = 0.7710994731576229e-2, (6, 2) = 0.6904143252000777e-2, (6, 3) = 0.597857727814529e-2, (7, 1) = 0.16073894873647625e-1, (7, 2) = 0.14106526134330544e-1, (7, 3) = 0.12082094482997605e-1, (8, 1) = 0.3299189209715141e-1, (8, 2) = 0.2843241617985271e-1, (8, 3) = 0.23798516748730524e-1, (9, 1) = 0.6644809848653324e-1, (9, 2) = 0.557806032742662e-1, (9, 3) = 0.4449468551972008e-1, (10, 1) = .12985832536945638, (10, 2) = .10392802789486236, (10, 3) = 0.7486075727856172e-1, (11, 1) = .2412044745700165, (11, 2) = .17615298103946722, (11, 3) = .10124788020766334, (12, 1) = .41268706362224855, (12, 2) = .253453131738646, (12, 3) = 0.8095725603441037e-1, (13, 1) = .6265186144300298, (13, 2) = .2808159637909694, (13, 3) = -0.23759577988088584e-1, (14, 1) = .8221742083467642, (14, 2) = .21447524251902275, (14, 3) = -.13799294360907377, (15, 1) = .9417404576451969, (15, 2) = .10313792024850699, (15, 3) = -.13549338699410124, (16, 1) = .9876095153044471, (16, 2) = 0.2975926876438234e-1, (16, 3) = -0.6054550454113386e-1, (17, 1) = .9983401352762191, (17, 2) = 0.5090030834695261e-2, (17, 3) = -0.14080417224400456e-1, (18, 1) = .9998619837867148, (18, 2) = 0.5172370314613091e-3, (18, 3) = -0.18083090735272773e-2, (19, 1) = .9999929502456798, (19, 2) = 0.31293924121315814e-4, (19, 3) = -0.13214520114987515e-3, (20, 1) = .999999778465008, (20, 2) = 0.1136651162622266e-5, (20, 3) = -0.562319239439797e-5, (21, 1) = .9999999959678437, (21, 2) = 0.23487052888232733e-7, (21, 3) = -0.1331803115273048e-6, (22, 1) = .9999999999485296, (22, 2) = 0.3247918498716828e-9, (22, 3) = -0.2076370678612273e-8, (23, 1) = 1.0000000000040963, (23, 2) = -0.2718938301063687e-10, (23, 3) = 0.19345641755166163e-9, (24, 1) = .9999999999982201, (24, 2) = 0.1158736274671473e-10, (24, 3) = -0.9081161559751732e-10, (25, 1) = 1.0000000000005311, (25, 2) = -0.5177298091061657e-11, (25, 3) = 0.44308794834529397e-10, (26, 1) = .9999999999995605, (26, 2) = 0.2484296054720659e-11, (26, 3) = -0.23052158999194988e-10, (27, 1) = 1.0, (27, 2) = -0.12649772507404847e-11, (27, 3) = 0.12649772507404847e-10}, 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(27, {(1) = -10.0, (2) = -9.232690925174811, (3) = -8.400948293356127, (4) = -7.568968590305362, (5) = -6.737727888889229, (6) = -5.907313057964531, (7) = -5.078204592534601, (8) = -4.251594150448721, (9) = -3.429391562535225, (10) = -2.6141287117372447, (11) = -1.8085471664745767, (12) = -1.0149443418684818, (13) = -.23460448462138273, (14) = .5322867133237555, (15) = 1.286245362574061, (16) = 2.0299285351172367, (17) = 2.7657746379629105, (18) = 3.4960588686445333, (19) = 4.222709593140783, (20) = 4.9471575201173215, (21) = 5.670371337996262, (22) = 6.392927283827333, (23) = 7.1151455505991725, (24) = 7.8371254600843265, (25) = 8.558285432902858, (26) = 9.279151313463982, (27) = 10.0}, datatype = float[8], order = C_order); Y := Matrix(27, 3, {(1, 1) = -0.8430463960375806e-10, (1, 2) = .0, (1, 3) = 0.6372865716221984e-11, (2, 1) = -0.7278369204253125e-10, (2, 2) = 0.1597206925078892e-10, (2, 3) = 0.24052089989356827e-10, (3, 1) = -0.53829300642653995e-10, (3, 2) = 0.3900449379747658e-10, (3, 3) = 0.4808201125108471e-10, (4, 1) = -0.4953774824535353e-10, (4, 2) = 0.45113087361062455e-10, (4, 3) = 0.5272889316948668e-10, (5, 1) = -0.12035755690995432e-9, (5, 2) = -0.27735040596012685e-10, (5, 3) = -0.2587682090031095e-10, (6, 1) = -0.3348786263069427e-9, (6, 2) = -0.2451273076111555e-9, (6, 3) = -0.25162616591548056e-9, (7, 1) = -0.5269956859085185e-9, (7, 2) = -0.42674669378525426e-9, (7, 3) = -0.4325686117487571e-9, (8, 1) = 0.20412785836243726e-9, (8, 2) = 0.296577932431626e-9, (8, 3) = 0.28017398315903054e-9, (9, 1) = 0.24926163090190593e-8, (9, 2) = 0.2127611241294755e-8, (9, 3) = 0.17897368469198382e-8, (10, 1) = 0.5632465773934123e-9, (10, 2) = -0.9260334379037399e-9, (10, 3) = -0.15545759151451124e-8, (11, 1) = -0.9881440672433423e-8, (11, 2) = -0.6775803922126295e-8, (11, 3) = -0.30684471502496324e-8, (12, 1) = 0.23811084513591336e-8, (12, 2) = 0.8938239872691941e-8, (12, 3) = 0.6981346856881195e-8, (13, 1) = 0.7361322206069884e-8, (13, 2) = -0.16006959885630092e-8, (13, 3) = -0.60954833295139e-8, (14, 1) = -0.8135577900607999e-9, (14, 2) = -0.5635256156486533e-8, (14, 3) = 0.6633628564464468e-8, (15, 1) = -0.7801726603638605e-8, (15, 2) = 0.8656969516773866e-8, (15, 3) = -0.8563833248490694e-8, (16, 1) = 0.5807906930124011e-8, (16, 2) = -0.4097918191400617e-8, (16, 3) = -0.33689146247052684e-9, (17, 1) = 0.45105985824776024e-9, (17, 2) = -0.5007083474759643e-8, (17, 3) = 0.18471029972130062e-7, (18, 1) = -0.29145669603870902e-8, (18, 2) = 0.8373496931255616e-8, (18, 3) = -0.26162309001022457e-7, (19, 1) = 0.1715569635685495e-8, (19, 2) = -0.535670525811051e-8, (19, 3) = 0.19030608416671472e-7, (20, 1) = -0.34977609660274384e-9, (20, 2) = 0.1474813130678135e-8, (20, 3) = -0.65215600082828415e-8, (21, 1) = 0.20443269195515585e-10, (21, 2) = -0.13984210716486318e-9, (21, 3) = 0.8613227459894349e-9, (22, 1) = -0.17789355570891232e-11, (22, 2) = 0.4736297077901396e-11, (22, 3) = -0.18903663780540682e-10, (23, 1) = 0.34855282413257125e-12, (23, 2) = -0.4570031050682057e-11, (23, 3) = 0.3051816454556092e-10, (24, 1) = -0.5988722303653391e-12, (24, 2) = 0.1853672208396266e-11, (24, 3) = -0.12062893365927572e-10, (25, 1) = -0.12720293841230909e-12, (25, 2) = -0.5532970050855484e-12, (25, 3) = 0.5393029277707399e-11, (26, 1) = -0.16725993567675494e-12, (26, 2) = 0.4580253496963815e-12, (26, 3) = -0.25878083851462074e-11, (27, 1) = .0, (27, 2) = .0, (27, 3) = 0.13176846361390094e-11}, 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[27] elif outpoint = "order" then return 8 elif outpoint = "error" then return HFloat(2.6162309001022457e-8) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [3, 27, [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[27] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[27] 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`(27, 3, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(3, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(27, 3, X, Y, outpoint, yout, L, V) end if; [eta = outpoint, seq('[f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), 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[27] elif outpoint = "order" then return 8 elif outpoint = "error" then return HFloat(2.6162309001022457e-8) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [3, 27, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[27] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[27] 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(3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(27, 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 ) = (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(3, {(1) = 0., (2) = 0., (3) = 0.}); `dsolve/numeric/hermite`(27, 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, 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), diff(diff(f(eta), eta), eta)]'[i] = res[i+1], i = 1 .. 3)] catch: error end try end proc

> 

plots:-odeplot(dsol);

If you change a = 10 to a = 20, you will see that we get pretty much the
same curve, so for all practical purposes, a = 10 is good enough approximation

for infinity.

 

Here are the solutions corresponding to the choices a=10 and a=20 shown together.

Download mw.mw

Solution...

Rouben Rostamian 7711
June 27 2019
1 2

The sign of the second derivative term is certainly in error.  It should be negative.  A typo?

The easiest way to handle the integral is to introduce a new unknown, let's say v(x,t) defined by diff(v(x,t),t) = u(x,t) and v(x,0)=0.  Then the integral of u is v.

Call Maple's pdsolve() to solve two PDEs in the two unknowns u and v.  This is not done through Crank-Nicolson, but it gets the answer.  If you insist on Crank-Nicolson, then you will have to write the finite difference formulation yourself, although there is not much point to it other than for a programming exercise.  In any case, the introduction of the unknown v gets rid of the integral so you won't have to worry about integration at all.

Here is my worksheet.

> 

restart;
> 

pde1 := diff(u(x,t),t) - diff(u(x,t),x,x) + u(x,t) + v(x,t) = f(x,t);

diff(u(x, t), t)-(diff(diff(u(x, t), x), x))+u(x, t)+v(x, t) = f(x, t)

> 

pde2 := diff(v(x,t),t) = u(x,t);

diff(v(x, t), t) = u(x, t)

> 

f := (x,t) -> (x^3 + t^2*x^2 - t^2*x + 4*t*x - 2*t + 1)*exp(-x*t) - x + 1;

proc (x, t) options operator, arrow; (x^3+t^2*x^2-t^2*x+4*t*x-2*t+1)*exp(-t*x)-x+1 end proc

> 

bc := u(0,t) = 0, u(1,t) = 0;

u(0, t) = 0, u(1, t) = 0

> 

ic := u(x,0) = x*(1-x),  v(x,0) = 0;

u(x, 0) = x*(1-x), v(x, 0) = 0

> 

pdsol := pdsolve({pde1,pde2}, {ic, bc}, numeric);

_m140212070567648

> 

pdsol:-plot3d(u(x,t), x=0..1, t=0..10);


Download mw.mw

Get a legal copy...

Rouben Rostamian 7711
June 27 2019
1 0

Where did you find your Maple?

Alternative...

Rouben Rostamian 7711
June 27 2019
0 3

Here's a slightly different approach.

> 

restart;
> 

interface(rtablesize=12):
> 

M := < seq('N[i],0', i=1..6);  
       seq('0,N[i]', i=1..6) >;

Matrix(2, 12, {(1, 1) = N[1], (1, 2) = 0, (1, 3) = N[2], (1, 4) = 0, (1, 5) = N[3], (1, 6) = 0, (1, 7) = N[4], (1, 8) = 0, (1, 9) = N[5], (1, 10) = 0, (1, 11) = N[6], (1, 12) = 0, (2, 1) = 0, (2, 2) = N[1], (2, 3) = 0, (2, 4) = N[2], (2, 5) = 0, (2, 6) = N[3], (2, 7) = 0, (2, 8) = N[4], (2, 9) = 0, (2, 10) = N[5], (2, 11) = 0, (2, 12) = N[6]})

> 

Z := < D[1], 0;  0, D[2];  D[1], D[2] >;

Matrix(3, 2, {(1, 1) = D[1], (1, 2) = 0, (2, 1) = 0, (2, 2) = D[2], (3, 1) = D[1], (3, 2) = D[2]})

> 

Matrix(3,12, (i,j)-> add(Z[i,k](M[k,j])(x,y),k=1..2)):
convert~(%,diff);

Matrix(3, 12, {(1, 1) = diff(N[1](x, y), x), (1, 2) = 0, (1, 3) = diff(N[2](x, y), x), (1, 4) = 0, (1, 5) = diff(N[3](x, y), x), (1, 6) = 0, (1, 7) = diff(N[4](x, y), x), (1, 8) = 0, (1, 9) = diff(N[5](x, y), x), (1, 10) = 0, (1, 11) = diff(N[6](x, y), x), (1, 12) = 0, (2, 1) = 0, (2, 2) = diff(N[1](x, y), y), (2, 3) = 0, (2, 4) = diff(N[2](x, y), y), (2, 5) = 0, (2, 6) = diff(N[3](x, y), y), (2, 7) = 0, (2, 8) = diff(N[4](x, y), y), (2, 9) = 0, (2, 10) = diff(N[5](x, y), y), (2, 11) = 0, (2, 12) = diff(N[6](x, y), y), (3, 1) = diff(N[1](x, y), x), (3, 2) = diff(N[1](x, y), y), (3, 3) = diff(N[2](x, y), x), (3, 4) = diff(N[2](x, y), y), (3, 5) = diff(N[3](x, y), x), (3, 6) = diff(N[3](x, y), y), (3, 7) = diff(N[4](x, y), x), (3, 8) = diff(N[4](x, y), y), (3, 9) = diff(N[5](x, y), x), (3, 10) = diff(N[5](x, y), y), (3, 11) = diff(N[6](x, y), x), (3, 12) = diff(N[6](x, y), y)})


 

Download mw.mw

Solution and perhaps a Maple bug...

Rouben Rostamian 7711
June 25 2019
1 10

In your worksheet you have defined many symbols as functions, some with improper syntax.  But there is no need to define a symbol as a function if defining it as an expression will do.  I have removed the unnecessary function definitions from your worksheet, and also have changed the three occurrences of pi to Pi, and 1.6 to 16/10.

After those fixes the worksheet executes but fails to find the solution of the trivially simple differential equation at the end.  That looks like a Maple bug to me.

However, that differential equation may be solved trivially direct integration.  I have done that at the end of the worksheet.

> 

restart;
> 

n := 3;

3

> 

my := R*(16/10+(1/2)*sin((1/10)*Pi*t))/`&lambda;opt`;

R*(8/5+(1/2)*sin((1/10)*Pi*t))/`&lambda;opt`

> 

p1op := 16/10+(1/2)*sin((1/10)*Pi*t);

8/5+(1/2)*sin((1/10)*Pi*t)

> 

f6 := diff(v1(t),t) + (1/2)*(p1-p1op)^2 = 0;

diff(v1(t), t)+(1/2)*(p1-8/5-(1/2)*sin((1/10)*Pi*t))^2 = 0

> 

dsolve({f6, v1(20) = 0}, v1(t)):
t2 := rhs(%);

-(537/400)*t+(8/5)*p1*t-(1/2)*p1^2*t-5*p1*cos((1/10)*Pi*t)/Pi+(5/16)*sin((1/5)*Pi*t)/Pi+8*cos((1/10)*Pi*t)/Pi+537/20-32*p1+10*p1^2+5*p1/Pi-8/Pi

> 

v[2] := t2;

-(537/400)*t+(8/5)*p1*t-(1/2)*p1^2*t-5*p1*cos((1/10)*Pi*t)/Pi+(5/16)*sin((1/5)*Pi*t)/Pi+8*cos((1/10)*Pi*t)/Pi+537/20-32*p1+10*p1^2+5*p1/Pi-8/Pi

> 

m1 := (1/2)*b11*ro*Pi*R^2*(my*cp2/(R*p1)-cp3)/p1;

(1/2)*b11*ro*Pi*R^2*((8/5+(1/2)*sin((1/10)*Pi*t))*cp2/(`&lambda;opt`*p1)-cp3)/p1

> 

d1 := diff(v[3-1], p1);

(8/5)*t-p1*t-5*cos((1/10)*Pi*t)/Pi-32+20*p1+5/Pi

> 

g := diff(v[3](t), t) + d1*(a11*p1+a13*p3+k11*my^2+m1) = 0;

diff(v[3](t), t)+((8/5)*t-p1*t-5*cos((1/10)*Pi*t)/Pi-32+20*p1+5/Pi)*(k11*R^2*(8/5+(1/2)*sin((1/10)*Pi*t))^2/`&lambda;opt`^2+a11*p1+a13*p3+(1/2)*b11*ro*Pi*R^2*((8/5+(1/2)*sin((1/10)*Pi*t))*cp2/(`&lambda;opt`*p1)-cp3)/p1) = 0

OK, we are ready to solve the differential equation but we see that Maple

fails to evaluate the integral.  This looks like a Maple bug to me.

> 

dsolve({g, v[3](20)=0}, v[3](t));

v[3](t) = Int(-(1/500)*(5*p1*_z1*Pi-100*p1*Pi-8*Pi*_z1+25*cos((1/10)*Pi*_z1)+160*Pi-25)*(50*Pi*R^2*b11*cp3*p1*ro*`&lambda;opt`^2-25*Pi*sin((1/10)*Pi*_z1)*R^2*b11*cp2*ro*`&lambda;opt`+25*R^2*cos((1/10)*Pi*_z1)^2*k11*p1^2-80*Pi*R^2*b11*cp2*ro*`&lambda;opt`-160*sin((1/10)*Pi*_z1)*R^2*k11*p1^2-100*a11*p1^3*`&lambda;opt`^2-100*a13*p3*`&lambda;opt`^2*p1^2-281*k11*R^2*p1^2)/(Pi*`&lambda;opt`^2*p1^2), _z1 = 20 .. t)

In fact, the differential equation is quite trivial since it is of the form
dv/dt = f(t) so solving it is a matter of integrating f(t)which is a
simple polynomial in t, sin(a*t), cos(a*t).  Here I will solve it "by hand":

> 

isolate(g, diff(v[3](t), t)):
int(rhs(%), t):
dsol := v[3](t) = simplify(% - eval(%,t=20));

v[3](t) = (1/6000)*(-37500*R^2*(cp2*`&lambda;opt`*Pi*ro*b11+(2/3)*k11*sin((1/10)*Pi*t)*p1^2+p1^3*k11+(24/5)*k11*p1^2)*cos((1/10)*Pi*t)^2-15000*(cp2*`&lambda;opt`*Pi*ro*b11+(1/2)*k11*sin((1/10)*Pi*t)*p1^2+(32/5)*k11*p1^2)*(Pi*(t-20)*p1-5+(-(8/5)*t+32)*Pi)*R^2*cos((1/10)*Pi*t)+150000*p1*(((32/5)*R^2*k11+2*`&lambda;opt`^2*a11)*p1^2+(2*`&lambda;opt`^2*p3*a13-(743/150)*R^2*k11)*p1+R^2*`&lambda;opt`*Pi*ro*b11*(-cp3*`&lambda;opt`+cp2))*sin((1/10)*Pi*t)+3000*`&lambda;opt`^2*Pi^2*a11*(t-20)^2*p1^4+(8055*(t-20)^2*(R^2*k11+(200/537)*(a13*p3-(8/5)*a11)*`&lambda;opt`^2)*Pi^2-30000*`&lambda;opt`^2*a11*(t-20)*Pi+37500*R^2*k11)*p1^3+(-1500*R^2*cp3*`&lambda;opt`^2*ro*b11*(t-20)^2*Pi^3-12888*((200/537)*`&lambda;opt`^2*p3*a13+R^2*k11)*(t-20)^2*Pi^2-80550*((200/537)*`&lambda;opt`^2*p3*a13+R^2*k11)*(t-20)*Pi-300000*R^2*k11)*p1^2+2400*b11*`&lambda;opt`*(t-20)*ro*((t-20)*(cp3*`&lambda;opt`+cp2)*Pi+(25/4)*cp3*`&lambda;opt`)*R^2*Pi^2*p1-3840*b11*`&lambda;opt`*ro*R^2*cp2*Pi*(25/8+Pi*(t-20))^2)/(Pi^2*p1^2*`&lambda;opt`^2)

Here we verify that we have obtained the correct solution:

> 

odetest(dsol, [g, v[3](20)=0]);

[0, 0]


 

Download mw.mw

Try this...

Rouben Rostamian 7711
June 23 2019
1 1 
restart;
eq := (x-2)^2*(x-3)+epsilon =0;
solve(eq, epsilon);
plot([%, x, x=1..3.5], view=[-1..1, 0..4], labels=[epsilon,x]);

Solution...

Rouben Rostamian 7711
June 21 2019
1 4

See if this is useful to you.

> 

restart;

You have x as a 17*1 matrix.  It will simplify things if you
present x as a list of 17 numbers.  Here I will do that conversion:

> 

<<0.556960605000000>,<3.39039994000000>,<2.09005200300000>,<0.645104568000000>,<5.31340491600000>,<3.32743462200000>,<0.635452131000000>,<1.56878297000000>,<0.282764039000000>,<1.02862059900000>,<3.14927606700000>,<0.654644768000000>,<1.30502450500000>,<2.13887537900000>,<2.11803658900000>,<7.29488570500000>,<0.478693554000000>>:
my_x := convert(%, list);

[.556960605000000, 3.39039994000000, 2.09005200300000, .645104568000000, 5.31340491600000, 3.32743462200000, .635452131000000, 1.56878297000000, .282764039000000, 1.02862059900000, 3.14927606700000, .654644768000000, 1.30502450500000, 2.13887537900000, 2.11803658900000, 7.29488570500000, .478693554000000]

Find the length of the list:

> 

n := nops(my_x);

17

Calculate xbar:

> 

my_xbar := add(my_x)/n;

2.116377234

Several comments regarding b:

• 

Since n is known, it is better to use add rather than sum in
your expression for b.
• 

Additionally, the exponent 1.5 is better presented as 3/2.
• 

Most important of all, the square brackets in the denominator
should be made into regular parentheses.  Square brackets
mean something quite different in Maple.
• 

For convenience, we define b as a function of x and xbar.

Putting all these together, we get

> 

(add((x[i]-xbar)^3, i = 1 .. n))/(n*((add((x[i]-xbar)^2, i = 1 .. n))/(n-1))^(3/2)):
b := unapply(%, [x,xbar]);

proc (x, xbar) options operator, arrow; (1/17)*((x[1]-xbar)^3+(x[2]-xbar)^3+(x[3]-xbar)^3+(x[4]-xbar)^3+(x[5]-xbar)^3+(x[6]-xbar)^3+(x[7]-xbar)^3+(x[8]-xbar)^3+(x[9]-xbar)^3+(x[10]-xbar)^3+(x[11]-xbar)^3+(x[12]-xbar)^3+(x[13]-xbar)^3+(x[14]-xbar)^3+(x[15]-xbar)^3+(x[16]-xbar)^3+(x[17]-xbar)^3)/((1/16)*(x[1]-xbar)^2+(1/16)*(x[2]-xbar)^2+(1/16)*(x[3]-xbar)^2+(1/16)*(x[4]-xbar)^2+(1/16)*(x[5]-xbar)^2+(1/16)*(x[6]-xbar)^2+(1/16)*(x[7]-xbar)^2+(1/16)*(x[8]-xbar)^2+(1/16)*(x[9]-xbar)^2+(1/16)*(x[10]-xbar)^2+(1/16)*(x[11]-xbar)^2+(1/16)*(x[12]-xbar)^2+(1/16)*(x[13]-xbar)^2+(1/16)*(x[14]-xbar)^2+(1/16)*(x[15]-xbar)^2+(1/16)*(x[16]-xbar)^2+(1/16)*(x[17]-xbar)^2)^(3/2) end proc

We are ready to compute c__1 which we view as a function of the
independent variables x, and xbar:

> 

diff(b(x,xbar), xbar):
c_1 := unapply(%, [x,xbar]):

As to c__2, it is a function of x, xbar, and i.  Note that here I am taking xbar and x
as independent variables as you have requested, therefore xbar does not get
differentiated even though behind the scenes it depends on x[i].

> 

c_2 := proc(x,xbar,i)
  local X;
  diff(b(X,xbar), X[i]);
  eval(%, X=x);
end proc:

OK, we are ready to try things out:

> 

c_1(my_x, my_xbar);

-1.479074045

> 

c_2(my_x, my_xbar, 1);

.1638601668

> 

c_2(my_x, my_xbar, 9);

.2054258169

> 

c_2(my_x, my_xbar, 17);

.1753357291

 

 


 

Download mw.mw

Solution...

Rouben Rostamian 7711
June 21 2019
2 2

The exponential function is written exp(...) in Maple.  I have fixed that in your worksheet.  Here is what it takes to evaluate your integral.

> 

restart;
> 

kernelopts(version);

`Maple 2019.1, X86 64 LINUX, May 21 2019, Build ID 1399874`

Your integrand is:

> 

P := exp(-h*Pi*v*(x+y)/beta) / cosh(Pi*v*(y-x)/beta)^(3+h);

exp(-h*Pi*v*(x+y)/beta)/cosh(Pi*v*(y-x)/beta)^(3+h)

and you wish to evaluate

> 

Int(P, y=-infinity..z, x=-infinity..w);

Int(exp(-h*Pi*v*(x+y)/beta)/cosh(Pi*v*(y-x)/beta)^(3+h), y = -infinity .. z, x = -infinity .. w)

I believe that this integral will converge provided that h < 0 (added in edit: and h>−3?)
Furthermore, I believe that to obtain an answer in terms of
elementary functions, h needs to be a negative integer.
Let's take h = -1

> 

Q := Int(eval(P, h=-1), y=-infinity..z, x=-infinity..w);

Int(exp(Pi*v*(x+y)/beta)/cosh(Pi*v*(y-x)/beta)^2, y = -infinity .. z, x = -infinity .. w)

> 

value(Q) assuming beta>0, v>0, w::real, z::real;

-(1/2)*beta^2*(2*exp(2*w*Pi*v/beta)*arctan(exp((-z+w)*Pi*v/beta))+2*arctan(exp(-(-z+w)*Pi*v/beta))*exp(2*Pi*v*z/beta)-exp(2*w*Pi*v/beta)*Pi-exp(2*Pi*v*z/beta)*Pi+2*exp((w+z)*Pi*v/beta))/(Pi^2*v^2)


 

Download mw.mw

Two issues...

Rouben Rostamian 7711
June 21 2019
1 0

In the first pass through the for-loop we have R=0, in which case EQ2 reduces to k=0.  Consequently EQ1 reduces to 0=0 and therefore wr cannot be determined.  That's what leads to the error message that you are seeing.

This brings us to the following question.  Regardless of the value of R, your system has k=0, wr=0 as a solution as you can see by doing eval(eqns, {k=0, wr=0}).  You don't need a for-loop (or Maple) to see that.  So what is it that you are really trying to do?

 

 

 

This may do what you want...

Rouben Rostamian 7711
June 17 2019
2 2
> 

restart;
> 

with(plots):

Vertices

> 

v[1], v[2], v[3] := <cos(Pi/6), sin(Pi/6),sqrt(2)/2>,
                    <-cos(Pi/6), sin(Pi/6),sqrt(2)/2>,
                    <0,-1,sqrt(2)/2>;
v[4], v[5], v[6] := <-cos(Pi/6), -sin(Pi/6),-sqrt(2)/2>,
                    <cos(Pi/6), -sin(Pi/6),-sqrt(2)/2>,
                    <0,1,-sqrt(2)/2>;

Let's seee what it looks like:

> 

display([
  pointplot3d([v[1],v[2],v[3],v[1]], connect),
  pointplot3d([v[4],v[5],v[6],v[4]], connect),
  pointplot3d([v[1],v[5],v[3],v[4],v[2],v[6],v[1]], connect)
], scaling=constrained);

The parameter t > 0 stretches the base.

> 

frame := proc(t)
  local V, i;
  V[1], V[2], V[3] := v[1], v[2], v[3];
  V[4], V[5], V[6] := (1+t)*v[4], (1+t)*v[5], (1+t)*v[6];
  display([
    pointplot3d([V[1],V[2],V[3],V[1]], connect),
    pointplot3d([V[4],V[5],V[6],V[4]], connect),
    pointplot3d([V[1],V[5],V[3],V[4],V[2],V[6],V[1]], connect)
  ], scaling=constrained);
end proc:
> 

animate(frame, [t], t=0..3, frames=50, paraminfo=false,
  scaling=constrained, orientation=[-90,0,0], axes=none);

Download stretch-octahedron.mw

Does this help?...

Rouben Rostamian 7711
June 16 2019
1 1

Add this at the end of your worksheet:

u__lower := (x,t) -> (-t - sqrt(t^2 - 4*x))/4;
lower_red := plot3d(u__lower(x,t), x=-140..40,t=-20..20);
x_t_curve := odeplot(p, color=red, thickness=3);
x_t_curve_lifted := plottools:-transform((t,x)->[x,t,u__lower(x,t)])(x_t_curve);
display([lower_red, x_t_curve_lifted]);
First 23 24 25 26 27 28 29 Last Page 25 of 53