TENTAMEN 2012-03-17   Korta svar - del 1

1. NR är  x(n+1)=x(n) - f(x(n)) / f'(x(n))

3. IP genom tre punkter = grad 2. Handräkning så Newtons ansats:

   y=c1+c2*(x-0)+c3*(x-0)*(x-1)

   x=0 ger  y=1=c1
   x=1 ger  y=2=c1+c2*(1-0) --> c2=1
   x=2 ger  y=4=c1+c2*(2-0)+c3*(2-0)*(2-1)=1+2+c3*2 --> c3=1/2

   y(0.2)=c1+c2*(0.2-0)+c3*(0.2-0)*(0.2-1)=1+1*0.2+1/2*0.2*(-0.8)=1.2-0.08=1.12

4. 1 st 2:a ordningen & 1 st 1:a ordning blir totalt 1*2+1*1=3.

6. Flytta över termer så att x=G(x). Endast översta och understa i
   vänstra kolumnen är OK.

   Fixpunktsmetoden fungerar bara om |G'(x)|<1. Endast den understa
   verkar ha liten derivata.  I den översta trillar en faktor 7 ut.

7. A=[x x.^2]=[0 0; 1 1; 2 4]. Lös med normalekvationerna.
   r=b-A*c

8. Vi startar i x=1 med h=0.5, så

   y'(1.5) = y'(1.0) + h*y''(1.0) = 2 + 0.5*1 = 2.5

   y(1.5) = y(1.0) + h*y'(1.0) = 3 + 0.5*2 = 4
   y(2.0) = y(1.5) + h*y'(1.5) = 4 + 0.5*2.5 = 5.25

9. Tre delar ger h=2.

   T(h=2)= h * (f(2)/2 + f(4) + f(6) + f(8)/2) =
         = 2 * (cos(2*pi)/2 + cos(4*pi) + cos(6*pi) + cos(8*pi)/2) =
	 = 2 * (1/2 + 1 + 1 + 1/2) = 6

---------------------------------------------------------------------------
Del 2

P1a 

Det är en 3:e ordningens DE, som skrivs om till 3st första ordningens DE.

Inför   u1=y,  u2=y', u3=y''  och stega med Eulers metod: y(n+1)=y(n)+h*y'(n)

	x        u1         u2          u3      u1'=u2     u2'=u3   u3'= DE.

    1.00000    3.00000    2.00000    1.00000    2.00000    1.00000   -1.00000
    1.50000    4.00000    2.50000    0.50000    2.50000    0.50000   -4.00000
    2.00000    5.25000    2.75000   -1.50000    2.75000   -1.50000  -11.56250
    2.50000    6.62500    2.00000   -7.28125    2.00000   -7.28125  -23.89062


y'(1.5)=2.5,  y''(1.5)=-1.5

P1b  Endast de rader utan markeringar.  (De med *** = P1c,   +++ = P1d)

     n=1000; x0=1; xN=3; u0=[3 2 1];
     yy=[];                                         %*** För att kunna spara
     for varv=1:2;				    %*** För att köra flera n
       x=x0; u=u0; 
       xx=x; uu=u;                                  %+++ För att kunna spara
       for i=1:n;
         h=(xN-x0)/n;
         uprim=[u(2), u(3), 8*x-u(1)^2];
         u=u+h*uprim;
         x=x+h;
         xx=[xx;x]; uu=[uu;u];                      %+++ Spara alla delsteg
       end;
       yslut=u(1)
       yy=[yy;yslut];                               %*** Spara alla yslut
       n=2*n;                                       %*** Halverat steg i nästa
     end;                                           %***
     etrunk=abs(diff(yy))                           %*** Etrunk är skillnaden
     % Om Etrunk för stort, öka n i början av prg.  %***
     y=uu(:,1); ybis=uu(:,3); w=ybis./y;            %+++ Plocka ut rätt vektor
     plot(xx,w);	      			    %+++

Ursprungsprogrammet kör en Euler med 1000 delsteg.
I P1c körs Eulern två gångar, med olika n och kollar skillnaden i slut-y.
I P1d sparas alla värden vi stegat igenom i xx & uu.


P2a

Plottar man de givna punkterna ser man att de ligger rätt jämnt
fördelade över ellipsen. Bra skattningar för medelunkten 
blir då tex medelvärdet eller värdet mittemellan min-o-max.
Halvaxlarna som avståndet från medelpunkten till min eller max.

