Items tagged with differential_equation

I would like to solve the following differential equation with a relatively complicated function that is best declared.

myfun := proc (x::float)

local output;

output := 4*x^2;

end proc;


de := diff(y(x), x)+myfun(x)*y(x) = 0.;

However, this gives the following message:

Error, invalid input: myfun expects its 1st argument, x, to be of type float, but received x

Any suggestions?

 

 

I wonder how to solve symbolically for like following coupled ODEs in Maple? 

On the other hand, I want to write the code for step by step solution of this problem. But I didn't find algorithm of the solution on the some books. Do you know some books including solving coupled ODEs ?

 

How to solve delay differential equation by method of steps in MAPLE software. 

I have been trying to find a solution for the equation below. Is there a non numerical explicit solution?
 

restart

with(DEtools)

[AreSimilar, Closure, DEnormal, DEplot, DEplot3d, DEplot_polygon, DFactor, DFactorLCLM, DFactorsols, Dchangevar, Desingularize, FunctionDecomposition, GCRD, Gosper, Heunsols, Homomorphisms, IVPsol, IsHyperexponential, LCLM, MeijerGsols, MultiplicativeDecomposition, ODEInvariants, PDEchangecoords, PolynomialNormalForm, RationalCanonicalForm, ReduceHyperexp, RiemannPsols, Xchange, Xcommutator, Xgauge, Zeilberger, abelsol, adjoint, autonomous, bernoullisol, buildsol, buildsym, canoni, caseplot, casesplit, checkrank, chinisol, clairautsol, constcoeffsols, convertAlg, convertsys, dalembertsol, dcoeffs, de2diffop, dfieldplot, diff_table, diffop2de, dperiodic_sols, dpolyform, dsubs, eigenring, endomorphism_charpoly, equinv, eta_k, eulersols, exactsol, expsols, exterior_power, firint, firtest, formal_sol, gen_exp, generate_ic, genhomosol, gensys, hamilton_eqs, hypergeomsols, hyperode, indicialeq, infgen, initialdata, integrate_sols, intfactor, invariants, kovacicsols, leftdivision, liesol, line_int, linearsol, matrixDE, matrix_riccati, maxdimsystems, moser_reduce, muchange, mult, mutest, newton_polygon, normalG2, ode_int_y, ode_y1, odeadvisor, odepde, parametricsol, particularsol, phaseportrait, poincare, polysols, power_equivalent, rational_equivalent, ratsols, redode, reduceOrder, reduce_order, regular_parts, regularsp, remove_RootOf, riccati_system, riccatisol, rifread, rifsimp, rightdivision, rtaylor, separablesol, singularities, solve_group, super_reduce, symgen, symmetric_power, symmetric_product, symtest, transinv, translate, untranslate, varparam, zoom]

(1)

``

sol := (JacobiCN((1/10)*sqrt(5)*sqrt(2)*t, (1/3)*sqrt(3))*sqrt(5)+2*sqrt(2))/(JacobiCN((1/10)*sqrt(5)*sqrt(2)*t, (1/3)*sqrt(3))*sqrt(5)+5*sqrt(2))

(JacobiCN((1/10)*5^(1/2)*2^(1/2)*t, (1/3)*3^(1/2))*5^(1/2)+2*2^(1/2))/(JacobiCN((1/10)*5^(1/2)*2^(1/2)*t, (1/3)*3^(1/2))*5^(1/2)+5*2^(1/2))

(2)

sol1 := diff(psi(t), t) = sol

diff(psi(t), t) = (JacobiCN((1/10)*5^(1/2)*2^(1/2)*t, (1/3)*3^(1/2))*5^(1/2)+2*2^(1/2))/(JacobiCN((1/10)*5^(1/2)*2^(1/2)*t, (1/3)*3^(1/2))*5^(1/2)+5*2^(1/2))

(3)

odeadvisor(sol1, psi(t))

[_quadrature]

(4)

sol2 := dsolve({sol1, psi(0) = 0})

