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));