Advanced Search

Maple 2024 Questions and Posts

These are Posts and Questions associated with the product, Maple 2024
View unanswered posts
Maple 2024 Questions and Posts Feed

Insights from My Recent Lecture on Euler Methods...

by: Paras31 245 Maple 2024
ode euler differential-equations
October 05 2024
1  1

Just finished an exciting lecture for undergraduate students on Euler methods for solving ordinary differential equations (ODEs)! 📝 We explored the Euler method and the Improved Euler method (a.k.a. Heun's method) and discussed how these fundamental techniques work for approximating solutions to ODEs.

To make things practical and hands-on, I demonstrated how to implement both methods using Maple—a fantastic tool for experimenting with numerical methods. Seeing these concepts come to life with visualizations helps understand the pros and cons of each method.

The Euler method is the simplest numerical approach, where we approximate the solution by stepping along the slope given by the ODE. But there's a catch: using a large step size can lead to large errors that accumulate quickly, causing significant deviation from the real solution.

Plot Results:

  • With a larger step size (h=0.5h = 0.5h=0.5), the solution tends to drift away significantly from the actual curve.
  • Reducing the step size improves accuracy, but it also means more steps and more computation.

Next, we discussed the Improved Euler method, which is like Euler's smarter sibling. Instead of blindly following the initial slope, it:

  1. Estimates the slope at the start.
  2. Uses that slope to predict an intermediate value.
  3. Then, it takes another slope at the intermediate point.
  4. Averages these two slopes for a better approximation.

This technique makes a big difference in accuracy and stability, especially with larger step sizes.

Plot Results:

  • Using the Improved Euler method, we found that even with a larger step size, the solution is smoother and closer to the true path.
  • The average slope helps mitigate the inaccuracies that arise from only using the beginning point's derivative, effectively reducing the local error.
  • The standard Euler method can produce solutions that oscillate or diverge, especially for larger step sizes.
  • The Improved Euler method follows the actual trend of the solution more faithfully, even with fewer steps. This makes it a more efficient choice for balancing computational effort and accuracy.
  • The Euler method is great for its simplicity but often requires very small step sizes to be accurate, leading to more computational effort.
  • The Improved Euler method—which we implemented and visualized on Maple—proved to be more reliable and accurate, especially for larger step sizes.

It was amazing to see the students engage with these foundational methods, and implementing them on Maple brought a deeper understanding of numerical analysis and the challenges of solving differential equations.


 


restart

with(plots)

with(DEtools)

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

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

 (1)

NULL

We calculate the general solution to the ODE

NULL

dsolve(ODE1, y(x))

y(x) = c__1/(x^2+1)

 (2)

NULL

Now let's solve the problem for the following two inital conditions

y(0) = 2 and y(0) = 1/2

dsolve({ODE1, y(0) = 2}, y(x))

y(x) = 2/(x^2+1)

 (3)

dsolve({ODE1, y(0) = 1/2}, y(x))

y(x) = 1/(2*x^2+2)

 (4)

Then, we are going to plot the solutions to the IVP's together with the slope field corresponding to the ODE.

 

NULLdfieldplot(ODE1, [y(x)], x = -5 .. 5, y = -5 .. 5, color = blue, scaling = constrained, axes = boxed)

  

 

DEplot(ODE1, y(x), x = -5 .. 5, y = -5 .. 5, [[y(0) = 2], [y(0) = 4]], linecolor = red, color = blue, scaling = constrained, axes = boxed)

  

Problem 2 (EULER METHOD)

NULL

NULL

ODE2 := diff(y(x), x) = sin(x*y(x))

diff(y(x), x) = sin(x*y(x))

 (5)

 

 

dsolve(ODE2, y(x))

NULL

Maple returned no output! That means Maple is unable to solve the equation.

NULL

If you are curious what steps Maple went through to  find a solution before failing
to do so, you can ask to see the steps using the command "infolevel". The levels
of information that you can request range from 0 to 5.

 

````

dsolve(ODE2, y(x))

soln := dsolve({ODE2, y(0) = 3}, y(x), numeric)

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 28, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..65, {(1) = 1, (2) = 1, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 21, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0, (64) = -1, (65) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = .16290418802201975, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..1, {(1) = 3.0}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..1, {(1) = .1}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8]), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..1, {(1) = 3.0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..1, {(1, 1) = .0, (2, 0) = .0, (2, 1) = .0, (3, 0) = .0, (3, 1) = .0, (4, 0) = .0, (4, 1) = .0, (5, 0) = .0, (5, 1) = .0, (6, 0) = .0, (6, 1) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = y(x)]`; YP[1] := sin(X*Y[1]); 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, 0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = y(x)]`; YP[1] := sin(X*Y[1]); 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 27 ) = (""), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0), ( 28 ) = (0) ] )) ] ); _y0 := Array(0..1, {(1) = 0.}); _vmap := array( 1 .. 1, [( 1 ) = (1) ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); _i := false; if _par <> [] then _i := `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then _i := `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) or _i end if; if _i then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch: end try else try procname(procname("right")) catch: end try end if end if end if; return elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch: end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch: end try end if end if end if; return elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if elif type(_xin, `=`) and lhs(_xin) = "setdatacallback" then if not type(rhs(_xin), 'nonegint') then error "data callback must be a nonnegative integer (address)" end if; _dtbl[1][28] := rhs(_xin) else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error end try; if _dat[17] <> _dtbl[1][17] then _dtbl[1][17] := _dat[17]; _dtbl[1][10] := _dat[10] end if; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; if type(_EnvDSNumericSaveDigits, 'posint') then _dat[4][26] := _EnvDSNumericSaveDigits else _dat[4][26] := Digits end if; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [x, y(x)], (4) = []}); _vars := _dat[3]; _pars := map(lhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol := 1; _ndsol := _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc

 (6)

[soln(1), soln(1.5), soln(2), soln(2.5), soln(3), soln(3.5), soln(4), soln(4.5), soln(5)]

[[x = 1., y(x) = HFloat(3.5280317971877286)], [x = 1.5, y(x) = HFloat(3.1226400064202693)], [x = 2., y(x) = HFloat(2.653970420106217)], [x = 2.5, y(x) = HFloat(2.323175669304077)], [x = 3., y(x) = HFloat(2.3019187052877816)], [x = 3.5, y(x) = HFloat(2.6587033583656714)], [x = 4., y(x) = HFloat(2.497876146260152)], [x = 4.5, y(x) = HFloat(2.220305976552359)], [x = 5., y(x) = HFloat(1.9758297601444128)]]

 (7)

p1 := odeplot(soln, [x, y(x)], x = 0 .. 5, color = magenta, thickness = 2, scaling = constrained, view = [0 .. 5, 0 .. 5])

  

p2 := dfieldplot(ODE2, [y(x)], x = 0 .. 5, y = 0 .. 5, color = blue, scaling = constrained)

display(p1, p2)

  

DEplot(ODE2, y(x), x = 0 .. 5, y = 0 .. 5, [[y(0) = 3]], linecolor = red, color = blue, scaling = constrained, axes = boxed)

  

Construct approximate solutions for x from 0 to 5 to the initial problem 2 using Euler's method with the three different step sizes Δx=0.5, 0.25, 0.125

f := proc (x, y) options operator, arrow; sin(x*y) end proc

proc (x, y) options operator, arrow; sin(y*x) end proc

 (8)

Eulermethod := proc (f, x_start, y_start, dx, n_total) local x, y, Y, i, k; x[1] := x_start; y[1] := y_start; for i to n_total do y[i+1] := y[i]+f(x[i], y[i])*dx; x[i+1] := x[1]+i*dx end do; Y := [seq([x[k], y[k]], k = 1 .. n_total)]; return Y end proc

proc (f, x_start, y_start, dx, n_total) local x, y, Y, i, k; x[1] := x_start; y[1] := y_start; for i to n_total do y[i+1] := y[i]+f(x[i], y[i])*dx; x[i+1] := x[1]+i*dx end do; Y := [seq([x[k], y[k]], k = 1 .. n_total)]; return Y end proc

 (9)

NULL

NULL

NULL

NULL

A1 := Eulermethod(f, 0, 3, .5, 10)

[[0, 3], [.5, 3.], [1.0, 3.498747493], [1.5, 3.323942476], [2.0, 2.842530099], [2.5, 2.560983065], [3.0, 2.620477947], [3.5, 3.120464063], [4.0, 2.621830593], [4.5, 2.185032319]]

 (10)

A2 := Eulermethod(f, 0, 3, .25, 20)

[[0, 3], [.25, 3.], [.50, 3.170409690], [.75, 3.420383740], [1.00, 3.556616077], [1.25, 3.455813236], [1.50, 3.224836021], [1.75, 2.976782400], [2.00, 2.757025820], [2.25, 2.583147563], [2.50, 2.469680146], [2.75, 2.442487816], [3.00, 2.547535655], [3.25, 2.791971512], [3.50, 2.877900373], [3.75, 2.727027361], [4.00, 2.547414295], [4.25, 2.374301829], [4.50, 2.219839448], [4.75, 2.086091178]]

 (11)

A3 := Eulermethod(f, 0, 3, .125, 40)

[[0, 3], [.125, 3.], [.250, 3.045784066], [.375, 3.132030172], [.500, 3.247342837], [.625, 3.372168142], [.750, 3.479586272], [.875, 3.542983060], [1.000, 3.548166882], [1.125, 3.498733738], [1.250, 3.409546073], [1.375, 3.297015011], [1.500, 3.174012084], [1.625, 3.049159854], [1.750, 2.927817142], [1.875, 2.813241461], [2.000, 2.707496816], [2.125, 2.612101612], [2.250, 2.528513147], [2.375, 2.458549923], [2.500, 2.404840953], [2.625, 2.371369082], [2.750, 2.364080535], [2.875, 2.391119623], [3.000, 2.460798019], [3.125, 2.572154041], [3.250, 2.695044010], [3.375, 2.772263417], [3.500, 2.780805371], [3.625, 2.742906337], [3.750, 2.680985435], [3.875, 2.607451728], [4.000, 2.528940353], [4.125, 2.449278434], [4.250, 2.370825619], [4.375, 2.295054886], [4.500, 2.222824117], [4.625, 2.154537647], [4.750, 2.090275081], [4.875, 2.029905439]]

 (12)

NULL

p3 := pointplot(A1, color = green, scaling = constrained, symbol = circle)

p4 := pointplot(A2, color = black, scaling = constrained, symbol = asterisk)

p5 := pointplot(A3, color = red, scaling = constrained, symbol = diamond)

NULL

display([p1, p3, p4, p5])

  

NULL

In this plot, the differences between the lines visually represent how using different step sizes affects the overall solution accuracy. The red points are closest to what a more accurate numerical solution would look like, while the green points are more of a rough approximation.

 

 

Problem 3 (IMPROVED EULER METHOD

``

 

 

 

ImprovedEulermethod := proc (f, x_start, y_start, dx, n_total) local x, y, Y, i, k, k1, k2; x[1] := x_start; y[1] := y_start; for i to n_total do k1 := f(x[i], y[i]); k2 := f(x[i]+dx, y[i]+k1*dx); y[i+1] := y[i]+(1/2)*dx*(k1+k2); x[i+1] := x[1]+i*dx end do; Y := [seq([x[k], y[k]], k = 1 .. n_total+1)]; return Y end proc

proc (f, x_start, y_start, dx, n_total) local x, y, Y, i, k, k1, k2; x[1] := x_start; y[1] := y_start; for i to n_total do k1 := f(x[i], y[i]); k2 := f(x[i]+dx, y[i]+k1*dx); y[i+1] := y[i]+(1/2)*dx*(k1+k2); x[i+1] := x[1]+i*dx end do; Y := [seq([x[k], y[k]], k = 1 .. n_total+1)]; return Y end proc

 (13)

A := ImprovedEulermethod(f, 0, 3, .5, 10); B := ImprovedEulermethod(f, 0, 3, .25, 20); C := ImprovedEulermethod(f, 0, 3, .125, 40)

[[0, 3], [.125, 3.022892033], [.250, 3.089335383], [.375, 3.190999573], [.500, 3.311460714], [.625, 3.426126738], [.750, 3.508312296], [.875, 3.539992283], [1.000, 3.518184043], [1.125, 3.451931748], [1.250, 3.354946855], [1.375, 3.240089444], [1.500, 3.117181611], [1.625, 2.992925708], [1.750, 2.871597591], [1.875, 2.755823574], [2.000, 2.647215450], [2.125, 2.546845751], [2.250, 2.455612702], [2.375, 2.374560360], [2.500, 2.305227222], [2.625, 2.250113264], [2.750, 2.213376654], [2.875, 2.201814459], [3.000, 2.225589190], [3.125, 2.295322044], [3.250, 2.406184220], [3.375, 2.516909932], [3.500, 2.583378006], [3.625, 2.599915321], [3.750, 2.579967129], [3.875, 2.537160528], [4.000, 2.481052245], [4.125, 2.417781798], [4.250, 2.351265142], [4.375, 2.284028571], [4.500, 2.217702881], [4.625, 2.153313287], [4.750, 2.091461577], [4.875, 2.032452273], [5.000, 1.976386380]]

 (14)

NULL

plot2 := pointplot(A, color = green, scaling = constrained, symbol = circle)

plot3 := pointplot(B, color = black, scaling = constrained, symbol = asterisk)

plot4 := pointplot(C, color = red, scaling = constrained, symbol = diamond)

NULL

display([p1, plot2, plot3, plot4])

  

The Improved Euler method, by averaging slopes, provides a significant improvement over the basic Euler method, particularly when the step size is relatively large.

Download Diff_equations.mw

square side wanted...

Asked by: Alfred_F 325
October 05 2024
1  1

Given a conventionally labeled triangle ABC with the two sides a=3 and b=4, c is not given. What is the side length of the largest inscribed square whose "base" lies on the triangle side c?

Little ODE for Puzzles...

Asked by: Alfred_F 325
October 04 2024
0  0

The usual ODE must be solved:
y´´*(y^3-y)+y´^2 *(y^2+1)=0
"Dangerous places" of the definition domain must be described: Where are the general solution y(x) and its derivatives continuous?

how to obtain the plot command used after updating...

Asked by: nm 11353
October 03 2024
3  5

This is probably asked before but can't find it. 

After making a plot, then RIGHT-CLICKing on it, option menu comes up that allows one to modify the plot (like adding gridlines, or change the line style).

How does one find the Maple command after doing such changes, so one can use the command in the code?

Here is an example. I modified a little acers plot in this answer and added the points to the plot

eq:=-0.0004*x^2 - 2.7680*10^(-28)*x^12 - 2.1685*10^(-43)*x^18 - 1.3245*10^(-37)*x^16 - 1.6917*10^(-32)*x^14 + 0.7650 + 6.6773*10^(-18)*x^8 - 2.5543*10^(-23)*x^10 - 8.0002*10^(-13)*x^6 + 3.6079*10^(-8)*x^4 = 0:
the_roots:=fsolve(eq):	
the_roots:=map(X->[X,0],[the_roots]):
p1:=plot(lhs(eq),x=-210..210,size=[500,200]):
p2:=plot(the_roots,style=point, color=red, symbol=solidcircle, symbolsize=20):
plots:-display(p1,p2)

 

I was lazy to look up the grid option syntax. So right clicked on the plot and found option to add gridline. So said, great. And the above is the result.

But now I'd like to see the command used so I know where the grid option goes and how its syntax is. There does not seem to be option in the interface to display the Maple command used.

How would one find it?

Maple 2024.1 on windows 10

 

Cycloidal curve

by: sija 165 Maple 2024
plotting rolling cycloid
October 01 2024
2  0

Hello,

Below I am attaching a Maple worksheet on rolling a circle on a flat curve and generating a cycloidal curve.

Regards!

Cycloidal_curve.mw



 

What is the correct way to check if one number is ...

Asked by: nm 11353
programming equality evalb
October 01 2024
2  9

What is the correct idiom in Maple to check that one number is greater than another or not?

For example if n=sqrt(3) and m=3 ?

Maple does not like to compare sqrt(3) with other number. It says 

   Error, cannot determine if this expression is true or false: 3 <= 3^(1/2)

I am doing this in code, so solution has to be such that it works for all cases of n and m, without the ability to look at screen and decide.

I found that using is works without the need to convert to float.  Should one then use is(n>=m) instead of evalb(n>=m) always?  

I can also always apply evalf() on each side before. But this seems like an overkill to me.

> 

restart;
> 

n:=sqrt(3);
m:=3;

3^(1/2)

3

> 

if evalb(n>=m) then
   "yes";
else
   "no";
fi;

Error, cannot determine if this expression is true or false: 3 <= 3^(1/2)

> 

if n>=m then
   "yes";
else
   "no";
fi;

Error, cannot determine if this expression is true or false: 3 <= 3^(1/2)

> 

if evalb(evalf(n)>=m) then
   "yes";
else
   "no";
fi;

"no"

> 

if is(n>=m) then
   "yes";
else
   "no";
fi;

"no"

 

 

Download checking_one_number_larger_than_another.mw

For reference, another software does compare these two numbers as is without the need to convert them to floating point;

But it looks the design of Maple in this aspect is different. Since the principal root of sqrt(3) is always the positive one, I did not think it will cause a problem.

In memory of IMO 1988...

Asked by: Alfred_F 325
September 30 2024
3  2

A task that was famous at the time is worth remembering:

If for whole numbers x and y the number N = (x^2+y^2)/(1+x*y) is a positive whole number, then it is also a square number.

It can be proven that the converse is also true. Therefore, here is the task:

If N is a square number, then the Diophantine equation has solutions. Solutions must be calculated for N = 9, 49 and 729.

Avoiding placeholders for symbolic differentiation...

Asked by: Paras31 245
differentiation differential-equations
September 28 2024
0  3

I’m currently working on a Maple project involving the symbolic differentiation of a user-defined potential energy function, which I then use in a system of differential equations. To implement this, I've been using placeholders like 'Un' when taking derivatives to substitute later during numerical calculations.

While this approach using placeholders works, I find it cumbersome and would prefer a more straightforward way to handle this symbolic differentiation without having to rely on intermediate placeholders like 'Un'.

Does anyone know if there is a more elegant way in Maple to achieve this functionality?

Specific Requirements:

  • I want to symbolically differentiate the user-defined function W directly.
  • The goal is to use this derivative in a system of differential equations for numerical solving.

Any suggestions on simplifying this process would be much appreciated!

Thanks in advance.

restart

NULL
W := proc (Un) options operator, arrow; -(1/4)*(b*d-e*q)^2/(b^2*(e+b*Un^2))-(2*b*c*m^2+(1/4)*q^2)*Un^2/b+c*Un^4 end proc

proc (Un) options operator, arrow; -(1/4)*(b*d-e*q)^2/(b^2*(e+b*Un^2))-(2*b*c*m^2+(1/4)*q^2)*Un^2/b+c*Un^4 end proc

 (1)

 

 

diff(W(Un), Un)

(1/2)*(b*d-e*q)^2*Un/(b*(Un^2*b+e)^2)-2*(2*b*c*m^2+(1/4)*q^2)*Un/b+4*c*Un^3

 (2)

b := -0.3070816340e-3; c := .4461893869; d := 0.1142581922e-1; e := 0.7675000601e-3; q := -0.3704049149e-2; m := 1.423206983; plot(W(Un), Un = -10 .. 10, title = "Graph of Potential W for set 1", color = blue, axes = boxed)

  

with(plots); with(DEtools); delta := .5; b := -0.222159366e-3; c := 1.088046084; d := 0.1308858509e-2; e := 0.394423507e-3; q := -0.64844084e-3; m := 1.518466566; W := proc (Un) options operator, arrow; -(1/4)*(b*d-e*q)^2/(b^2*(e+b*Un^2))-(2*b*c*m^2+(1/4)*q^2)*Un^2/b+c*Un^4 end proc; Un_placeholder := 'Un'; Wprime := diff(W(Un_placeholder), Un_placeholder); sys := {diff(u(t), t) = v(t), diff(v(t), t) = -delta*v(t)+subs(Un_placeholder = u(t), Wprime)}; ics1 := {u(0) = 0.1e-1, v(0) = 0.1e-1}; sol := dsolve(`union`(sys, ics1), numeric, output = listprocedure); xrange := -0.10e-1 .. 0.11e-1; yrange := -0.3e-1 .. 0.3e-1; phase_plot := odeplot(sol, [u(t), v(t)], 0 .. 100, title = "Phase Plot (Velocity vs Displacement)", numpoints = 8000, color = red, labels = ["u(t)", "v(t)"], axes = boxed, labeldirections = [horizontal, vertical]); vectorfield_plot1 := fieldplot([v, -delta*v+subs(Un_placeholder = u, Wprime)], u = xrange, v = yrange, arrows = medium, color = blue, axes = boxed); display([vectorfield_plot1, phase_plot], title = "Combined Phase Plot and Vector Field", labels = ["u(t)", "v(t)"], axes = boxed)

  

NULL

Download proteins.mw

Check if two vectors are equal?...

Asked by: Ronan 1341
parallel vector linear-algebra
September 28 2024
0  6

What on earth am I doing wrong here? I am trying to check if two vectors are equal for parallelism.Yes, thinking about it I could use CrossProduct but that is the same problem because I would be checking if that new vector =<0,0,0>

l1:=<a,b,c>:l2:=<a,b,c>:
if l1 = l2 then 
print("Equal Vectors");
end if;

"isolve" does not solve...

Asked by: Alfred_F 325
isolve diophantine_equation
September 28 2024
1  7

While cleaning up old documents, I found the following Diophantine equation and its solution. I tried to solve it using Maple´s "isolve" but it didn't work. Please give me some advice.
Equation:
x^2 - 12*x*y + 6*y^2 + 4*x + 12*y - 3=0

ODESteps gives improper op or subscript selector...

Asked by: nm 11353
ode step-by-step differential-equations
September 26 2024
2  2

Why Maple gives

           improper op or subscript selector

On this ode? I was not expecting to see steps for this, but was just trying to see what it will show.

> 

restart;
> 

interface(version);

`Standard Worksheet Interface, Maple 2024.1, Windows 10, June 25 2024 Build ID 1835466`

> 

Physics:-Version();

`The "Physics Updates" version in the MapleCloud is 1810 and is the same as the version installed in this computer, created 2024, September 18, 18:16 hours Pacific Time.`

> 

libname;

"C:\Users\Owner\maple\toolbox\2024\Physics Updates\lib", "C:\Program Files\Maple 2024\lib"

> 

restart;
> 

ode:=sin(x)*diff(y(x),x$2)+cos(x)*diff(y(x),x)+(sin(x)-cos(x))*y(x)=0;

sin(x)*(diff(diff(y(x), x), x))+cos(x)*(diff(y(x), x))+(sin(x)-cos(x))*y(x) = 0

> 

dsolve(ode);

y(x) = c__1*(tan((1/2)*x)+I)^(1/2+(1/2)*(5-4*I)^(1/2))*(tan((1/2)*x)-I)^((1/2)*(5+4*I+(5+4*I)^(1/2))/(5+4*I)^(1/2))*HeunG(2, 2*(((3/4)*(5-4*I)^(1/2)+7/4-I)*(5+4*I)^(1/2)+5/4+I+(5/4+I)*(5-4*I)^(1/2))/(5+4*I)^(1/2), 1+(1/2)*(10+2*41^(1/2))^(1/2), (1/2)*((5+4*I)^(1/2)*(5-4*I)^(1/2)+5+4*I+2*(5+4*I)^(1/2))/(5+4*I)^(1/2), 1+(5-4*I)^(1/2), 1, 1-I*tan((1/2)*x))+c__2*(tan((1/2)*x)+I)^(1/2-(1/2)*(5-4*I)^(1/2))*(tan((1/2)*x)-I)^((1/2)*(5+4*I+(5+4*I)^(1/2))/(5+4*I)^(1/2))*HeunG(2, -2*(((3/4)*(5-4*I)^(1/2)+(-7/4+I))*(5+4*I)^(1/2)+(-5/4-I)+(5/4+I)*(5-4*I)^(1/2))/(5+4*I)^(1/2), -(1/2)*(5-4*I)^(1/2)+1+(1/2)*(5+4*I)^(1/2), (1/2)*(-(5+4*I)^(1/2)*(5-4*I)^(1/2)+5+4*I+2*(5+4*I)^(1/2))/(5+4*I)^(1/2), 1-(5-4*I)^(1/2), 1, 1-I*tan((1/2)*x))

> 

Student:-ODEs:-ODESteps(ode)

Error, (in Student:-ODEs:-OdeSolveOrder2) improper op or subscript selector

> 

 


 

Download internal_error_ODEsteps_sept_26_2024.mw

Differences in Domain Results Between solve(0 < x ...

Asked by: Paras31 245
solve inequality
September 25 2024
2  2

I recently encountered an interesting issue when solving for the domain of a function in Maple using two different approaches, and I am curious about the differences in results between the following commands:

  1. domain := solve(0 < x and x <> 2, x);
  2. domain := solve({0 < x, x <> 2}, x);

I wonder if anyone could explain why Maple might treat these two syntaxes differently.

> 

restart
> 

``
> 

h := proc (x) options operator, arrow; ln(x)/(x-2) end proc

proc (x) options operator, arrow; ln(x)/(x-2) end proc

 (1)
> 

domain_ln := solve(x > 0, x)

RealRange(Open(0), infinity)

 (2)
> 

denom_undefined := solve(x-2 = 0, x)

2

 (3)
> 

domain := solve(`and`(x > 0, x <> 2), x)

RealRange(Open(0), infinity)

 (4)
> 

solve({x > 0, x <> 2}, x)

{2 < x}, {0 < x, x < 2}

 (5)
> 

``
> 

plot(h(x), x = 0 .. 10, discont = true, color = blue)

 
> 

``


 

Download domain.mw

fsolve within loop returns unevaluated...

Asked by: smart2718 0
September 25 2024
1  5

I am new to Maple will appreciate any advice.

I am trying to write a procedure to solve a cubic equation using fsolve, print the root at each stage, and store results in a table for later use.  

I have restricted the range of fsolve, but I know from other sources that there is a solution within this range.

The instructions work if entered manually line by line, but not as a procedure.

The code is copied below.  Thanks for any suggestions.

Stephen Martin

Reg_IV := proc (beta, n) local alpha, b, c, d, k, `&alpha;L`, 

   `&alpha;R`, delta, eqn, x; `&alpha;L` := 3+beta-6*beta^(1/2);\

   `&alpha;R` := 1-beta; delta := (`&alpha;R`-`&alpha;L`)/n; 

   print(`&beta; = `, evalf(beta, 5), `n = `, evalf(n, 5)); 

   print(`&alpha;L = `, evalf(`&alpha;L`, 5), ` &alpha;R = `, 

   evalf(`&alpha;R`, 5), `&delta; = `, evalf(delta, 5)); 

   print(``); for k from 0 to n+1 do alpha := `&alpha;L`+k*delta\

  ; b := -5+alpha-3*beta; c := 8-4*alpha-6*beta+2*beta^2-2*alpha\

  *beta; d := -(4*beta+4)*(1-alpha-beta); eqn := proc (b, c, d, 

   x) options operator, arrow; {x^3+b*x^2+c*x+d} end proc; x[k] 

   := fsolve(eqn, 0 .. 2); print(evalf(x[k], 5)) end do; 

   print(``); print(evalf(x, 5)) end proc


Reg_IV(0.2, 4);
                    &beta; = , 0.2, n = , 4.

 &alpha;L = , 0.51672,  &alpha;R = , 0.8, &delta; = , 0.070820


                      fsolve(eqn, 0 .. 2)

                      fsolve(eqn, 0 .. 2)

                      fsolve(eqn, 0 .. 2)

                      fsolve(eqn, 0 .. 2)

                      fsolve(eqn, 0 .. 2)

                      fsolve(eqn, 0 .. 2)


   TABLE([0 = fsolve(eqn, 0 .. 2), 1 = fsolve(eqn, 0 .. 2), 

     2 = fsolve(eqn, 0 .. 2), 3 = fsolve(eqn, 0 .. 2), 

     4 = fsolve(eqn, 0 .. 2), 5 = fsolve(eqn, 0 .. 2)])


 

Next one evalc and abs problem...

Asked by: specter123 25
complex syntax abs
September 24 2024
1  2

Verification_24_09_2024.mw

In the attached file, in the subsection named #Transmittance calculation I had an issue with the absolute value calculation. The calculation in the first part for the expression a1 was carried out although with addition of some strange t parameter in the exponent power. The second one just didn`t. I will be more than grateful for some help and advice. In addition, I have been struggling with that for a long time so if you have in mind some book or youtube lesson that clarifies all of that I will be more than happy. Nothing valuable was found on the web so far :(

Plotting intersection point of two equations fails...

Asked by: Paras31 245
syntax substitution pointplot
September 23 2024
1  1

Hello everyone,

I successfully plotted the intersection point of two linear equations using a matrix-based approach. Still, I’m facing an issue when trying to achieve the same result using the second method — solving the system symbolically and plotting without matrices.

The Matrix-Based Approach (Working): I used the following steps to plot both equations and mark the intersection point (solution). However, when I switch to solving the system symbolically without matrices, I can plot the two lines, but I'm unable to plot the intersection point correctly.

It seems that the symbolic solution isn't being properly converted to a numeric form for pointplot, even though I’m using evalf.

Any help would be greatly appreciated!

Download Linear_System.mw

First 22 23 24 25 26 27 28 Last Page 24 of 43