psi(t) = Int((JacobiCN((1/10)*_z1*10^(1/2), (1/3)*3^(1/2))*5^(1/2)+2*2^(1/2))/(JacobiCN((1/10)*_z1*10^(1/2), (1/3)*3^(1/2))*5^(1/2)+5*2^(1/2)), _z1 = 0 .. t)

(5)

``

``

sol3 := convert(sol, parfrac)

1-(3/5)*2^(1/2)*5^(1/2)/(JacobiCN((1/10)*5^(1/2)*2^(1/2)*t, (1/3)*3^(1/2))+2^(1/2)*5^(1/2))

(6)

sol4 := diff(psi(t), t) = sol3

diff(psi(t), t) = 1-(3/5)*2^(1/2)*5^(1/2)/(JacobiCN((1/10)*5^(1/2)*2^(1/2)*t, (1/3)*3^(1/2))+2^(1/2)*5^(1/2))

(7)

dsolve({sol4, psi(0) = 0})

psi(t) = Int(-(3/5)*10^(1/2)/(JacobiCN((1/10)*_z1*10^(1/2), (1/3)*3^(1/2))+10^(1/2)), _z1 = 0 .. t)+t

(8)

sol5 := diff(psi(t), t) = 3*sqrt(2)*sqrt(5)/(5*(JacobiCN((1/10)*sqrt(5)*sqrt(2)*t, (1/3)*sqrt(3))+sqrt(5)*sqrt(2)))

diff(psi(t), t) = 3*2^(1/2)*5^(1/2)/(5*2^(1/2)*5^(1/2)+5*JacobiCN((1/10)*5^(1/2)*2^(1/2)*t, (1/3)*3^(1/2)))

(9)

odeadvisor(sol5)

[_quadrature]

(10)

dsolve({sol5})

{psi(t) = Int((3/5)*10^(1/2)/(JacobiCN((1/10)*t*10^(1/2), (1/3)*3^(1/2))+10^(1/2)), t)+_C1}

(11)

``


 

Download Jacobi_Diff_eqn_Mapleprime_post.mw

restart;

Digits := 18;
with(LinearAlgebra);
f := proc (n) 3*sin(x[n]) end proc;

g := proc (n) 3*cos(x[n])

end proc;

#problem call.
for n from 0 to 0 do

e1 := expand(-y[n+3/2]+y[n]-3*y[n+1/2]+3*y[n+1]+1/11612160*(5856*h^4*g(n+1/2)-19968*h^4*g(n+3/2)+2343*h^4*g(n)-76356*h^4*g(n+1)-7058*h^4*g(n+2)+608864*h^3*f(n+1/2)+104864*h^3*f(n+3/2)+28489*h^3*f(n)+702864*h^3*f(n+1)+6439*h^3*f(n+2)));

e2 := expand(-y[n+2]+3*y[n]-8*y[n+1/2]+6*y[n+1]+1/5806080*(18768*h^4*g(n+1/2)-32880*h^4*g(n+3/2)+3867*h^4*g(n)-76356*h^4*g(n+1)-2229*h^4*g(n+2)+965728*h^3*f(n+1/2)+461728*h^3*f(n+3/2)+45953*h^3*f(n)+1405728*h^3*f(n+1)+23903*h^3*f(n+2)));

e3 := expand(-z[n]+(1/383201280*(-4207440*h^4*g(n+1/2)-930192*h^4*g(n+3/2)+371973*h^4*g(n)-3631932*h^4*g(n+1)-41259*h^4*g(n+2)+16136096*h^3*f(n+1/2)+3866720*h^3*f(n+3/2)+5752543*h^3*f(n)+5810400*h^3*f(n+1)+367681*h^3*f(n+2))+4*y[n+1/2]-3*y[n]+y[n+1])/h);

e4 := expand(-z[n+1/2]+(1/191600640*(376320*h^4*g(n+1/2)+118896*h^4*g(n+3/2)-29469*h^4*g(n)+532764*h^4*g(n+1)+5079*h^4*g(n+2)-5812112*h^3*f(n+1/2)-508016*h^3*f(n+3/2)-381553*h^3*f(n)-1236168*h^3*f(n+1)-45511*h^3*f(n+2))-y[n]+y[n+1])/h);

