Advanced Search

Rouben Rostamian

8586 Reputation

23 Badges

19 years, 140 days
Contact Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

Solution and perhaps a Maple bug...

Rouben Rostamian 8586
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))/`λopt`;

R*(8/5+(1/2)*sin((1/10)*Pi*t))/`λ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/(`λ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/`λopt`^2+a11*p1+a13*p3+(1/2)*b11*ro*Pi*R^2*((8/5+(1/2)*sin((1/10)*Pi*t))*cp2/(`λ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*`λopt`^2-25*Pi*sin((1/10)*Pi*_z1)*R^2*b11*cp2*ro*`λopt`+25*R^2*cos((1/10)*Pi*_z1)^2*k11*p1^2-80*Pi*R^2*b11*cp2*ro*`λopt`-160*sin((1/10)*Pi*_z1)*R^2*k11*p1^2-100*a11*p1^3*`λopt`^2-100*a13*p3*`λopt`^2*p1^2-281*k11*R^2*p1^2)/(Pi*`λ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*`λ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*`λ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*`λopt`^2*a11)*p1^2+(2*`λopt`^2*p3*a13-(743/150)*R^2*k11)*p1+R^2*`λopt`*Pi*ro*b11*(-cp3*`λopt`+cp2))*sin((1/10)*Pi*t)+3000*`λopt`^2*Pi^2*a11*(t-20)^2*p1^4+(8055*(t-20)^2*(R^2*k11+(200/537)*(a13*p3-(8/5)*a11)*`λopt`^2)*Pi^2-30000*`λopt`^2*a11*(t-20)*Pi+37500*R^2*k11)*p1^3+(-1500*R^2*cp3*`λopt`^2*ro*b11*(t-20)^2*Pi^3-12888*((200/537)*`λopt`^2*p3*a13+R^2*k11)*(t-20)^2*Pi^2-80550*((200/537)*`λopt`^2*p3*a13+R^2*k11)*(t-20)*Pi-300000*R^2*k11)*p1^2+2400*b11*`λopt`*(t-20)*ro*((t-20)*(cp3*`λopt`+cp2)*Pi+(25/4)*cp3*`λopt`)*R^2*Pi^2*p1-3840*b11*`λopt`*ro*R^2*cp2*Pi*(25/8+Pi*(t-20))^2)/(Pi^2*p1^2*`λ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 8586
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 8586
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 8586
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 8586
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 8586
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 8586
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]);

Take y to be a function of x.  Then...

Rouben Rostamian 8586
June 15 2019
5 0

Take y to be a function of x. Then differentiate the equation once to get rid of B. Then isolate A, and then differentiate again to get rid of A. The resulting differential equation will not have A or B in it. Here is how.

> 

restart;
> 

t1 := x^2 + A*y(x)^2 = B;

x^2+A*y(x)^2 = B

> 

t2 := diff(t1, x);

2*x+2*A*y(x)*(diff(y(x), x)) = 0

> 

t3 := isolate(t2, A);

A = -x/(y(x)*(diff(y(x), x)))

> 

t4 := diff(t3, x);

0 = -1/(y(x)*(diff(y(x), x)))+x/y(x)^2+x*(diff(diff(y(x), x), x))/(y(x)*(diff(y(x), x))^2)

That's the differential equation you are looking for, but let's simplify it:

> 

de := factor(isolate(t4, diff(y(x),x,x)));

diff(diff(y(x), x), x) = -(x*(diff(y(x), x))-y(x))*(diff(y(x), x))/(y(x)*x)

OK.  Let's see if Maple can recover the original family of curves from this:

> 

dsolve(de, y(x), implicit);

_C1*x^2-_C2+(1/2)*y(x)^2 = 0

Yes!  That's equivalent to the original equation, with some re-named constants.


 

Download mw.mw

PS: The English verb for computing a derivative is "to differentiate", not "to derive", but don't feel bad for having mixed up the two.  Some native English speakers make the same mistake.

 

Geodesics...

Rouben Rostamian 8586
June 12 2019
4 3

Some time ago when you asked about a related question, I made worksheets for calculating geodesics on all five regular polyhedra but didn't post them because I felt that there ought to be a better approach than what I have right now.  Now that you are asking again, let me show how I do geodesics on a cube.  I have attempted to document the approach in cube.mw in detail, but let me know if you have questions.  I have not documented the work on the other four regular polyhedra so I am not posting those worksheets.  I do, however, post a sample of a geodesic on a dodecahedron just to show that I have done it.

 

Download worksheet: cube.mw

 

Animation of the Geneva Drive...

Rouben Rostamian 8586
June 07 2019
3 0

I saw this question a few days ago and was unable to read the original undecipherable code, so I wrote my own code from scratch which I did not post then since I was waiting for an update from the original poster.

Checking the thread today I see that there has been quite a bit of discussion in the meantime, so perhaps this is a good time to post my code which is entirely independent of that presented in the discussion.

Here are a few sample animations for your enjoyment!  See Geneva Drive in Wikipedia for more info.

 

Download geneva.mw

 

Missing colon...

Rouben Rostamian 8586
June 04 2019
1 3

Insert missing colon in
if flag=1 then  dumbindex:=dumbindex+1  xdumb[dumbindex]:=x[i]: ...
to make it
if flag=1 then  dumbindex:=dumbindex+1:  xdumb[dumbindex]:=x[i]: ...

Solution with fixes...

Rouben Rostamian 8586
June 02 2019
2 6
> 

restart;
> 

with(plots):
> 

with(CurveFitting):
> 

pts := <0, 3.016315859, 5.753536424, 8.489713643, 11.22797756, 13.96519812, 16.09414745, 18.83136802, 21.26445296, 24.00167353, 26.73889409, 29.47611466, 31.60506399, 34.34228455, 37.07950512, 39.51259007, 41.6415394, 44.37875996, 46.81184491, 49.54906547, 51.6780148, 54.11109975, 56.84832032, 58.97726964, 62.32276145, 64.45171078, 66.88479572, 69.01374505, 71.44578665, 74.1814422, 76.91814109, 79.0455254, 81.47704533, 84.21165752, 86.33695514, 88.76691004, 90.89116431, 93.3195542, 96.04894966, 98.17059556, 100.5958554, 103.3231642, 105.4432451, 107.8679832, 110.5947703, 113.3236441, 116.0514746, 118.1746855, 120.9066893, 125.4640285, 127.8960701, 130.634334, 132.763805, 135.1994984, 137.9382839, 140.06932, 142.5034483, 145.2411905, 147.976846, 150.7119799, 153.4450271, 156.1770309, 158.9069481, 161.9399575, 164.6693529, 167.0951345, 169.2178237, 171.9472192, 174.3756091, 176.4998633, 179.2334322, 181.9680444, 184.4011293, 186.531122, 189.2719943, 192.0149532, 194.4542983, 197.2014306, 201.7822452, 203.921628, 207.5889923, 210.3397763, 213.699875, 217.3651526, 220.1133283, 222.859939, 224.9961918, 227.7412375, 230.1795391, 232.3131835, 234.7504419, 237.4913141, 240.5368437, 242.9720153, 245.7108009, 247.8313153, 247.8413153, 247.8513153, 250>,
<-0.041150436, -0.041150436, -0.04113729, -0.038722776, -0.041110998, -0.041097851, -0.041087627, -0.04107448, -0.041062795, -0.041049649, -0.041036503, -0.041023356, -0.041013132, -0.040999985, -0.040986839, -0.040975154, -0.040964929, -0.040951783, -0.040940097, -0.040926951, -0.040916726, -0.040905041, -0.040891895, -0.04088167, -0.040865602, -0.040855378, -0.040843692, -0.040833467, -0.038420415, -0.034805218, -0.033591388, -0.029979112, -0.026365376, -0.020348812, -0.011933802, -0.004718015, 0.006098363, 0.0169162, 0.034939601, 0.051759396, 0.069781335, 0.09260747, 0.113029316, 0.132251939, 0.156278757, 0.175502841, 0.197128292, 0.210346036, 0.222366019, 0.233194081, 0.235607134, 0.233218913, 0.232028454, 0.226036722, 0.222447817, 0.217655307, 0.215265626, 0.214078088, 0.217693285, 0.222509166, 0.232127781, 0.244147763, 0.26097048, 0.280196024, 0.298219424, 0.31504068, 0.329459108, 0.347482508, 0.358300346, 0.369116723, 0.377534655, 0.383551219, 0.383562904, 0.381171762, 0.372780123, 0.35958575, 0.345189232, 0.32238939, 0.27918669, 0.255183243, 0.214377529, 0.183172901, 0.149569828, 0.113566848, 0.088365639, 0.06676648, 0.049967135, 0.031970027, 0.019974876, 0.009178949, -0.000414835, -0.008806474, -0.018397336, -0.023188385, -0.026777289, -0.030369115, -0.030369115, -0.030369115, -0.030369115>;

"pts:=[[[[[0],[3.016315859],[5.753536424],[8.489713643],[11.22797756],[13.96519812],[16.09414745],[18.83136802],[21.26445296],[24.00167353],[&vellip;]]]],["99 element Vector[column]"]],[[[[[-0.041150436],[-0.041150436],[-0.04113729],[-0.038722776],[-0.041110998],[-0.041097851],[-0.041087627],[-0.04107448],[-0.041062795],[-0.041049649],[&vellip;]]]],["99 element Vector[column]"]]"

> 

A := unapply(Spline(pts[1], pts[2], x, degree = 3), x):
> 

plot(A(x), x=0..250);

> 

PDE := diff(C(x,t),t) = diff(C(x,t),x,x) - diff(C(x,t),x) + 5*diff(C(x,t), x$4);

diff(C(x, t), t) = diff(diff(C(x, t), x), x)-(diff(C(x, t), x))+5*(diff(diff(diff(diff(C(x, t), x), x), x), x))

> 

IBC := C(x,0)=A(x), D[1](C)(0,t)=0, D[1](C)(250,t)=0, D[1,1](C)(0,t)=0, D[1,1](C)(250,t)=0:
> 

pdsol := pdsolve(PDE, {IBC}, numeric);

_m140347094656320

> 

pdsol:-plot3d(x=0..250, t=0..600, orientation=[120,65,0]);

In the long run the solution converges to a constant.
Here is what it looks like at t=600.

> 

pdsol:-plot(t=600, view=-0.05..0, gridlines=false);


 

Download mw.mw

Another way...

Rouben Rostamian 8586
May 30 2019
1 0

As an alternative to what tomleslie has suggested, you may also want to try:

restart;
with(plots):
P, Q := <1,1>, <2,2>:
frame := t -> pointplot([P, (1-t)*P + t*Q], connect, color="Green"):
nf := 50:
display([seq(frame(i/(nf-1)), i=0..nf-1)], insequence);

 

Solution via finite differences...

Rouben Rostamian 8586
May 24 2019
0 2

As far as I can tell, Maple's numerical solver is not equipped to solve a heat conduction problem with discontinuous conductivity (but I will be happy to be corrected on this.)  For now, I have written a finite differences solver with a specific goal of handling such cases which you will find in
https://www.mapleprimes.com/posts/210623-A-Finite-Difference-Scheme-For-The-Heat

Get the worksheet from there and then solve your problem by entering:
k := x -> piecewise(x < 10, 10, 20);
sol := heat_solve(c=k, a=0, b=20, alpha=20, beta=10, u__0=0, T=15, n=50, m=50):

Then animate the solution:
display(convert(sol, list), insequence);
 

How about this?...

Rouben Rostamian 8586
May 20 2019
2 3
> 

restart;

The equation of heat condition in a medium of
variable conductivity k(x) is
"(&PartialD; u)/(&PartialD; t)=(&PartialD;)/(&PartialD; x)(k(x) (&PartialD; u)/(&PartialD; x))."

In your case the conductivity is a piecewise constant
function which jumps from k__1 to k__2 at the interval's

midpoint. In principle we should be able to take

"k(x) = `k__1`+(`k__2`-`k__1`) Heaviside(x-L),"

or equivalently

"k(x)=`k__1`+(`k__2`-`k__1`) piecewise(x<L, `k__1`, `k__2`)."

I am unable to make either of these representations
work with Maple's numeric solver, so I do the next
best thing which is to approximate the Heaviside function
with a smooth function such as

> 

myHeaviside := x -> (1+erf(4*x))/2;

proc (x) options operator, arrow; 1/2+(1/2)*erf(4*x) end proc

Let's compare myHeaviside with the real thing:

> 

plot([Heaviside(x), myHeaviside(x)], x=-10..10,
  color=[red,blue], axes=frame, gridlines=false,
  legend=[Heavisie, myHeaviside]);

OK, then, let's start:

> 

pde := Diff(u(x,t),t) = Diff(k(x)*Diff(u(x,t),x), x);

Diff(u(x, t), t) = Diff(k(x)*(Diff(u(x, t), x)), x)

> 

ibc := u(x,0)=0, u(0,t)=v1,  u(2*L,t)=v2;

u(x, 0) = 0, u(0, t) = v1, u(2*L, t) = v2

> 

k := x -> k1 + (k2-k1)*myHeaviside(x-L);

proc (x) options operator, arrow; k1+(k2-k1)*myHeaviside(x-L) end proc

> 

L := 10: v1 := 20: v2 := 10: k1 := 10: k2 := 20:

To handle the sharp change in conductivity, we refine the spacestep size from

its default value of 2*L/20 by dividing it by 5.

> 

pdsol := pdsolve(pde, {ibc}, numeric, spacestep=(2*L/20)/5);

_m139639715851360

> 

pdsol:-animate(t=15, frames=50, labels=[x,u(x,t)],
  title="time = %f");

Download mw.mw

PS:  Looking at the animation closely, we see that the graph wobbles a bit near the initial time. That's an artifact due to our insufficiently small step size. Reducing the spacestep from (2*L)/20/5 to (2*L)/20/20 makes the wobble go away.

First 28 29 30 31 32 33 34 Last Page 30 of 58