Items tagged with ode ode Tagged Items Feed

Hi, Im now trying to run my code. But it took like years to even getting the results. may I know any solutions on how to get faster results? Because I have run this code for almost 4 hours yet there is still 'Evaluating...' at the corner left. And when I tried to stop the program, it will stop at 'R1...'.

 

Digits := 18;
with(plots):n:=1.4: mu(eta):=(diff(U(eta),eta)^(2)+diff(V(eta),eta)^(2))^((n-1)/(2)):
Eqn1 := 2*U(eta)+(1-n)*eta*(diff(U(eta), eta))/(n+1)+diff(W(eta), eta) = 0;
Eqn2 := U(eta)^2-(V(eta)+1)^2+(W(eta)+(1-n)*eta*U(eta)/(n+1))*(diff(U(eta), eta))-mu(eta)*(diff(U(eta), eta, eta))-(diff(U(eta), eta))*(diff(mu(eta), eta)) = 0;
Eqn3 := 2*U(eta)*(V(eta)+1)+(W(eta)+(1-n)*eta*U(eta)/(n+1))*(diff(V(eta), eta))-mu(eta)*(diff(V(eta), eta, eta))-(diff(V(eta), eta))*(diff(mu(eta), eta)) = 0;
bcs1 := U(0) = 0, V(0) = 0, W(0) = 0;
bcs2 := U(4) = 0, V(4) = -1;
R1 := dsolve({Eqn1, Eqn2, Eqn3, bcs1, bcs2}, {U(eta), V(eta), W(eta)}, initmesh = 30000, output = listprocedure, numeric);
Warning, computation interrupted
for l from 0 by 2 to 4 do R1(l) end do;
plot1 := odeplot(R1, [eta, U(eta)], 0 .. 4, numpoints = 2000, color = red);

 

Thankyou in advance :)

how can i improve doing the graph with various parameter

do anyone have abother numerical method in maple rather than rk45 felhberg

like keller box/homotopy/ or anything

i have attchassignment.mwsassignment.pdf

restart

with(plots)

%?

Eq1 := diff(f(eta), `$`(eta, 3))+(diff(f(eta), `$`(eta, 2)))*f(eta)-(diff(f(eta), eta))^2+4 = 0

diff(diff(diff(f(eta), eta), eta), eta)+(diff(diff(f(eta), eta), eta))*f(eta)-(diff(f(eta), eta))^2+4 = 0

(1)

%?

Eq2 := diff(theta(eta), `$`(eta, 2))+Pr*(diff(theta(eta), eta))*f(eta) = 0

diff(diff(theta(eta), eta), eta)+Pr*(diff(theta(eta), eta))*f(eta) = 0

(2)

%?

VPr := [0.1e-1, 0.2e-1, 0.3e-1]

etainf := 27

bcs := (D(f))(0) = 0, f(0) = 0, (D(theta))(0) = -1, (D(f))(etainf) = 2, theta(etainf) = 0

(D(f))(0) = 0, f(0) = 0, (D(theta))(0) = -1, (D(f))(27) = 2, theta(27) = 0

(3)

dsys := {Eq1, Eq2, bcs}

for i to 3 do Pr := VPr[i]; dsol[i] := dsolve(dsys, numeric); print(Pr); print(dsol[i](0)) end do

0.1e-1

[eta = 0., f(eta) = HFloat(0.0), diff(f(eta), eta) = HFloat(0.0), diff(diff(f(eta), eta), eta) = HFloat(3.4862842650940435), theta(eta) = HFloat(9.305856096452466), diff(theta(eta), eta) = HFloat(-0.9999999999999998)]

0.2e-1

[eta = 0., f(eta) = HFloat(0.0), diff(f(eta), eta) = HFloat(0.0), diff(diff(f(eta), eta), eta) = HFloat(3.4862842653392216), theta(eta) = HFloat(6.7064688361073745), diff(theta(eta), eta) = HFloat(-1.0)]