e5 := expand(-z[n+1]+(1/383201280*(-31920*h^4*g(n+1/2)-433776*h^4*g(n+3/2)+71547*h^4*g(n)-2519748*h^4*g(n+1)-17493*h^4*g(n+2)+18565216*h^3*f(n+1/2)+1933216*h^3*f(n+3/2)+885665*h^3*f(n)+10391328*h^3*f(n+1)+158015*h^3*f(n+2))-5*y[n+1/2]+y[n]+3*y[n+1])/h);

e6 := expand(-z[n+3/2]+(1/95800320*(250224*h^4*g(n+1/2)-730680*h^4*g(n+3/2)+61266*h^4*g(n)-1526256*h^4*g(n+1)-22044*h^4*g(n+2)+15680504*h^3*f(n+1/2)+4712456*h^3*f(n+3/2)+735469*h^3*f(n)+22576428*h^3*f(n+1)+203623*h^3*f(n+2))-8*y[n+1/2]+3*y[n]+5*y[n+1])/h);

e7 := expand(-z[n+2]+(1/383201280*(3873264*h^4*g(n+1/2)+332976*h^4*g(n+3/2)+497649*h^4*g(n)-1407564*h^4*g(n+1)-720255*h^4*g(n+2)+114710816*h^3*f(n+1/2)+93716192*h^3*f(n+3/2)+5705827*h^3*f(n)+191366496*h^3*f(n+1)+9635389*h^3*f(n+2))-12*y[n+1/2]+5*y[n]+7*y[n+1])/h);

e8 := expand(-p[n]+(1/191600640*(13423440*h^4*g(n+1/2)+3068304*h^4*g(n+3/2)-1621317*h^4*g(n)+11615292*h^4*g(n+1)+137451*h^4*g(n+2)-32503712*h^3*f(n+1/2)-12664928*h^3*f(n+3/2)-32539039*h^3*f(n)-16869600*h^3*f(n+1)-1223041*h^3*f(n+2))-8*y[n+1/2]+4*y[n]+4*y[n+1])/h^2);

e9 := expand(-p[n+1/2]+(1/191600640*(-3053856*h^4*g(n+1/2)-213216*h^4*g(n+3/2)+98049*h^4*g(n)-509436*h^4*g(n+1)-10191*h^4*g(n+2)-1045120*h^3*f(n+1/2)+831104*h^3*f(n+3/2)+1331083*h^3*f(n)-1207008*h^3*f(n+1)+89941*h^3*f(n+2))-8*y[n+1/2]+4*y[n]+4*y[n+1])/h^2);

e10 := expand(-p[n+1]+(1/63866880*(194160*h^4*g(n+1/2)-373968*h^4*g(n+3/2)+52329*h^4*g(n)-2514924*h^4*g(n+1)-14727*h^4*g(n+2)+14006304*h^3*f(n+1/2)+1695712*h^3*f(n+3/2)+634955*h^3*f(n)+15463008*h^3*f(n+1)+133461*h^3*f(n+2))-8*y[n+1/2]+4*y[n]+4*y[n+1])/h^2);

e11 := expand(-p[n+3/2]+(1/191600640*(1491168*h^4*g(n+1/2)-4758240*h^4*g(n+3/2)+190977*h^4*g(n)-509436*h^4*g(n+1)-103119*h^4*g(n+2)+46274944*h^3*f(n+1/2)+48151168*h^3*f(n+3/2)+2215307*h^3*f(n)+93985056*h^3*f(n+1)+974165*h^3*f(n+2))-8*y[n+1/2]+4*y[n]+4*y[n+1])/h^2);

e12 := expand(-p[n+2]+(1/191600640*(4772688*h^4*g(n+1/2)+11719056*h^4*g(n+3/2)+338619*h^4*g(n)+11615292*h^4*g(n+1)-1822485*h^4*g(n+2)+59770976*h^3*f(n+1/2)+79609760*h^3*f(n+3/2)+3528289*h^3*f(n)+109647648*h^3*f(n+1)+34844287*h^3*f(n+2))-8*y[n+1/2]+4*y[n]+4*y[n+1])/h^2) end do;
M := {e || (1 .. 12)};