xm:  medel=4.4,  mitt=(2+8)/2=5   
ym:  medel=3.0,  mitt=(1+5)/2=3	  vilket ger   a=3 b=2 (lite grovt)

P2b

Om xm och ym är fixerade är problemet linjärt i 1/a^2 resp 1/b^2.
Lös med MKV det linjära problemet:

   c1*(x-xm)^2 + c2*(y-ym)^2 = 1

i Matlab:

  x=[2 2 4 6 8]';  y=[2 4 1 5 3]'; xm=...; ym=...; 
  A=[(x-xm).^2 (y-ym).^2]; b=ones(size(x));
  c=A\b;
  a=sqrt(1/c(1))
  b=sqrt(1/c(2))

P2c

Nu är det ickelinjärt, så det måste lösas med Gauss-Newton.
Startvärden enligt a-uppgiften.

c0=[5 3 3 2]; % startvarden [xm a ym b] 
t=1; c=c0(:); it=0;
while (norm(t)>1e-6)&(it<12);
  f=((x-c(1))/c(2)).^2+((y-c(3))/c(4)).^2-1;
  J=[-2*(x-c(1))/(c(2))^2,   -2*(x-c(1)).^2/(c(2))^3, ...
     -2*(y-c(3))/(c(4))^2,   -2*(y-c(3)).^2/(c(4))^3];
  t=J\f;
  tn=norm(t)
  c=c-t;
  it=it+1;
end;
c'

xm=c(1), ym=c(3), a=c(2), b=c(4);


P3a  Integralskattning och ekvationslösning.

P3b  Trapetsregeln för integralen och sekantmetoden för ekavtionslösningen.

P3c


  k0=2;  k1=3;

  h=0.1;
  x=0:h:4;
  y0=100*cos(k0*x)./(x+k0);
  int0=h*(sum(y0)-(y0(1)+y0(end))/2);
  f0=int0-k0^2;

  t=1;
  while abs(t)>1e-5;

    y1=100*cos(k1*x)./(x+k1);
    int1=h*(sum(y1)-(y1(1)+y1(end))/2);
    f1=int1-k1^2;
  
    t=f1*(k1-k0)/(f1-f0)
    k0=k1; f0=f1;
    k1=k1-t;
  end;
k=k1


P4a

  X=[0.0 0.3 ....]';  Y=[0.2 0.7 ....]';  % mätdata
  ij=1:4; x=X(ij); y=Y(ij);
  c=polyfit(x,y,3);
  xx=0.0:0.001:0.9;
  yy=polyval(c,xx);
  plot(xx,yy,x,y,'o');

P4b   fortsätter programmet:

  ij=4:9;
  plot(X(ij),Y(ij),'-o');

P4c   fortsätter programmet:

  baklykta=polyval(c,0.15)

  ij=8:9;
  k=polyval(X(ij),Y(ij),1);  % beräkna sista linjen
  framlykta=polyval(k,3.75)


P5a

  y'' + xy + x = y'  skrivs om till   y'' - y' + xy = -x

  Diskretisera med h=0.5: x0=2, x5=4.5; 
  Ersätt alla derivator med differenser:

  [y(i-1)-2*y(i)+y(i+1)]/h^2 - [y(i+1)-y(i-1)]/(2*h) + x(i)*y(i) = -x(i), i=1:4

  Samla ihop på index hos y:

  [1/h^2 + 1/(2*h)]*y(i-1) + [x(i)-2/h^2]*y(i) + [1/h^2 - 1/(2*h)]*y(i+1) =-x(i)

  med 1/h^2 = 4,  1/(2*h)=1

  Ger:

   [ x1-8   4-1   0      0    ]      [  y1 ]     [ -x1  - (4+1)*y0   ]
   [ 4+1    x2-8  4-1    0    ]      [  y2 ]     [ -x2               ]
   [ 0      4+1   x3-8   4-1  ]   *  [  y3 ]   = [ -x3               ]
   [ 0      0     4+1    x4-8 ]      [  y4 ]     [ -x4  - (4-1)*y5   ]

P5b

x0=2, x1=2.5, x2=3.0, x3=3.5, x4=4.0 och x5=4.5 så sökt y är y2.

P5c

FDM är av noggrannhetsordning 2, dvs felet är proportionellt mot h^2,
så det blir en fjärdedel så stort om steglängden halveras.

Den givna formeln är bara av ordning 1, dvs onödigt grov.
(En bättre formel finns tex i exempelsamlingens uppgift 1.3)