0.3e-1

[eta = 0., f(eta) = HFloat(0.0), diff(f(eta), eta) = HFloat(0.0), diff(diff(f(eta), eta), eta) = HFloat(3.4862842648459234), theta(eta) = HFloat(5.552583608770695), diff(theta(eta), eta) = HFloat(-0.9999999999999998)]

(4)

SDf1 := odeplot(dsol[1], [eta, f(eta)], 0 .. etainf, color = green, axes = box); SDf2 := odeplot(dsol[2], [eta, f(eta)], 0 .. etainf, color = red); SDf3 := odeplot(dsol[3], [eta, f(eta)], 0 .. etainf, color = blue)

display([SDf1, SDf2, SDf3], labels = ["η", "f (η)"], labeldirections = [horizontal, vertical], labelfont = [italic, 16, bold], axes = boxed, axesfont = [times, 14], thickness = 3)
%?

 

%?

SDfd1 := odeplot(dsol[1], [eta, diff(f(eta), eta)], 0 .. etainf, color = green, axes = box); SDfd2 := odeplot(dsol[2], [eta, diff(f(eta), eta)], 0 .. etainf, color = red); SDfd3 := odeplot(dsol[3], [eta, diff(f(eta), eta)], 0 .. etainf, color = blue)

%?

display([SDfd1, SDfd2, SDfd3], labels = ["η", "f  ' (η)"], labeldirections = [horizontal, vertical], labelfont = [italic, 16, bold], axes = boxed, axesfont = [times, 14], thickness = 3)

 

%?

`Sθ1` := odeplot(dsol[1], [eta, theta(eta)], 0 .. etainf, color = green, axes = box); `Sθ2` := odeplot(dsol[2], [eta, theta(eta)], 0 .. etainf, color = red); `Sθ3` := odeplot(dsol[3], [eta, theta(eta)], 0 .. etainf, color = black)

display([`Sθ1`, `Sθ2`, `Sθ3`], labels = ["η", "θ (η)"], labeldirections = [horizontal, vertical], labelfont = [italic, 16, bold], axes = boxed, axesfont = [times, 14], thickness = 3)

 

%?

`Sθd1` := odeplot(dsol[1], [eta, diff(theta(eta), eta)], 0 .. etainf, color = green, axes = box); `Sθd2` := odeplot(dsol[2], [eta, diff(theta(eta), eta)], 0 .. etainf, color = red); `Sθd3` := odeplot(dsol[3], [eta, diff(theta(eta), eta)], 0 .. etainf, color = black)

display([`Sθd1`, `Sθd2`, `Sθd3`], labels = ["η", "θ '(η)"], labeldirections = [horizontal, vertical], labelfont = [italic, 16, bold], axes = boxed, axesfont = [times, 14], thickness = 3)

 

%?

 

 

Download assignment.mws

 

Hello Everybody,

I was trying to apply the Newton-Raphon method("Newton" in maple) for the following problem to obtain the constants (c1...cM) without success.

eq[1] := -0.0139687 c[2] - 0.0132951 c[1] = 0
eq[2] := 24806.4 c[2] - 0.0139687 c[1] = 0

If anybody has experience with it or knows how to use it I would really appreciate your help.

 

Just to introduce the previous steps of calculations leading to the problem:To calculate the critical buckling force N and the shape of a rectangular uniformly loaded plate the governing diff. equation is the following

D11*w''''(x)+(2*D12+4*D66)*w''(x)''(y)+D22*w''''(y)+N*w''(x)=0

For the follwing solution of a boundary value problem(Boundary conditions:clamped-clamped: w(x=0,x=a)=0 & w'(x=0,x=a)=0) I applied the Ritz method:

w:=sum(c[i]*(cos(2*i*Pi*x/a)-1)),i=1..M)

Thus the Potential Energy P is:

P:=16537.6*c[2]^2-0.0139687*c[2]*c[1]+1033,60*c[1]^2-9.86960*N*c[2]^2-2.46740*N*c[1]^2

 deriving P to each constant c and setting them =0 leads to:

eq[1] := -0.0139687 c[2] + 2067.20 c[1] - 4.93480 N c[1] = 0
eq[2] := 33075.2 c[2] - 0.0139687 c[1] - 19.7392 N c[2] = 0

After calculating N=418,902  and feeding eq1 and eq 2,the follwing equations if two terms are considered:

eq[1] := -0.0139687 c[2] - 0.0132951 c[1] = 0
eq[2] := 24806.4 c[2] - 0.0139687 c[1] = 0

Everything I tried resulted in any c1 = c2 = 0 which is not  realistic. Maybe I made a mistake earlier.

Thanks a lot in advance.

Sam

 Plate_Buckling.mw

I have an ODE with L(x,y), which includes partial derivatives of L(x,y) and coefficients y,_y1 and _y1^2. I want to split this equation into seperate equations of these coefficients. So I can solve for the unknown variable L(x,y).

 

ODE:=diff(L(x,y),x,x)+2*_y1(x)*diff(L(x,y),x,y)+_y1(x)^2*diff(L(x,y),y,y)+y(x)/x*diff(L(x,y),x)+(y(x)*_y1(x)/x+_y1(x))*diff(L(x,y),y)-y(x)/x^2*L(x,y)=0;

Could someone help me? 

 

I am currently trying to use the coeffs command but not really getting anywhere.

 

Thank you.

I'm trying to find lypunov exponent for this  system of ODEs. I know I need to take the Jacobian but not sure if it's possible the way it's currently defined. If anyone could provide insight it would be much appreciated.

Eqns:= 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(theta(t),t)=omega(t), diff(xH(t),t)=vxH(t),diff(vxH(t),t)=-(G*M*xH(t))/((xH(t)^(2)+yH(t)^(2))^(1.5)),diff(yH(t),t)=vyH(t),diff(vyH(t),t)=-(G*M*yH(t))/((xH(t)^(2)+yH(t)^(2))^(1.5)): ;

ICs := omega(0) = omega0, theta(0) = theta0, xH(0) = a*(1+e), yH(0) = 0, vxH(0) = 0, vyH(0) = sqrt(G*M*(1-e)/(a*(1+e)));

I want the exponent for omega, I  procedure that takes some initial conditions, changing just w0, and computes the long term value of omega. This plots a sort of bifurcation diagram. I'd like an estimate of the exponent to compare what I see. 

thanks for the help

ft

I don't know what is the problem with this. Please take a look at it. Any help is appreciated. Thanks.

restart

with(RealDomain)

[Im, Re, `^`, arccos, arccosh, arccot, arccoth, arccsc, arccsch, arcsec, arcsech, arcsin, arcsinh, arctan, arctanh, cos, cosh, cot, coth, csc, csch, eval, exp, expand, limit, ln, log, sec, sech, signum, simplify, sin, sinh, solve, sqrt, surd, tan, tanh]

(1)

with(Slode)

[DEdetermine, FPseries, FTseries, candidate_mpoints, candidate_points, dAlembertian_formal_sol, dAlembertian_series_sol, hypergeom_formal_sol, hypergeom_series_sol, mhypergeom_formal_sol, mhypergeom_series_sol, msparse_series_sol, polynomial_series_sol, rational_series_sol, series_by_leastsquare]

(2)

ODE := x*(diff(y(x), x, x))-(diff(y(x), x))-x*y(x)

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

(3)

SS := convert(sin(x), FormalPowerSeries)

Sum((-1)^k*x^(2*k+1)/factorial(2*k+1), k = 0 .. infinity)

(4)

FPseries(ODE = SS, y(x), v(n))

Error, the right-hand side of the differential equation must be a polynomial or a series in x

 

``

 

Download AA.mw

Dear all;

Good morning everyone.

