<Text-field style="Heading 1" layout="Heading 1">Karakterligningen</Text-field>

Line\303\246re differentialligninger af 2. orden

Preben Alsholm 14/8 2007

Vi betragter i kurset kun line\303\246re differentialligninger med konstante koefficienter.

Vi har dog alligevel medtaget et eksempel p\303\245 en line\303\246r differentialligning med variable koefficienter (eksempel 4).

Medtaget er ogs\303\245 et eksempel p\303\245 en uline\303\246r differentialligning af 2. orden (pendulligningen).

with(plots): with(plottools):

<Text-field style="Heading 1" layout="Heading 1">Karakterligningen</Text-field>

xligning:=a*diff(x(t),t,t)+b*diff(x(t),t)+c*x(t)=0;

subs( x(t)=exp(R*t), xligning);

factor(%);

I DEtools-pakken ligger de2diffop der kan levere karakterpolynomiet s\303\245ledes:

DEtools[de2diffop](xligning,x(t),[R,t]);

<Text-field style="Heading 1" layout="Heading 1">Eksempel 1 To reelle r\303\270dder i karakterpolynomiet</Text-field>

En homogen ligning:

ligning1:=2*diff(x(t),t,t)+7*diff(x(t),t)-4*x(t)=0;

Den fuldst\303\246ndige l\303\270sning til denne homogene ligning ved brug af dsolve

dsolve(ligning1);

Mellemregning:

Karakterpolynomiet er

DEtools[de2diffop](ligning1,x(t),[R,t]);

der har r\303\270dderne

solve(%=0,R);

En inhomogen ligning:

ligning1i:=2*diff(x(t),t,t)+7*diff(x(t),t)-4*x(t)=-4*t^2+14*t+8;

Den fuldst\303\246ndige l\303\270sning til denne inhomogene ligning ved brug af dsolve

dsolve(ligning1i);

Er der givet begyndelsesbetingelserne x(0) = 1, x'(0) = -1, kan l\303\270sningen bestemmes ved brug af dsolve s\303\245ledes:

res:=dsolve({ligning1i,x(0)=1,D(x)(0)=-1});

Vi gemmer h\303\270jresiden i variablen X (stort!):

X:=subs(res,x(t));

plot(X,t=-1/2..1);

