Items tagged with odeplot odeplot Tagged Items Feed

I'm trying to plot the varying results to a second degree differential function with different values for one constant in one graph. Here is what I have so far, which is already working.


> with(plots);
> m := 0.46e-1; d := 0.42e-1; v := 60; alpha0 := convert(12*degrees, radians); g := 9.81; pa := 1.205; cd := .2; n := 4000; omega := 2*Pi*(1/60);
> p := 6*m/(Pi*d^3);
> k1 := (3/4)*cd*pa/(d*p); k2 := (3/8)*omega*n*pa/p;
> gl1 := vx(t) = diff(x(t), t);
> gl2 := vy(t) = diff(y(t), t);
> gl3 := diff(vx(t), t) = -k1*vx(t)*(vx(t)^2+vy(t)^2)^(1/2)-k2*vy(t);
> gl4 := diff(vy(t), t) = -g-k1*vy(t)*(vx(t)^2+vy(t)^2)^(1/2)+k2*vx(t);
> init1 := x(0) = 0;
> init2 := y(0) = 0;
> init3 := vx(0) = v*cos(alpha0);
> init4 := vy(0) = v*sin(alpha0);
> sol := dsolve({gl1, gl2, gl3, gl4, init1, init2, init3, init4}, {vx(t), vy(t), x(t), y(t)}, type = numeric);

> sol(.5);
> odeplot(sol, [x(t), y(t)], t = 0 .. 6.7);

What I'd like to do now is, for example, plot the solutions in one graph (preferably as a gif) for when n=1500, n=3000, n=4500 etc. Is there a simple way to achieve this? I've tried various methods so far without success.

I want to plot the position of Spherical pendulum. there are differential equation for spherical pendulum in spherical coordinates

sys := {((D@@2)(phi))(t) = -2*(D(phi))(t)*(D(theta))(t)*cos(theta(t))/sin(phi(t)),
((D@@2)(theta))(t) = (D(phi))(t)^2*cos(theta(t))*sin(theta(t))-9.8*sin(theta(t))}


with initial conditions

theta(0) = (1/2)*Pi, (D(theta))(0) = 0, phi(0) = (1/2)*Pi, (D(phi))(0) = 1

I tried:

eq := dsolve([((D@@2)(theta))(t) = (D(phi))(t)^2*cos(theta(t))-9.8*sin(theta(t)),
((D@@2)(phi))(t) = -2*(D(phi))(t)*(D(theta))(t)*cos(theta(t))/sin(phi(t)),
theta(0) = (1/2)*Pi, (D(theta))(0) = 0, phi(0) = (1/2)*Pi, (D(phi))(0) = 1],numeric)

how to change coordinates

x(t) = sin(theta(t))*cos(phi(t))
y(t) = sin(theta(t))*sin(phi(t))
z(t) = cos(theta(t))

and how to plot it from t=0 to 10?

I have a set of differential equations on 3 variables, B[1],B[2] and C. Its not physically meaningful for B[1]+B[2]>0.5 so i would ideally like to replace the cube that the solutions are displayed on (the axis take the limits B[1]=0...0.5,B[2]=0...0.5,C=0...100 ) with a triangular prism (the axis take the limits B[1]=0...0.5,B[2]=0...0.5,B[1]+B[2]<0.5,C=0...100 ).

Failing that i'd like the plot to display a plane showing where the meaningful values for the variables end.

here is the code I use to put the ODEplot together

Model := [diff(B[1](t), t) = k[a1]*C(t)*(R-B[1](t)-B[2](t))-k[d1]*B[1](t), diff(B[2](t), t) = k[a2]*C(t)*(R-B[1](t)-B[2](t))-k[d2]*B[2](t), diff(C(t), t) = (-(k[a1]+k[a2])*C(t)*(R-B[1](t)-B[2](t))+k[d1]*B[1](t)+k[d2]*B[2](t)+k[m]*((I)(t)-C(t)))/h];
DissMod := subs((I)(t) = 0, Model);
AssMod := subs((I)(t) = C[T], Model);

Pars1a := [k[a1] = 6*10^(-4), k[d1] = 7*10^(-3), k[a2] = 6*10^(-4), k[d2] = (7/5)*10^(-3), R = .5, k[m] = 10^(-4), C[T] = 100, h = 10^(-6)];