I solve a simple ode, i want how can I write this program as procedure with output the two coefficient involved in the solution, after solving with ics. i.e coef:=dsolve(ode);  in my example I want as output [4, -1] .

 

restart:
ode:=diff(y(x),x,x)=3*diff(y(x),x)-2*y(x);
coef:=dsolve(ode);
ics:=y(0)=3, D(y)(0)=2;
dsolve({ics,ode});

with best regards

 

Hello,

I have two regimes. Each regime is characterized by system of state differential equations (diff(S(t), t), diff(K(t), t)) and co-state differential equations (diff(psi[S](t), t), diff(psi[Iota](t), t)) as follows: 

(1)

diff(S(t), t) = -eta*K(t)*S(t)/(w*N*(S(t)+K(t))), diff(K(t), t) = eta*K(t)*S(t)/(w*N*(S(t)+K(t)))-upsilon

diff(psi[S](t), t) = eta*K(t)^2*(psi[S](t)-psi[Iota](t))/(w*N*(S(t)+K(t))^2), diff(psi[Iota](t), t) = eta*S(t)^2*(psi[S](t)-psi[Iota](t))/(w*N*(S(t)+K(t))^2)

(2)

diff(S(t), t) = -eta*K(t)*S(t)/(w*N*(S(t)+K(t))), diff(K(t), t) = eta*K(t)*S(t)/(w*N*(S(t)+K(t)))

diff(psi[S](t), t) = eta*K(t)^2*(psi[S](t)-psi[Iota](t))/(w*N*(S(t)+K(t))^2), diff(psi[Iota](t), t) = eta*S(t)^2*(psi[S](t)-psi[Iota](t))/(w*N*(S(t)+K(t))^2)

The first regime is employed from 0 to t1 (where t1 is unknown) and then regime 2, from t1 to T. I know the initial values for the state variables of the first system at t=0, that is, S(0)=S0 and K(0)=K0, as well as boundary conditions for the co-state variables for regime 2 at t=T, that is, psi_S(T)=a and psi_I(T)=b. 

I also know that my unknown t1 should satisfy the algebraic equation: -c - psi_S(t1)=0, where psi_S(t1) is the solution of the co-state diff equation of the regime 2 at t=t1. 

My question is: If I assume that the systems have not analytical solution, how can I found unknown t1 numerically? Moreover, asssume that all other parameters, such as T, eta, upsilon and others are given.

Hello,

My code records the values I need, however, I need to implement a modulo of 2*Pi on my result for theta. But this leads to a graph with no plots and I'm not sure how to fix it. Any help is greatly aprreciated! Thank you in advance!

Kind regards,

Gam

with(plots):

a := 1.501*10^9:

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