y_init := 1;

z_init := 0;

p_init := -2;

x_init := 0; A := 0; B := 1; N := 40;

h := evalf((B-A)/N); count := 1;

X := y[k], y[k+1/2], y[k+1], y[k+3/2], z[k], z[k+1/2], z[k+1], z[k+3/2], p[k], p[k+1/2], p[k+1], p[k+3/2];

step := seq(eval(x, x = n*h), n = 1 .. N);

y_exact := ([seq])(eval(3*cos(x)+(1/2)*x^2-2, x = n*h), n = 1 .. N);

z_exact := ([seq])(eval((1/3*(3*x^2+6*x+3))/(x^3+3*x^2+3*x+1), x = n*h), n = 1 .. N);

p_exact := ([seq])(eval((1/3*(6*x+6))/(x^3+3*x^2+3*x+1)-(1/3)*(3*x^2+6*x+3)^2/(x^3+3*x^2+3*x+1)^2, x = n*h), n = 1 .. N);
vars := seq(X, k = 1);
printf("\n%4s%13s%15s%15s\n", "@", "y_Num", "y_Exact", "y_Error");

for q to N do

for ix to 4 do

x[ix] := h*ix+x_init end do;

result := eval(`<,>`(vars), fsolve(eval(M, [x[0] = x_init, x[1/2] = x_init, x[3/2] = x_init, y[0] = y_init, y[1/2] = y_init, y[3/2] = y_init, z[0] = z_init, z[1/2] = z_init, z[3/2] = z_init, p[0] = p_init, p[1/2] = p_init, p[3/2] = p_init]), {vars}));

for k to 4 do

printf("%5.2f %14.15f", step[count], result[k]);

printf("%20.15f %10.18G \n", y_exact[count], abs(result[k]-y_exact[count]));

count := count+1;

P := [result[k]]

end do;

x_init := x[ix-1];

y_init := result[4];

z_init := result[8];

p_init := result[12]

end do;

 

please that is the code i write to solve the problem after using the matrix form to generate the value but is given me error of the form


   @        y_Num        y_Exact        y_Error
