18) Differential Equations

Direction Fields and Solution Curves

Find the direction field of y'=y-t^2.

> with(DEtools):with(plots):

> diff(y(t),t)=y(t)-t^2;

> dfieldplot(diff(y(t),t)=y(t)-t^2,y(t),t=-3..3,y=-1..4);

> DEplot(diff(y(t),t)=y(t)-t^2, y(t), t=-3..3,{[0,0],[0,1],[0,2],[0,3],[0,-1]} ,y=-1..4);

Direction Fields and Isoclines(curves with constant slope)

> dfield:=dfieldplot(diff(y(t),t)=y(t)-t^2, y(t),t=-3..3,y=-1..4):

> iso:=implicitplot(y(t)-t^2=-1,t=-3..3,y=-1..4):

> display({dfield,iso});

Inflection curves (Curves given by y''=0);

> incurve:=implicitplot((y(t)-t^2)-2*t=0,t=-3..3,y=-1..4):

> display({dfield,incurve});

Laplace Transforms

Consider the following function:

> restart;

> g:=(S)->(S^3+2*S+1)/((S^2+1)*(S-2));

> YS:=convert(g(S),parfrac,S);

If Ys represents the Laplace transform of an IVP, the solution will be given by

> with(inttrans):

> Yt:=invlaplace(YS,S,t);

We can check the solution by taking the laplace transform

Resonance in an undamped vibrating system

Find the solution to the IVP y''+(omega^2)y=Ksin((gamma)t) , y(0)=0 , y'(0)=0.

> eq:=diff(y(t),t$2)+omega^2*y(t)=F*sin(gamma*t);

> Lap:=laplace(eq,t,s);

> Lap1:=subs(y(0)=0,D(y)(0)=0,Lap);

> Fs:=solve(Lap1,laplace(y(t),t,s));

> y:=invlaplace(Fs,s,t);

> Y:=simplify(y);

> subs(F=-12,gamma=2,omega=1,Y);

> plot(%,t=0..40);

> subs(F=-12,gamma=2,omega=2,Y);

Unit step function and Dirac's delta function.

Solve the IVP y''+3y'+2y=r(t) , y(0)=0 , y'(0)=0 where r(t) is:

a) a unit step in 0<t<1 ( r(t)= u(t)-u(t-1))

> restart;

> readlib(Heaviside):readlib(Dirac):

> r:=Heaviside(t)-Heaviside(t-1);

> eq:=diff(y(t),t$2)+3*diff(y(t),t)+2*y(t)=r;

> yp:=dsolve({eq,y(0)=0,D(y)(0)=0},y(t),method=laplace);

> with(plots):

> P1:=plot(r,t=0..3): P2:=plot(rhs(yp),t=0..3):

> display({P1,P2});

b) a hammer blow at t=1 (r(t)= delta(t-1))

> restart;

> r1:=Dirac(t-1);

> eq1:=diff(y(t),t$2)+3*diff(y(t),t)+2*y(t)=r1;

> yp1:=dsolve({eq1,y(0)=0,D(y)(0)=0},y(t),method=laplace);

> plot(rhs(yp1),t=0..4);

Bessel's DE

> x^2*diff(y(x),x,x) +x*diff(y(x),x)+x^2* y(x) = 0;

> dsolve(x^2*diff(y(x),x,x) +x*diff(y(x),x)+x^2* y(x) = 0, y(x));

> J1:=series(BesselJ(0,x),x=0,9);

> J1P:=convert(J1,polynom);

> Y1:=series(BesselY(0,x),x=0,9);

> Y1P:=convert(Y1,polynom);

> plot({BesselJ(0,x),BesselY(0,x)},x=0..10);

Partial Differential Equations

The Tricomi Equation

> y:='y':u:='u':y*diff(u(x,y),x,x)+diff(u(x,y),y,y)=0;

> u(x,y):=F(x)*G(y);

> eq:=y*diff(u(x,y),x,x)=-diff(u(x,y),y,y);

> eq0:=expand(eq/(y*u(x,y)));

> simplify(rhs(eq0)=-p^2);

> %*y*G(y);

> dsolve(%,G(y));