HyperionOrbit := proc (`θIC`, `ωIC`, n) local a, Mh, Msat, G, e, beta, M, Eqns, ICs, soln; option remember; global `ω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) = `θIC`, omega(0) = `ωIC`; soln := dsolve({Eqns, ICs}, numeric, maxfun = 0, output = array([seq(i, i = 0 .. n*Th, Th)])); plots:-odeplot(soln, [modp(theta(t), 2*Pi), omega(t)/`ωH`], 0 .. n*Th, labels = ["θ(t)","ω(t)/ωH"], axes = boxed, style = plottools:-point, size = [.25, .75]) end proc:

plots:-display(HyperionOrbit(.5, 1.8*`ωH`, 10));

Download Poincare_section_Boyd_plot_fixing_theta.mw

bia Man

Hello,

I have a question about poincare sections. I have this piece of code i need to analyse and I want to use a poincare section in order to so. How could I do it? I am interested in theta and omega. Any help is greatly appreciated! Thank you in advance!

Kind regards,

Gambia Man

with(plots):

a := 1.501*10^9:

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

1876321.326

 

0.3348672330e-5

(1)

HyperionOrbit := proc (`θIC`, `ωIC`) local a, Mh, Msat, G, e, beta, M, Eqns, ICs; global `ωH`, Th, soln; 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) = `θIC`, omega(0) = `ωIC`; soln := dsolve({Eqns, ICs}, numeric); odeplot(soln, [theta(t), omega(t)/`ωH`], 0 .. 5*Th, numpoints = 2000, labels = ["θ(t)","ω(t)/ωH"], axes = boxed, size = [.25, .75]) end proc

``

 

Download New_Poincare_section.mw

Hi everyone, I'm working a problem

Given a dog and a man. At t=0, the position of the man and his dog is (0,0) and (0,h), respectively. Then, the man start moving along Ox with constant speed vm. The dog keep running toward its master with constant speed vd. Describe the movement of the dog?

Therefore, let say the position of the dog is (f1(x),f2(x)), I wrote these:

restart;

with(plots);

MD := proc (h, vm, vd) local x, y, ode, ics, sol;

ode := {(diff(f2(t), t))*(vm*t-f1(t))+(diff(f1(t), t))*f2(t), diff(f2(t), t) = -sqrt(vd^2-(diff(f1(t), t))^2)};

ics := {f1(0) = h, f2(0) = 0};

sol := dsolve(ode union ics);

x := unapply(eval(f1(t), sol), t);

y := unapply(eval(f2(t), sol), t);

plot([x(t), y(t), t = 0 .. 10], scaling = constrained) end proc;

MD(10, 1, 5);

However, "sol" is returned NULL, which means the equation has no solution. I supposed I completed the motion part of the problem correctly. Please help me or point out an another method

 

Hello,

I have an error in my code. I don't know where it came from. Earlier today I loaded this and it worked fine and now an error comes up. Any help is greatly appreciated! Thank you in advance!

Kind regards,

Gambia Man

with(plots):

``

HyperionOrbit := proc (`θIC`, `ωIC`) local a, Mh, Msat, G, e, beta, M, Eqns, Th, ICs, soln, `ωH`; 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); Th := sqrt(4*Pi^2*a^3/(G*(Mh+Msat))); `ωH` := 2*Pi/Th; ICs := xH(0) = a*(1+e), yH(0) = 0, vxH(0) = 0, vyH(0) = sqrt(G*M*(1-e)/(a*(1+e))), theta(0) = `θIC`, omega(0) = `ωIC`; soln := dsolve({Eqns, ICs}, numeric); odeplot(soln, [theta(t), omega(t)/`ωH`], 0 .. 5*Th, numpoints = 2000, labels = ["θ(t)","ω(t)/ωH"], axes = boxed) end proc

HyperionOrbit(.5, 1.8*`ωH`)

Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)

 

Error, (in dsolve/numeric/SC/IVPpoints) unable to store '11.3097335529232/Th' when datatype=float[8]

 

``

``

``

 

Download Poincare_section.mw

Hello,

I have an issue with this code. The issue lies in the warning which is outputed. However, I don't know how to resolve it. Any help would be greatly appreciated! Thank you in advance!

Kind regards,

Gambia Man

with(plots):

a := 1.501*10^9:

Eqns := diff(xH(t), t) = vxH(t), diff(vxH(t), t) = -G*(Mh+Msat)*xH(t)/(xH(t)^2+yH(t)^2)^(3/2), diff(yH(t), t) = vyH(t), diff(vyH, t) = -G*(Mh+Msat)*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*(Mh+Msat)*(1-e)/(a*(1+e))):

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

soln := dsolve({Eqns, ICs}, numeric):

plots:-odeplot(soln, [xH(t)/a, yH(t)/a], 0 .. 10*T, scaling = constrained, labels = ["x/a", "y/a"], numpoints = 2000);

Warning, cannot evaluate the solution further right of 453574.15, probably a singularity

 

 

``

 

Download Hyperion_Orbit.mw

 

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!

 

This is the code hw2_final.mw

This is the warning 

Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters )"

 

How to solve it?

3 4 5 6 7 8 9 Last Page 5 of 29