## 1350 Reputation

15 years, 9 days
University of Twente (retired)
Enschede, Netherlands

My "website" consists of a Maple Manual in Dutch

## Alternative...

@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)");`

## Newton-Raphson...

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

## Newton-Raphson...

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

## Function of x and omega...

@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.

## Function of x and omega...

@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.

## Subtle distinction.....

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);
 (1)
 > b := 4^(1/2);
 (2)
 > c := 4^.5;
 (3)
 > 4.^(1/2);
 (4)
 > if a=b then "a=b" else "a<>b" end if;
 (5)
 > if a=c then "a=c" else "a<>c" end if;
 (6)
 > if b=c then "b=c" else "b<>c" end if;
 (7)
 >

## Subtle distinction.....

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);
 (1)
 > b := 4^(1/2);
 (2)
 > c := 4^.5;
 (3)
 > 4.^(1/2);
 (4)
 > if a=b then "a=b" else "a<>b" end if;
 (5)
 > if a=c then "a=c" else "a<>c" end if;
 (6)
 > if b=c then "b=c" else "b<>c" end if;
 (7)
 >

## Perhaps it is a bug?...

`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
```

## Perhaps it is a bug?...

`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
```

## No error......

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

## No error......

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

## Of course I did!...

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

## More transparent code...

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:`

## More transparent code...

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:`

## Compare with exact solution:...

(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 sulutionplot( Y1-Y, 0..1 );`

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