Error, invalid input: eval received fsolve({-6398.00004614630940+6400.00000000000000*y[1], -6397.99992849910140+6400.00000000000000*y[1], -6397.99909739580050+6400.00000000000000*y[1], -199.999989717789185+200.000000000000000*y[1], -40.0000000791700798+40.0000000000000000*y[1], -2.99999993737911015+3*y[1], 39.999999768462113+40.0000000000000000*y[1], -p[1]-6399.99961623646730+6400.00000000000000*y[1], -p[2]-6399.99798489466010+6400.00000000000000*y[1], -y[2]-4.99999972458202552+6*y[1], -z[1]-159.999999048856193+120.000000000000000*y[1], -z[2]-279.999973921987948+280.000000000000000*y[1]}, {p[1], p[2], p[3/2], p[5/2], y[1], y[2], y[3/2], y[5/...
 

This question explores the family of differential equations dy/dx=sqrt(􏰐 1 +􏰏( a*x )+ 􏰏 (2 *y)) for various values of the parameter a.  

For the case a = 􏰐 0 find the analytical solution that passes through the point (0, 1) and verify that this is a solution to the differential equation. Use this solution to find the value of y correct to 4 decimal placeswhen x=􏰐1. 

In maple i did

y:=(1/2)*x^2+sqrt(3)*x+1:
diff(y,x)
                             
i got the answer x + sqrt(3)

as shown in the markscheme. please cluld anyone help how to get y before this step and what to do after.

    

 

 

I'm new using maple and I trying to solve this equation:

u*(diff(u*(diff(R(u), u)), u))-(-m^2+u^2+fu)*R(u) = 0

Maple give this solution:

R(u) = _C1*BesselI(sqrt(-m^2+fu), u)+_C2*BesselK(sqrt(-m^2+fu), u)

But I don't know how I can find de constant C1 and C2.

Thanks in advance.

 

Dear All,

I am working on ODEs and have obtained the plot for "variable vs time". I would like to know if it is possible and how to analyze those data in the frequency domain.ODEs.mw

Thank you.

Very kind wishes,

Wang Zhe

Hi guys,

I'm trying to solve this system I have but the solution doesn't display:

-I have two second degree differential equations with two functions.

-I have a set of two boundary conditions per function.

 

Thank you!
 

``

restart;

eq2:=2*diff(y(x), x$2)+diff(z(x), x$2)=0;

2*(diff(diff(y(x), x), x))+diff(diff(z(x), x), x) = 0

(1)

eq3:=2*diff(z(x), x$2)+diff(y(x), x$2)=0;

2*(diff(diff(z(x), x), x))+diff(diff(y(x), x), x) = 0

(2)

SOL:=dsolve({eq2, eq3, D(y)(0)=0, D(y)(1)=1, D(z)(0)=0, D(z)(1)=1}, {y(x), z(x)});

"SOL := "

(3)

``


 

Download trial.mwtrial.mw

Hello,

I have been trying to solve a simple nonlinear equation. Im interested in the solution per say rather than the plot but when I browsed about the commands to use, this came up. I tried it in my case and it is giving me the following errors:

ode.mw
 

``

restart;

``

with(plots);

[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fieldplot, fieldplot3d, gradplot, gradplot3d, implicitplot, implicitplot3d, inequal, interactive, interactiveparams, intersectplot, listcontplot, listcontplot3d, listdensityplot, listplot, listplot3d, loglogplot, logplot, matrixplot, multiple, odeplot, pareto, plotcompare, pointplot, pointplot3d, polarplot, polygonplot, polygonplot3d, polyhedra_supported, polyhedraplot, rootlocus, semilogplot, setcolors, setoptions, setoptions3d, shadebetween, spacecurve, sparsematrixplot, surfdata, textplot, textplot3d, tubeplot]

(1)

eq5:=C*sqrt(y(x)*((diff(y(x),x))^2+1))-y(x)=0;

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

(2)

C:=1;

1

(3)

bcs:=y(-1)=1, y(1)=1;

y(-1) = 1, y(1) = 1

(4)

dsys:={eq5,bcs};

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

(5)

dsol:=dsolve(dsys, numeric); odeplot(dsol,[x,y(x)],0..1,color=red,axes=box);

Error, (in dsolve/numeric/bvp/convertsys) unable to convert to an explicit first-order system

 

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

 

``


 

Download ode.mw

 

A family of curves has polar equation r=cos^n (theta/n), 0<=theta,n*pi, where n is a positive even integer.

Using t = theta as the parameter, find a parametric form of the equation of the family of curves and show that 

dy/dx = (sin(t)sin(t/n)-cos(t)cos(t/n)) /( sin(t)cos(t/n)+cos(t)sin(t/n))

on maple i tried

x:=((cos(t/n))^n)*cos(t):

y:=((cos(t/n))^n)*sin(t):

w:=diff(x,t)

z:=diff(y,t)

z/w

and i never got the above answer so i did

simplify(z/w)

and still never got the answer instead i got 

(cos(t/n)*sin(t/n)-sin(t)*cos(t))/(cos(t/n)^2-cos(t)^2)

 

 

 

 

Hi,

I have a Maple code which produces an output plot for a first order differential equation,

nde := evalf(subs(npar, de));
nds := dsolve({nde, sigma(0) = -1e-8}, sigma(t), type = numeric);
acc_nds := (sig0, ae, re) -> dsolve({nde, sigma(0) = sig0}, sigma(t), type = numeric,
method = lsode[backfull], abserr = ae, relerr = re, maxfun = 0, ctrl=Ctrl);

odeplot(acc_nds(-0e-7, 1e-13, 1e-13), [t, sigma(t)], t = 0..2);

This produces the outputplot that I need for sigma(t). I need to produce a outut for d(sigmat)/dt, and how can this be done? what is the command I should use?

Additionally, how can i get the data set out of the polt in to a excel file or a text file?

I am quite new to maple, so i expect your kind support

Thanks

Hello all

 

I am new to Maple, and I am solving a system of two coupled partial differential equations using pdsolve, but I am having a hard time retrieving the solution evaluated at some point from the output. The output of pdsolve is a module, which appears to have different "methods" on it, including "plot3d" and "value". I can easily get a plot of my solution by using plot3d, but I don't know how to get a meaningful value out. For instance, if my solution is (f(x,y), g(x,y)), I would like to define H(x,y) = (f(x,y), g(x,y)), and be able to type H(10,10) into Maple to have my solution evaluated at that point. The result should be (1,1).

Here is a toy example:

firstEq := diff(f(x, y), x)+diff(f(x, y), y) = f(x, y)+g(x, y);
secondEq := diff(g(x, y), x)+diff(g(x, y), y) = 2*f(x, y)+g(x, y);
pdsystem := {firstEq, secondEq};

bv11 := f(10, y) = 1;
bv12 := f(x, 10) = 1;
bv21 := g(10, y) = 1;
bv22 := g(x, 10) = 1;
bvs := {bv11, bv12, bv21, bv22};
 
pdsolution := pdsolve(pdsystem, bvs, numeric, time = x, range = 0 .. 10);
 
pdsolution:-plot3d(x = 1 .. 10, y = 0 .. 10);
pdsolution:-value(10, 10);
Error, (in pdsolve/numeric/value) got additional unknown arguments {2}
 

Best regards.

 

I'm trying to plot the velocity of a ball thrown upwards with air resistance proportional to v^2 and also some simpler forms of this.

But the solution to v^2 returns root of and the plot stops for some specific time value. How can I proceed this plot to let's say 10 sec?

Staffan


 

``

``

restart

``

deq1 := m*(diff(v(t), t)) = -m*g:

``

sol := dsolve({deq1, v(0) = v__0}, v(t))

v(t) = -g*t+v__0

(1)

V := unapply(rhs(sol), t):

``

``

``

deq2 := m*(diff(v2(t), t)) = -m*g-k*v2(t):

``

sol2 := dsolve({deq2, v2(0) = v__0}, v2(t))

v2(t) = -g*m/k+exp(-k*t/m)*(v__0+g*m/k)

(2)

V2 := unapply(rhs(sol2), t):

``

deq3 := m*(diff(v3(t), t)) = -m*g-k*v3(t)*abs(v3(t))

m*(diff(v3(t), t)) = -m*g-k*v3(t)*abs(v3(t))

(3)

sol3 := dsolve({deq3, v3(0) = v__0}, v3(t))

v3(t) = RootOf(t+m*piecewise(_Z <= 0, arctanh(k*_Z/(k*m*g)^(1/2))/(k*m*g)^(1/2), 0 < _Z, arctan(k*_Z/(k*m*g)^(1/2))/(k*m*g)^(1/2))-m*piecewise(v__0 <= 0, arctanh(k*v__0/(k*m*g)^(1/2))/(k*m*g)^(1/2), 0 < v__0, arctan(k*v__0/(k*m*g)^(1/2))/(k*m*g)^(1/2)))

(4)

V3 := unapply(rhs(sol3), t):

``

m := 0.258e-2:

``

plot([V(t), V2(t), V3(t)], t = 0 .. 5, color = [blue, red, black], gridlines = true)

 

``


 

Download tal_3.9_sid_66_b.mw

The first half of this work sheet deals with SHM of pendulum. In the second half of the work sheet I attempt to solve for the general case of a swinging pendulum. Maple introduces a place holder (correct me if I have used the incorrect termonology) " __a" which I do not understand. What variable(s) should I replace it with and is there an automatic way of doing so?
 

restart

NULL

``

``

Simple*Harmonic*Motion*of*a*Pendulum

NULL

Eq1 := diff(Theta(t), t, t) = -omega^2*Theta(t)

diff(diff(Theta(t), t), t) = -omega^2*Theta(t)

(1)

ics := Theta(0) = 0, (D(Theta))(0) = Vmax

Theta(0) = 0, (D(Theta))(0) = Vmax

(2)

SHM := dsolve({Eq1, ics})

Theta(t) = Vmax*sin(omega*t)/omega

(3)

diffSHM := diff(SHM, t)

diff(Theta(t), t) = Vmax*cos(omega*t)

(4)

convert(diffSHM, D)

(D(Theta))(t) = Vmax*cos(omega*t)

(5)

eval[recurse](%, {t = 0, (D(Theta))(0) = Vmax})

Vmax = Vmax

(6)

solve(%, {_C1})

{_C1 = _C1}

(7)

assign(%); _C1

_C1

(8)

``

SHM

Theta(t) = Vmax*sin(omega*t)/omega

(9)

``

(General*Equation*of*Motion*of)*a*Pendulum

restart

diff(Theta(t), t, t) = -omega^2*sin(Theta(t))

diff(diff(Theta(t), t), t) = -omega^2*sin(Theta(t))

(10)

ics := Theta(0) = 0, (D(Theta))(0) = Vmax

Theta(0) = 0, (D(Theta))(0) = Vmax

(11)

Sol := dsolve(diff(Theta(t), t, t) = -omega^2*sin(Theta(t)))

Intat(1/(2*omega^2*cos(_a)+_C1)^(1/2), _a = Theta(t))-t-_C2 = 0, Intat(-1/(2*omega^2*cos(_a)+_C1)^(1/2), _a = Theta(t))-t-_C2 = 0

(12)

Sol[1]

Intat(1/(2*omega^2*cos(_a)+_C1)^(1/2), _a = Theta(t))-t-_C2 = 0

(13)

_C2 := 0

0

(14)

Sol[1]

Intat(1/(2*omega^2*cos(_a)+_C1)^(1/2), _a = Theta(t))-t = 0

(15)

``

dffSol[1] := diff(Sol[1], t)

(diff(Theta(t), t))/(2*omega^2*cos(Theta(t))+_C1)^(1/2)-1 = 0

(16)

``convert(dffSol[1], D)

(D(Theta))(t)/(2*omega^2*cos(Theta(t))+_C1)^(1/2)-1 = 0

(17)

 

eval[recurse](%, {t = 0, Theta(0) = 0, (D(Theta))(0) = Vmax})

Vmax/(2*omega^2+_C1)^(1/2)-1 = 0

(18)

solve(%, {_C1})

{_C1 = Vmax^2-2*omega^2}

(19)

assign(%); 1; _C1

Vmax^2-2*omega^2

(20)

dffSol[1]

(diff(Theta(t), t))/(2*omega^2*cos(Theta(t))+Vmax^2-2*omega^2)^(1/2)-1 = 0

(21)

``

dsolve(dffSol[1]); 1; SOL1 := int((diff(Theta(t), t))/sqrt(2*omega^2*cos(Theta(t))+Vmax^2-2*omega^2)-1, t = 0 .. Theta(t)) = 0

int((diff(Theta(t), t))/(2*omega^2*cos(Theta(t))+Vmax^2-2*omega^2)^(1/2)-1, t = 0 .. Theta(t)) = 0

(22)

Sol[1]

Intat(1/(2*omega^2*cos(_a)+Vmax^2-2*omega^2)^(1/2), _a = Theta(t))-t = 0

(23)

"Using  ( 1-cos(theta))/(2)=sin(theta/(2))^(2)and substituting by hand"

NotsoSHM := Intat(1/(Vmax*sqrt(1-2*omega^2*sin((1/2)*_a)^2/Vmax^2)), _a = (1/2)*Theta(t))-t = 0

Intat(1/(Vmax*(1-2*omega^2*sin((1/2)*_a)^2/Vmax^2)^(1/2)), _a = (1/2)*Theta(t))-t = 0

(24)

``

``


 

Download SHM_and_not_so_SHM.mw

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