Pars := Pars1a; thing11 := subs(Pars, AssMod[1]), subs(Pars, AssMod[2]); thing12 := diff(C(t), t) = piecewise(t <= 100, subs(Pars, rhs(AssMod[3])), subs(Pars, rhs(DissMod[3]))); sol := dsolve({thing11, thing12, C(0) = 0, B[1](0) = 0, B[2](0) = 0}, {C(t), B[1](t), B[2](t)}, numeric, output = listprocedure, maxstep = 2, maxfun = 1000000); ParsPlot1a := odeplot(sol, [B[1](t), B[2](t), C(t)], t = 0 .. 700, color = blue, view = [0 .. .5, 0 .. .5, 0 .. 100], tickmarks = [[0 = 0, .5 = R], [0 = 0, .5 = R], [0 = 0, 25 = (1/4)*C[T], 50 = (1/2)*C[T], 100 = C[T]]]);

But I can't see a way of either making the plane or making the ODEplot axis into something other than a cube.

I have some differential equations that I want to plot on the same axis (as i have below), but would like to plot with a log scale to illustrate the fact that one is simply a logarithmic decay (solid line) and all the others are not.

I had used deplot and display to make the above graph,

but if you want more detail, here is a worksheet with the relevant differntial equations:




Hi everyone!

I tried to plot the solution of the following ode, but I only got the message error:

Warning, cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

(file attached)


Please, help me!


Thank you so much!


thanks. I played around, and had problems implementing your ideas for one of the systems I'm interested in.I don't see a difference between this and what you had advised me on, but it gets an error.

any idea why?
or how to fix it?

thing1 := diff(B[1](t), t) = piecewise(t <= 500, 0.3e-2-(63/10000)*B[1](t)-(3/500)*B[2](t), -(3/10000)*B[1](t)):
thing2 := diff(B[1](t), t) = piecewise(t <= 500, 0.1e-1-(1/50)*B[1](t)-(13/625)*B[2](t), -(1/1250)*B[2](t)):
sol := dsolve({thing1, thing2, B[1](0) = 0, B[2](0) = 0}, {B[1](t), B[2](t)}, numeric, output = listprocedure); plots:-odeplot(sol, [B[1](t), B[2](t)], t = 450 .. 550);

Error, (in dsolve/numeric/DAE/explicit) unable to obtain the standard form of the DAE system due to the presence of leading dependent variables/derivatives in the piecewise: piecewise(t <= 500, 1/100-(1/50)*B[1](t)-(13/625)*B[2](t), -(1/1250)*B[2](t))-piecewise(t <= 500, 3/1000-(63/10000)*B[1](t)-(3/500)*B[2](t), -(3/10000)*B[1](t))
Error, (in plots/odeplot) curve is not fully specified in terms of the ODE solution, found additional unknowns {B[1](t), B[2](t)}