Et faseplansplot af l\303\270sningen er et plot af kurven (x(t), x'(t) ):

Xm:=diff(X,t);

p1:=plot([[1,-1]],style=point,symbol=circle,symbolsize=20,color=blue): p2:=plot([X,Xm,t=-1/2..2],labels=["x","x'"],title="Faseplansplot"): display(p1,p2);

Faseplansplottet animeret:

animate(plot,[[X,Xm,t=-1/2..T],labels=["x","x'"],caption="Faseplansplot"],T=-1/2..2,background=p1);

<Text-field style="Heading 2" layout="Heading 2">Samme ved simuleret h\303\245ndkraft</Text-field>

F\303\270rst skal karakterligningen nedskrives:

karlign:=DEtools[de2diffop](ligning1,x(t),[lambda,t])=0;

S\303\245 skal karakterligningen l\303\270ses:

lambda1,lambda2:=solve(karlign,lambda);

Den fuldst\303\246ndige l\303\270sning til den homogene ligning er nu, hvor r\303\270dderne er reelle og forskellige:

fuldhom:=c1*exp(lambda1*t)+c2*exp(lambda2*t);

Ansats til partikul\303\246r l\303\270sning til den inhomogene ligning

xp:=a*t^2+b*t+c;

Vi inds\303\246tter i den inhomogene ligning:

subs(x(t)=xp,ligning1i);

collect(%,t);

Dette skal v\303\246re en identitet i t, dvs. at ligningen skal g\303\246lde for alle t:

solve(identity(%,t),{a,b,c});

Hermed er en partikul\303\246r l\303\270sning til den inhomogene differentialligning givet ved

xpkonk:=subs(%,xp);

Den fuldst\303\246ndige l\303\270sning til den inhomogene ligning er dermed givet ved

fuldinhom:=x(t)=xpkonk+fuldhom;

Er der givet begyndelsesbetingelserne x(0) = 1, x'(0) = -1, kan l\303\270sningen bestemmes ved brug af simuleret h\303\245ndkraft s\303\245ledes:

F\303\270rst inds\303\246t t = 0 og x(0) = 1 i formlen for x(t):

L1:=subs(t=0,x(0)=1,fuldinhom);

L1;

Differenti\303\251r derefter udtrykket for x(t):

diff(fuldinhom,t);

Det er i det f\303\270lgende en fordel at bruge D:

convert(%,D);

Inds\303\246t t = 0 og x'(0) = -1:

L2:=subs(t=0,D(x)(0)=-1,%);

L2;

L\303\270s nu ligningerne L1 og L2 for c1 og c2:

solve({L1,L2},{c1,c2});

Inds\303\246t de fundne c-v\303\246rdier i den fuldst\303\246ndige l\303\270sning:

subs(%,fuldinhom);

<Text-field style="Heading 1" layout="Heading 1">Eksempel 2 En reel dobbeltrod i karakterpolynomiet</Text-field>

F\303\270rst en homogen ligning:

ligning2:=diff(x(t),t,t)+4*diff(x(t),t)+4*x(t)=0;

Den fuldst\303\246ndige l\303\270sning:

dsolve(ligning2);

En inhomogen ligning:

ligning2i:=diff(x(t),t,t)+4*diff(x(t),t)+4*x(t)=10*exp(-t);

Den fuldst\303\246ndige l\303\270sning:

dsolve(ligning2i);

Er der givet begyndelsesbetingelserne x(0) = -2, x'(0) = 0, kan l\303\270sningen bestemmes s\303\245ledes

res:=dsolve({ligning2i,x(0)=-2, D(x)(0)=0});

Vi gemmer h\303\270jresiden i variablen X (stort!):

X:=subs(res,x(t));

plot(X,t=-0.4..6);

Et faseplansplot af l\303\270sningen er et plot af kurven (x(t), x'(t) ):

Xm:=diff(X,t);

p1:=plot([[-2,0]],style=point,symbol=circle,symbolsize=20,color=blue): p2:=plot([X,Xm,t=-0.4..6],labels=["x","x'"],title="Faseplansplot"): display(p1,p2);

Faseplansplottet animeret:

animate(plot,[[X,Xm,t=-0.4..T],labels=["x","x'"],caption="Faseplansplot"],T=-0.4..6,frames=50,background=p1);

<Text-field style="Heading 2" layout="Heading 2">Simuleret h\303\245ndkraft</Text-field>

F\303\270rst skal karakterligningen nedskrives:

karlign:=DEtools[de2diffop](ligning2,x(t),[lambda,t])=0;

S\303\245 skal karakterligningen l\303\270ses:

lambda1,lambda2:=solve(karlign,lambda);

Vi ser, at der er en reel dobbeltrod. Den fuldst\303\246ndige l\303\270sning til den homogene ligning er nu:

fuldhom:=c1*exp(lambda1*t)+c2*t*exp(lambda1*t);

Ansats til partikul\303\246r l\303\270sning til den inhomogene ligning

xp:=a*exp(-t);

Vi inds\303\246tter i den inhomogene ligning:

subs(x(t)=xp,ligning2i);

Dette skal v\303\246re en identitet i t, dvs. at ligningen skal g\303\246lde for alle t:

solve(identity(%,t),{a});

Hermed er en partikul\303\246r l\303\270sning til den inhomogene differentialligning givet ved

xpkonk:=subs(%,xp);

Den fuldst\303\246ndige l\303\270sning til den inhomogene ligning er dermed givet ved

fuldinhom:=x(t)=xpkonk+fuldhom;

Er der givet begyndelsesbetingelserne x(0) = -2, x'(0) = 0, kan l\303\270sningen bestemmes ved brug af simuleret h\303\245ndkraft s\303\245ledes:

F\303\270rst inds\303\246t t = 0 og x(0) = -2 i formlen for x(t):

L1:=subs(t=0,x(0)=-2,fuldinhom);

L1;

Differenti\303\251r derefter udtrykket for x(t):

diff(fuldinhom,t);

Det er i det f\303\270lgende en fordel at bruge D:

convert(%,D);

Inds\303\246t t = 0 og x'(0) = 0:

L2:=subs(t=0,D(x)(0)=0,%);

L2;

L\303\270s nu ligningerne L1 og L2 for c1 og c2:

solve({L1,L2},{c1,c2});

Inds\303\246t de fundne c-v\303\246rdier i den fuldst\303\246ndige l\303\270sning:

subs(%,fuldinhom);

<Text-field style="Heading 1" layout="Heading 1">Eksempel 3 To imagin\303\246re r\303\270dder i karakterpolynomiet</Text-field>

En homogen differentialligning

ligning3:=diff(x(t),t,t)+4*diff(x(t),t)+13*x(t)=0;

Den fuldst\303\246ndige l\303\270sning:

dsolve(ligning3);

En inhomogen differentialligning:

ligning3i:=diff(x(t),t,t)+4*diff(x(t),t)+13*x(t)=4*cos(t);

Den fuldst\303\246ndige l\303\270sning:

dsolve(ligning3i);

Med begyndelsesbetingelserne x(0) = 1, x'(0) = 2:

res:=dsolve({ligning3i,x(0)=1,D(x)(0)=2});

Vi gemmer h\303\270jresiden i variablen X (stort!):

X:=subs(res,x(t));

plot(X,t=0..20);

Et faseplansplot af l\303\270sningen er et plot af kurven (x(t), x'(t) ):

Xm:=diff(X,t);

p1:=plot([[1,2]],style=point,symbol=circle,symbolsize=20,color=blue): p2:=plot([X,Xm,t=0..20],labels=["x","x'"],caption="Faseplansplot",scaling=constrained): display(p1,p2);

Faseplansplottet animeret:

animate(plot,[[X,Xm,t=0..T],labels=["x","x'"],caption="Faseplansplot"],T=0..20,frames=100,background=p1,scaling=constrained);

<Text-field style="Heading 2" layout="Heading 2">Simuleret h\303\245ndkraft</Text-field>

F\303\270rst skal karakterligningen nedskrives:

karlign:=DEtools[de2diffop](ligning3,x(t),[lambda,t])=0;

S\303\245 skal karakterligningen l\303\270ses:

lambda1,lambda2:=solve(karlign,lambda);

Vi ser, at der er 2 imagin\303\246re r\303\270dder, og at de er hinandens kompleks konjugerede. Den fuldst\303\246ndige l\303\270sning til den homogene ligning er nu:

fuldhom:=c1*Re(exp(lambda1*t))+c2*Im(exp(lambda1*t)) assuming t::real;

Ansats til partikul\303\246r l\303\270sning til den inhomogene ligning

xp:=a*cos(t)+b*sin(t);

Vi inds\303\246tter i den inhomogene ligning:

subs(x(t)=xp,ligning3i);

collect(%,[cos(t),sin(t)]);

Dette skal v\303\246re en identitet i t, dvs. at ligningen skal g\303\246lde for alle t:

solve(identity(%,t),{a,b});

Hermed er en partikul\303\246r l\303\270sning til den inhomogene differentialligning givet ved

xpkonk:=subs(%,xp);

Den fuldst\303\246ndige l\303\270sning til den inhomogene ligning er dermed givet ved

fuldinhom:=x(t)=xpkonk+fuldhom;

Er der givet begyndelsesbetingelserne x(0) = 1, x'(0) = 2, kan l\303\270sningen bestemmes ved brug af simuleret h\303\245ndkraft s\303\245ledes:

F\303\270rst inds\303\246t t = 0 og x(0) = 1 i formlen for x(t):

L1:=subs(t=0,x(0)=1,fuldinhom);

L1;

Hermed er c1 jo faktisk bestemt.

Differenti\303\251r derefter udtrykket for x(t):

diff(fuldinhom,t);

Det er i det f\303\270lgende en fordel at bruge D:

convert(%,D);

Inds\303\246t t = 0 og x'(0) = 2:

L2:=subs(t=0,D(x)(0)=2,%);

L2;

L\303\270s nu ligningerne L1 og L2 for c1 og c2:

solve({L1,L2},{c1,c2});

Inds\303\246t de fundne c-v\303\246rdier i den fuldst\303\246ndige l\303\270sning:

subs(%,fuldinhom);

<Text-field style="Heading 1" layout="Heading 1">Eksempel 4 En differentialligning med variable koefficienter</Text-field>

Differentialligninger med ikke-konstante koefficienter har ofte l\303\270sninger, der ikke kan udtrykkes ved de element\303\246re funktioner.

Lad LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVElJm51O0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRic= v\303\246re et reelt tal. Differentialligningen

ligning4:=t^2*diff(x(t),t,t)+t*diff(x(t),t)+(t^2-nu^2)*x(t)=0;

kaldes Bessels differentialligning. L\303\270sningerne er givet ved Besselfunktioner:

dsolve(ligning4);

plot([seq(BesselJ(k,t),k=0..5)],t=0..10,legend=[seq("BesselJ"||k,k=0..5)]);

plot([seq(BesselY(k,t),k=0..5)],t=0..10,-2..0.6,legend=[seq("BesselY"||k,k=0..5)]);

Maple ved en del om Besselfunktioner:

FunctionAdvisor(BesselJ);

<Text-field style="Heading 1" layout="Heading 1">Eksempel 5 En uline\303\246r differentialligning af 2. orden (pendulet)</Text-field>

Matematisk pendul. LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEoJnRoZXRhO0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRic= er udsvinget fra lodret stilling.

Newtons anden lov:

m*diff(L*theta(t),t,t)=-m*g*sin(theta(t));

Vi s\303\246tter LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2I1EhRictRiM2JS1GLDYlUSJrRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUkmbXNxcnRHRiQ2Iy1JJm1mcmFjR0YkNigtRiM2Iy1GLDYlUSJnRidGNEY3LUYjNiMtRiw2JVEiTEYnRjRGNy8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGY28vJSliZXZlbGxlZEdGQkYr . Dermed har vi:

ligning5:=diff(theta(t),t,t)=-k^2*sin(theta(t));

Den fuldst\303\246ndige l\303\270sning kommer ikke p\303\245 noget, der blot tiln\303\246rmelsesvis ligner eksplicit form:

dsolve(ligning5);

Vi fors\303\270ger alligevel at at finde den l\303\270sning, der svarer til at pendulet starter med et udsving p\303\245 30 grader fra lodlinien. Vi s\303\246tter k til 1:

dsolve({subs(k=1,ligning5),theta(0)=Pi/6,D(theta)(0)=0});

value([%]);

Det er nok en bedre id\303\251 at l\303\270se differentialligningen numerisk:

res:=dsolve({subs(k=1,ligning5),theta(0)=Pi/6,D(theta)(0)=0},type=numeric,output=listprocedure,range=0..20);

Vi gemmer den procedure, der beregner vinklen LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEoJnRoZXRhO0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRic= i Theta, og den der beregner vinkelhastigheden LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2I1EhRictRiM2Ji1JJm1mcmFjR0YkNigtRiM2Iy1JI21vR0YkNi1RMCZEaWZmZXJlbnRpYWxEO0YnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGPy8lKXN0cmV0Y2h5R0Y/LyUqc3ltbWV0cmljR0Y/LyUobGFyZ2VvcEdGPy8lLm1vdmFibGVsaW1pdHNHRj8vJSdhY2NlbnRHRj8vJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZOLUYjNiVGKy1GIzYkRjYtRiw2JVEidEYnLyUnaXRhbGljR1EldHJ1ZUYnL0Y7USdpdGFsaWNGJ0YrLyUubGluZXRoaWNrbmVzc0dRIjFGJy8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0Zcby8lKWJldmVsbGVkR0Y/LUY3Ni1RMSZJbnZpc2libGVUaW1lcztGJ0Y6Rj1GQEZCRkRGRkZIRkpGTEZPLUYjNiUtRiw2JVEoJnRoZXRhO0YnL0ZZRj9GOi1GNzYtUTAmQXBwbHlGdW5jdGlvbjtGJ0Y6Rj1GQEZCRkRGRkZIRkpGTEZPLUkobWZlbmNlZEdGJDYkLUYjNiNGVUY6RitGKw== gemmer vi i Thetam:

Theta,Thetam:=op(subs(res,[theta(t),diff(theta(t),t)]));

plot(Theta,0..20,title="Vinkeludsving som funktion af tiden");

Et faseplansplot af l\303\270sningen:

p1:=plot([[Pi/6,0]],style=point,symbol=circle,symbolsize=20,color=blue): p2:=plot([Theta,Thetam,0..20],title="Faseplansplot",scaling=constrained): display(p1,p2);

Faseplansplottet animeret:

animate(plot,[[Theta,Thetam,0..T],title="Faseplansplot"],T=0..7,frames=50,scaling=constrained,view=[-0.6..0.6,-0.6..0.6],background=p1);

Her f\303\270lger s\303\245 en animering af pendulet, hvor vi af hensyn til senere \303\246ndringer har taget alt med:

T0:=20: N:=150: theta0:=Pi/6: thetam0:=0: vw:=[-1..1,-1.2..0]: res:=dsolve({subs(k=1,ligning5),theta(0)=theta0,D(theta)(0)=thetam0},type=numeric,output=listprocedure,range=0..T0): Theta:=subs(res,theta(t)): TS:=seq(Theta(T0*j/N)-Pi/2,j=0..N): pp1:=seq(disk([cos(TS[j]),sin(TS[j])],.1,color=red),j=1..N+1): pp2:=seq(plot([[0,0],[1,TS[j]]],coords=polar,view=vw,thickness=2),j=1..N+1): pp:=seq(display(pp1[j],pp2[j]),j=1..N+1): display(pp,insequence=true,scaling=constrained);

I ovenst\303\245ende program kan man blot \303\246ndre T0, N, theta0, thetam0 og vw.

Det hele kan lidt mere bekvemt indpakkes i en procedure:

pendul:=proc(theta0::realcons,thetam0::realcons,T0::positive,N::posint,vw::list(range)) local res,Theta,TS,pp1,pp2,pp; global k,ligning5,theta,t; res:=dsolve({subs(k=1,ligning5),theta(0)=theta0,D(theta)(0)=thetam0}, type=numeric,output=listprocedure,range=0..T0); Theta:=subs(res,theta(t)); TS:=seq(Theta(T0*j/N)-Pi/2,j=0..N); pp1:=seq(disk([cos(TS[j]),sin(TS[j])],.1,color=red),j=1..N+1); pp2:=seq(plot([[0,0],[1,TS[j]]],coords=polar,view=vw,thickness=2),j=1..N+1); pp:=seq(display(pp1[j],pp2[j]),j=1..N+1); display(pp,insequence=true,scaling=constrained) end proc;

Vi afpr\303\270ver proceduren med et startudsving p\303\245 45 grader, en vinkelhastighed p\303\245 0, en sluttid p\303\245 20, med antal billeder lig 100 og med vindue givet ved intervallerne [-1..1,-1.2..0]:

pendul(Pi/4,0,20,100,[-1..1,-1.2..0]);

Med st\303\270rre udsving (90 grader) til at begynde med:

pendul(Pi/2,0,20,150,[-1.2..1.2,-1.2..0.1]);

Startudsving t\303\246t p\303\245 180 grader:

pendul(3,0,40,300,[-1.2..1.2,-1.2..1.2]);

Startudsving 180 grader, men positiv begyndelseshastighed stor nok til at rotere rundt:

pendul(3,-.145,40,300,[-1.2..1.2,-1.2..1.2]);

<Text-field style="Heading 1" layout="Heading 1">Eksempel 1 som system af 2 differentialligninger af f\303\270rste orden</Text-field>

En differentialligning af 2. orden kan omskrives til et system af 2 differentialligninger af f\303\270rste orden. Id\303\251en er den enkle, at indf\303\270re x'(t) som en ny ubekendt f.eks. ben\303\246vnt y(t) som i eksemplet nedenfor.

with(LinearAlgebra):

ligning1:=2*diff(x(t),t,t)+7*diff(x(t),t)-4*x(t)=0;

L1:=diff(x(t),t)=y(t);

Vi erstatter x'(t) med y(t) i ligning1:

L2:=subs(diff(x(t),t)=y(t),ligning1);

Det nye system er nu:

sys:=[L1,isolate(L2,diff(y(t),t))];

Nu kan systemet l\303\270ses som vi tidligere har gjort vha. egenv\303\246rdierne og egenvektorerne for matricen A givet ved

A:=Matrix([[0,1],[2,-7/2]]);

CharacteristicPolynomial(A,lambda);

EV,P:=Eigenvectors(A);

Den fuldst\303\246ndige l\303\270sning:

c1*exp(EV[1]*t)&*Column(P,1)+c2*exp(EV[2]*t)&*Column(P,2): %=eval(%,`&*`=`*`);

De andre eksempler ovenfor kan p\303\245 ganske tilsvarende m\303\245de omkrives til systemer af differentialligninger af f\303\270rste orden.