Adri van der Meer

Adri vanderMeer

1345 Reputation

16 Badges

14 years, 32 days
University of Twente (retired)
Enschede, Netherlands

My "website" consists of a Maple Manual in Dutch

MaplePrimes Activity


These are replies submitted by Adri van der Meer

@Bendesarts Use output=listprocedure, and plot as a parametric curve:

sol:=dsolve({sys,Cinit},numeric,output=listprocedure):
Sol := subs( sol, [x(t),z(t)] ): X,Z := Sol[]:
plot( [X(t),Z(t),t=0..10],numpoints=200,color=blue,legend="z(x)");

@sunit Why Newton-Raphson and not just fsolve? What did you try? Submit a new question if your attempt doesn't work properly.

@sunit Why Newton-Raphson and not just fsolve? What did you try? Submit a new question if your attempt doesn't work properly.

@sunit Notice that I start with omega>0 to prevent division by zero:

pnts := [seq( [omega, fsolve(f)], omega=0.001..10,0.1 )]:
y := (omega,x) -> [x,x^2/omega]:
ypts := (y@op)~(pnts):
plot(ypts);

Explanation: I make y a function of two variables, returning a list consisting of the second variable and the calculated y-value.
Now you can not apply y directly to all elements of the list pnts by y~(pnts), because these elements are lists, so we make a composition of y and op, which extract the elements of the sublists of pnts.

@sunit Notice that I start with omega>0 to prevent division by zero:

pnts := [seq( [omega, fsolve(f)], omega=0.001..10,0.1 )]:
y := (omega,x) -> [x,x^2/omega]:
ypts := (y@op)~(pnts):
plot(ypts);

Explanation: I make y a function of two variables, returning a list consisting of the second variable and the calculated y-value.
Now you can not apply y directly to all elements of the list pnts by y~(pnts), because these elements are lists, so we make a composition of y and op, which extract the elements of the sublists of pnts.

@Alejandro Jakubi and @Carl Love 

This is not the first time that you teach me something I didn't know! I was aware of the distinction between x^(1/2) and x^.5. A demonstration:

 

a := sqrt(4);

2

(1)

b := 4^(1/2);

4^(1/2)

(2)

c := 4^.5;

2.000000000

(3)

4.^(1/2);

2.000000000

(4)

if a=b then "a=b" else "a<>b" end if;

"a<>b"

(5)

if a=c then "a=c" else "a<>c" end if;

"a=c"

(6)

if b=c then "b=c" else "b<>c" end if;

"b<>c"

(7)

 

 

 

Download Sqrt.mw

 

@Alejandro Jakubi and @Carl Love 

This is not the first time that you teach me something I didn't know! I was aware of the distinction between x^(1/2) and x^.5. A demonstration:

 

a := sqrt(4);

2

(1)

b := 4^(1/2);

4^(1/2)

(2)

c := 4^.5;

2.000000000

(3)

4.^(1/2);

2.000000000

(4)

if a=b then "a=b" else "a<>b" end if;

"a<>b"

(5)

if a=c then "a=c" else "a<>c" end if;

"a=c"

(6)

if b=c then "b=c" else "b<>c" end if;

"b<>c"

(7)

 

 

 

Download Sqrt.mw

 

@N00bstyle I added the command

print(subs(lengte = _lengte, eigenfrequentie));

in your procedure reductiefactorvoetgangerslengte. Now I get

lengte := 22;
reductiefactorvoetgangerslengte(lengte);
                0.0077 √ 234256     
Error, (in reductiefactorvoetgangerslengte) cannot determine if this expression is true or false: 0.769201511143595e-2*234256^(1/2) <= 1.25

but

evalb(0.77e-2*sqrt(234256)<=1.25);
                             false

@N00bstyle I added the command

print(subs(lengte = _lengte, eigenfrequentie));

in your procedure reductiefactorvoetgangerslengte. Now I get

lengte := 22;
reductiefactorvoetgangerslengte(lengte);
                0.0077 √ 234256     
Error, (in reductiefactorvoetgangerslengte) cannot determine if this expression is true or false: 0.769201511143595e-2*234256^(1/2) <= 1.25

but

evalb(0.77e-2*sqrt(234256)<=1.25);
                             false

... because f(1) = a+b, so perfectly well defined. But the right derivative only exists if a+b=1.

... because f(1) = a+b, so perfectly well defined. But the right derivative only exists if a+b=1.

Of course. I only wanted to show an inexpected pecularity!

My answer focuses on the question why youre original code didn't work. This results in a rather obscure program. This can be made more transparent whwn you realize that eq1 and eq2 are intended to be dependent on the variable a, so make these functions of a:

restart; B := 53*Pi*(1/180): e := 3:
eq1 := a -> -2.005689708*a:
eq2 := a -> -2.005689708*a+5.369606170:
f := proc (a)
  if a <= 0 then 0
  elif a <= evalf(e*sin((1/2)*B)) then eq1(a)
  elif a <= evalf(2*e*sin((1/2)*B)) then eq2(a)
  else 0
  end if
end proc:

My answer focuses on the question why youre original code didn't work. This results in a rather obscure program. This can be made more transparent whwn you realize that eq1 and eq2 are intended to be dependent on the variable a, so make these functions of a:

restart; B := 53*Pi*(1/180): e := 3:
eq1 := a -> -2.005689708*a:
eq2 := a -> -2.005689708*a+5.369606170:
f := proc (a)
  if a <= 0 then 0
  elif a <= evalf(e*sin((1/2)*B)) then eq1(a)
  elif a <= evalf(2*e*sin((1/2)*B)) then eq2(a)
  else 0
  end if
end proc:

@Heeka 

(1) Y_vals produces a list of y(x), for x=0, 0.1, 0.2, ..., 1.0.

(2) Comparison exact and numeric solution:

restart;
DE := (x+b*y(x))*diff(y(x),x) + y(x) = 0;
BC := y(1)=1:
b := 1/10:
# exact solution:
sol1 := dsolve( {DE,BC}, y(x) ): Y1 := unapply( subs( sol1, y(x) ), x );
# numeric solution:
sol := dsolve( {DE,BC}, y(x), numeric, output=listprocedure ):
Y := subs( sol, y(x) ):
# comparison exact and numeric sulution
plot( Y1-Y, 0..1 );

1 2 3 4 5 6 7 Last Page 3 of 11