Harjoitus 10 -- ratkaisut

1

Clear[m,g,k,h]

Solve[m*(fp1-2*f+fm1)/h^2==m*g-k*(fp1-fm1)/(2*h),fp1]

Remove[x,y]

eq1=m*y''[t]==m*g-k*y'[t]
eq2=m*x''[t]==-k*x'[t]

n=1000;t[0]=0.0; t[n]=5;h=(t[n]-t[0])/(n-1);

Do[t[i]=t[0]+i*h,{i,1,n}];

y[0]=0;y[1]=y[0]+h*20;

x[0]=0;x[1]=x[0]+h*10;

m=0.5;g=-9.81;k=0.1;

Do[
y[i]=(h*k-2*m)/(h*k+2*m)*y[i-2]+4*m/(h*k+2*m)*y[i-1]+2*g*m*h^2/(h*k+2*m);
x[i]=(h*k-2*m)/(h*k+2*m)*x[i-2]+4*m/(h*k+2*m)*x[i-1];,{i,2,n}]

ListPlot[Table[{x[i],y[i]},{i,0,n}],PlotRange->All]

2

eq=y''[x]+4x y'[x]+(102+4 x^2)*y[x]==0
r[x_]=y[x]/.Flatten[DSolve[{eq,y[0]==0,y'[0]==1},y[x],x]]

r[x_]=ComplexExpand[r[x]] //Simplify

Plot[r[x],{x,-3,3}, PlotRange->All]

3.

Remove[x,y]

n=100; x[0]=0.0; x[n]=3.0; x[-n]=-3.0; h=(x[n]-x[0])/(n-1);

Do[x[i]=x[0]+i*h; x[-i]=x[0]-i*h;,{i,1,n-1}]

y0=0;dy0=1;
y[0]=y0;
y[1]=y[0]+dy0*h;
y[-1]=y[0]-dy0*h;

Do[y[i]=1/(1/h^2+2*x[i-1]/h)*(y[i-1]*(2/h^2-102-4*x[i-1]^2)+y[i-2]*(-1/h^2+2*x[i-1]/h)); 
y[-i]=1/(1/h^2-2*x[-i+1]/h)*(y[-i+1]*(2/h^2-102-4*x[-i+1]^2)-y[-i+2]*(1/h^2+2*x[-i+1]/h)); ,{i,2,n}]

ListPlot[Table[{x[i],y[i]},{i,-n,n}],PlotRange->All]

exa[m_]:=(
Clear[x,y];
n=m;
x[0]=0.0; x[n]=3.0; x[-n]=-3.0; h=(x[n]-x[0])/(n-1);
Do[x[i]=x[0]+i*h; x[-i]=x[0]-i*h;,{i,1,n-1}];
y0=0;dy0=1;y[0]=y0;y[1]=y[0]+dy0*h;y[-1]=y[0]-dy0*h;
Do[y[i]=1/(1/h^2+2*x[i-1]/h)*(y[i-1]*(2/h^2-102-4*x[i-1]^2)+y[i-2]*(-1/h^2+2*x[i-1]/h)); 
y[-i]=1/(1/h^2-2*x[-i+1]/h)*(y[-i+1]*(2/h^2-102-4*x[-i+1]^2)-y[-i+2]*(1/h^2+2*x[-i+1]/h)); ,{i,2,n}];Show[ListPlot[Table[{x[i],y[i]},{i,-n,n}],PlotRange->All],Plot[r[t],{t,-3,3},PlotRange->All]])

exa[1000]

exa[500]

exa[100]

exa[50]

exa[30]

exa[20]