I would like to plot gait diagrams (the lines you can see on the picture belowà from the solutions obtained with a NL oscillator (composed with 8 coupled odes). Here the result that I would like to obtain.

Initial plot:


Desired plot



I would like to obtain 4 lines corresponding to the 4 elliptic trajectories obtained with the NL oscillator. The four lines should be done like this. When the trajectory is above 0, the line should be colored in green. When the trajectory is below 0, the line should be colored in black. 

May you help me to define this kind of graph called gait diagrams from the solution of the NL oscillator ?

Here you can find my maple code:

K:=Matrix([<0, -1, 1, -1>,<-1, 0, -1, 1>,<-1, 1, 0,-1>,<1, -1, -1,0>]);

for i to 4
end do:


Differential system 
ic:=[u[1](0)=0, v[1](0)=0,u[2](0)=0, v[2](0)=-0.1,u[3](0)=0, v[3](0)=0.1,u[4](0)=0, v[4](0)=0.1];
Initial boundaries
ic2:=[seq(u[i](0)=eval(u[i](t), res(tcalc)),i=1..4),seq(v[i](0)=eval(v[i](t), res(tcalc)),i=1..4)];

tmax:= 40:
plots:-odeplot(res,[t,v[1](t)],0..tmax,thickness=2, view=[0..5, -1.5..1.5],numpoints = numpts);
plots:-odeplot(res,[t,v[2](t)],0..tmax,thickness=2, view=[0..5, -1.5..1.5],numpoints = numpts);
plots:-odeplot(res,[t,v[3](t)],0..tmax,thickness=2, view=[0..5, -1.5..1.5],numpoints = numpts);
plots:-odeplot(res,[t,v[4](t)],0..tmax,thickness=2, view=[0..5, -1.5..1.5],numpoints = numpts);
plots:-odeplot(res,[seq([t,v[i](t)+i*5],i=1..4)],0..tmax,thickness=2,view=[0..5,0..25], numpoints = numpts);

Thanks a lot for your help

Dear all;

Thank you for helping me to solve this  question.

I solve an ode, but I have an error when I would like to plot the solution.

uanble to achieve continuous solution with requested accuracy of 0.1e-5 with maximum 128 point mesh (was able to get 0.14e-5), consider increasing `maxmesh` or using larger `abserr`
Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

I try to increase the point mesh or take a large abserr but always I have the same problem.



ode := diff(y(x), x, x) = x*y(x)+sqrt(x);

ics := y(0) = 0, y(1) = 1;
sol:=dsolve({ode,ics}, numeric):
odeplot( sol,[x, y(x)], x=0..1, maxmesh=1000);

hi guys,


i have a first order differential equation and i want to plot it with odeplot in polar coordinates .



thanks in

I am exploring the onset of chaos in the driven damped pendulum. I have set up the non-linear initial value problem that describes the situation in terms of phi(t), the angle of the pendulum from the vertical. I can solve the differential equation numerically for certain initial conditions using dsolve, and then graph that solution using odeplot. I can then explore what happens when I change the initial conditions ever so slightly. Again I can solve this new initial value problem numerically with dsolve, and then graph it with odeplot, Suppose the solution to the first problem is phi1(t), and the solution to the second initial value problem phi2(t). I would like to plot log | phi2(t) - phi1(t) | versus t, that is, the logarithm of the absolute value of the difference of the two solutions versus t. I’m interested in how sensitive the system is to initial conditions.

So, how can I disentangle the two solutions, phi1(t) and phi2(t), from dsolve and then plot the graph I want?


I've been trying to collect points at certain intervals in my plot, however, the sample option doesn't seem to be working. Can anyone lend me a hand? Any help is greatly appreciated! Thank you in advance!

Kind regards,



a := 1.501*10^9:

Th := sqrt(4*Pi^2*a^3/(G*(Mh+Msat))):

HyperionOrbit := proc (`&theta;IC`, `&omega;IC`, n) local a, Mh, Msat, G, e, beta, M, Eqns, ICs, soln; global `&omega;H`, Th; a := 1.501*10^9; Mh := 5.5855*10^18; Msat := 5.6832*10^26; G := 6.67259/10^11; e := .232; beta := .89; M := Mh+Msat; Eqns := diff(theta(t), t) = omega(t), diff(omega(t), t) = -G*Msat*beta^2*(xH(t)*sin(theta(t))-yH(t)*cos(theta(t)))*(xH(t)*cos(theta(t))+yH(t)*sin(theta(t)))/(xH(t)^2+yH(t)^2)^2.5, diff(xH(t), t) = vxH(t), diff(vxH(t), t) = -G*M*xH(t)/(xH(t)^2+yH(t)^2)^(3/2), diff(yH(t), t) = vyH(t), diff(vyH(t), t) = -G*M*yH(t)/(xH(t)^2+yH(t)^2)^(3/2); ICs := xH(0) = a*(1+e), yH(0) = 0, vxH(0) = 0, vyH(0) = sqrt(G*M*(1-e)/(a*(1+e))), theta(0) = `&theta;IC`, omega(0) = `&omega;IC`; soln := dsolve({Eqns, ICs}, numeric, maxfun = 0); plots:-odeplot(soln, [theta(t), omega(t)/`&omega;H`], 0 .. n*Th, numpoints = 2000, labels = ["&theta;(t)","&omega;(t)/&omega;H"], axes = boxed, style = point, sample = [seq(i, i = 0 .. n*Th, Th)]) end proc:

HyperionOrbit(.5, 1.8*`&omega;H`, 10)

Error, (in plot/options2d) unexpected option: sample = [0, 1876321.326, 3752642.652, 5628963.978, 7505285.304, 9381606.630, 11257927.96, 13134249.28, 15010570.61, 16886891.93, 18763213.26]






bia Man

Suppose we have a differential equation that it's order is 2. For example

diff(x(t), t, t)+(1000*(x(t)^2-1))*(diff(x(t), t))+x(t) = 0,

and we want to solve it numerically. When we use some initial values and dsolve it, the maple solves it very quickly. We can plot each new equations of x(t) without any problem respect of "t" or else.

But when we want to show the variation of "diff(x(t), t, t)" respect of "t" (when our derivation order is equal or upper than our differential equation order) our answer is 0 however dx(t)/dt is not 0 or constant!

What's happen here?

I tested it with all kind of methods as "rkf45","rosenbrock","lsode" etc.  but non of them showed me a correct answer.

How can I solve it correctly?
Please help me!



I'm writing a code and I seem to have an issue when trying to implement a procedure. Here is the code:



Z := 75; A := 189; k := 14.6; Rm := 8*R; r0 := 10^(-8)*R; c := 137.036; ms := 105.66/(.51100)

fmtoau := 10^(-15)/(0.529177e-1*10^(-9)):

R := 1.1*fmtoau*A^(1.0/(3.0)):

f := proc (x) options operator, arrow; 1/(1+exp(k*(x-1))) end proc:

n0 := 3*Z*k^2/(4*Pi*(Pi^2+k^2)*R^3):

n := proc (r) options operator, arrow; 4*Pi*n0*f(r/R) end proc:

int(r^2*n(r), r = 0 .. Rm);



plot(n/n0, 0 .. 2*R);


v1 := unapply(int(x^2*f(x), x), x):

Vfermi := proc (r) options operator, arrow; -4*Pi*n0*R^2*(R*(v1(r/R)-v10)/r+v2Rm-v2(r/R)) end proc:

Vuniform := proc (r) options operator, arrow; piecewise(r < R, -Z*(3/2-(1/2)*r^2/R^2)/R, -Z/r) end proc:

plot([Vuniform(r), Vfermi(r), -Z/r], r = r0 .. 2*R, V = 1.2*Vfermi(r0) .. 0, legend = ["uniform charge", "Fermi distribution", "point charge"]);



plotsol1s := proc (E, K, r0, S, col) local Eqns, ICs, fnl, gnl, r, soln; global ms, c; Eqns := diff(fnl(r), r) = gnl(r)*[E+2*ms*c^2-Vuniform(r)]/c-(K+1)*fnl(r)/r, diff(gnl(r), r) = -fnl(r)*[E-Vuniform(r)]/c-(1-K)*fnl(r)/r; ICs := fnl(r0) = 1, gnl(r0) = 0; soln := dsolve({Eqns, ICs}, numeric); plots:-odeplot(soln, [r, fnl(r)], r0 .. S, color = col) end proc:

plotsol1s(-3*10^5, -1, 10^(-10), Rm, red)

Error, (in f) unable to store '[HFloat(0.0)]' when datatype=float[8]




Any help would be greatly appreciated.


Gambia Man


sys := [(diff(c(t), t))*(diff(a(t), t))*(diff(b(t), t))+(diff(c(t), t))*t^2, (diff(c(t), t))*(diff(a(t), t))*t+(diff(c(t), t))*t*(diff(b(t), t)), (diff(c(t), t))*(diff(a(t), t))*(diff(b(t), t))+(diff(c(t), t))*(diff(a(t), t))*t+(diff(c(t), t))*t^2]
DEplot(sys, [a(t), b(t), c(t)], t = 0 .. 2, a = -15 .. 15, b = -15 .. 15, c = -15 .. 15, color = magnitude, title = `Stable Limit Cycles`, arrows = curve, dirfield = 800, axes = none);
odeplot(sys, [t, a(t), b(t), c(t)], -4 .. 4, color = orange);

local eqs,yp2,yp4;
eqs:=[yp2*Y[3]+yp4*Y[2]*sin(Y[1]^2)+cos(yp4*Y[3]) = sin(t), Y[2]*yp4*sin(Y[1]*Y[3])+5*yp2*Y[4]*cos(Y[1]^2)+t^2*Y[1]*Y[3]^2 = exp(-Y[3]^2)];
end proc:

It's a long time to wait for the odeplot.

Any advice is appreciated.

1 2 3 Page 1 of 3