Fortran Numerical Analysis Programs

Kurdistan Region-Iraq Sulaimani University College of Science Physics Department

Numerical Analysis Programs Using Fortran 90 (2009)

Prepared by Dr. Omed Gh. Abdullah

1

Problem: Write a program to find the roots of the equation y = e x − 3 x using Bisection method in the interval [1,2] , within −3 the tolerance 10 .

f(x)=exp(x)-3*x tol=.001 a=1 b=2 10 c=(a+b)/2 print*,a,b,c if (f(c)*f(a))20,30,30 20 b=c goto 40 30 a=c 40 if (abs (a-b).lt.tol) goto 50 goto 10 50 print*,'The Root is',c end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1.000000 1.500000 1.500000 1.500000 1.500000 1.500000 1.500000 1.507813 1.511719 1.511719

2.000000 2.000000 1.750000 1.625000 1.562500 1.531250 1.515625 1.515625 1.515625 1.513672

1.500000 1.750000 1.625000 1.562500 1.531250 1.515625 1.507813 1.511719 1.513672 1.512695

The Root is 1.512695

2

Problem: Write a program to find the roots of the equation f ( x ) = x sin( x ) − 1 using False position method in the interval

[1,2] , within the tolerance 10−4 . f(x)=x*sin(x)-1 tol=.0001 a=1 b=2 10 c=(a*f(b)-b*f(a))/(f(b)-f(a)) print*,a,b,c if (f(a)*f(c))20,30,30 20 b=c goto 40 30 a=c 40 if (abs (a-b).lt.tol) goto 50 goto 10 50 print*,'The Root is',c end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1.000000 2.000000 1.162241 1.000000 1.162241 1.114254 1.000000 1.114254 1.114157 The Root is 1.114157

3

Problem: Write a program to find the roots of the equation f ( x ) = 3 x + sin x − e x using Secant method in the interval [0,1] , −4 within the tolerance 10 .

f(x)=3*x+sin(x)-exp(x) tol=.00001 x1=0 x2=1 10 x3=(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1)) print*,x1,x2,x3 x1=x2 x2=x3 if (abs (x1-x2).lt.tol) goto 20 goto 10 20 print*,'The Root is',x3 end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 0.000000 1.000000 4.709896E-01 3.075085E-01 3.626132E-01 3.604615E-01

1.000000 4.709896E-01 3.075085E-01 3.626132E-01 3.604615E-01 3.604217E-01

4.709896E-01 3.075085E-01 3.626132E-01 3.604615E-01 3.604217E-01 3.604217E-01

The Root is 3.604217E-01

4

Problem: Write a program to find the roots of the equation f ( x ) = x 2 − 2 x − 3 by iterative method, within the tolerance 10−4 . read*,xo tol=0.0001 5 x=sqrt(2*xo+3) print*,xo,x if(abs(x-xo)
3.316625 3.103748 3.034385 3.011440 3.003811 3.000423 3.000141 3.000047

The Root is 3.00047

5

Problem: Write a program using the Newton-Raphson procedure to determine the roots of the equation

f ( x) = x 3 − x 2 − 2 x + 1 ,

−6 within the tolerance 10 .

f(x)=x**3-x**2-2*x+1 df(x)=3*x**2-2*x-2 tol=1e-6 xo=2 k=0 10 k=k+1 x=xo-f(xo)/df(xo) if(abs(x-xo).lt.tol) goto 20 print *,k,x xo=x goto 10 20 print *,'The Root is',x end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1 2 3 4

1.833333 1.802935 1.801939 1.801938

The Root is 1.801938

If xo=0 The Root is 4.450419E-01

If xo=-1 The Root is -1.246980

6

Problem: Write a program to determine the roots of the equation x 4 − 1.99 x 3 − 1.76 x 2 + 5.22 x − 2.23 = 0 that are close to x = 1.5 by Newton-Raphson method. f(x)=x**4-1.99*x**3-1.76*x**2+5.22*x-2.23 df(x)=4*x**3-5.97*x**2-3.52*x+5.22 ddf(x)=12*x**2-11.94*x-3.52 tol=0.00001 x=1.5 5 rat=df(x)/ddf(x) x=x-rat if(abs(rat)-tol)10,10,5 10 eps=sqrt(-2*f(x)/ddf(x)) sign=1 do i=1,2 k=0 y=x+sign*eps 15 k=k+1 print *,y rat=f(y)/df(y) y=y-rat if (abs(rat)-tol)20,15,15 20 sign=-1 print*,'++++++++++++++++++++++++++++++++' enddo end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1.569134 1.565973 1.565888 ++++++++++++++++++++++++++ 1.428165 1.424016 1.424112 ++++++++++++++++++++++++++

7

Problem: Write a program to determine the roots of the equation x 4 − 1.99 x 3 − 1.76 x 2 + 5.22 x − 2.23 = 0 that are close to x = 1.5 by Newton-Raphson method. f(x)=x**4-1.99*x**3-1.76*x**2+5.22*x-2.23 df(x)=4*x**3-5.97*x**2-3.52*x+5.22 ddf(x)=12*x**2-11.94*x-3.52 tol=0.00001 xo=1.5 10 x=xo-df(xo)/ddf(xo) if(abs(x-xo).lt.tol) goto 20 xo=x goto 10 20 eps=sqrt(-2*f(x)/ddf(x)) y1=x+eps 30 y=y1-f(y1)/df(y1) if(abs(y-y1).lt.tol) goto 40 print *,y y1=y goto 30 40 print *,'The First Root is',y print*,'************************' y1=x-eps 50 y=y1-f(y1)/df(y1) if(abs(y-y1).lt.tol) goto 60 print *,y y1=y goto 50 60 print *,'The Second Root is',y end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1.565973 1.565888 The First Root is 1.565888 ************************************ 1.424016 1.424112 The Second Root is 1.424112 ++++++++++++++++++++++++++ 8

Problem: Write a program to determine the roots of the equation 8 x 3 − 20 x 2 − 26 x + 33 = 0 that are close to x = 0.8 using BirgeVieta method. dimension a(0:10),b(0:10),c(0:10) a(0)=8 a(1)=-20 a(2)=-26 a(3)=33 m=3 xo=0.8 tol=1e-6 5 k=k+1 b(0)=a(0); c(0)=b(0) do j=1,m b(j)=a(j)+xo*b(j-1) c(j)=b(j)+xo*c(j-1) enddo x=xo-b(m)/c(m-1) print *,k,x if(abs(x-xo)
2.860506 3.049551 2.983405 3.005704 2.998054 3.000665 2.999773 3.000078 2.999973 3.000009

The Root is 1.999989

3.000009

11

Problem: Write a program to solve the following system of equations, by using Newton-Raphson Method x 2 − y − 0 .2 = 0 ; y 2 − x − 0 .3 = 0 Initiate the computation with guesses of x = 1.2 and y = 1.2 . f(x,y)=x**2-y-.2 g(x,y)=y**2-x-.3 fx(x,y)=2*x fy(x,y)=-1 gx(x,y)=-1 gy(x,y)=2*y xo=1.2 yo=1.2 tol=0.0001 5 xj=fx(xo,yo)*gy(xo,yo)-fy(xo,yo)*gx(xo,yo) hh=(g(xo,yo)*fy(xo,yo)-f(xo,yo)*gy(xo,yo))/xj hk=(f(xo,yo)*gx(xo,yo)-g(xo,yo)*fx(xo,yo))/xj x=xo+hh y=yo+hk print*,x,y if((abs(x-xo)
1.221849 1.221601 1.221601

The Root is 1.192309

1.221601

12

Problem: Write a program to find the interpolated value for x = 1.15 , using Newton forward method, for these tabulated data. X

1

1.2

1.4

1.6

1.8

2.0

f

2.317

2.425

2.522

2.609

2.689

2.762

program Newton_Forward_Iterpolation dimension x(20),y(20),d(20,20) n=6 xx=1.15 data (x(i), i=1,6) /1,1.2,1.4,1.6,1.8,2/ data (y(i), i=1,6)/2.317,2.425,2.522,2.609,2.689,2.762/ print*,(x(i), i=1,6) print*,(y(i), i=1,6) s=y(1) do i=1,n-1 do j=1,n-i y(j)=y(j+1)-y(j) d(i,j)=y(j) enddo enddo r=(xx-x(1))/(x(2)-x(1)) do i=1,n-1 p=1 do j=0,i-1 p=p*(r-j) enddo f=1 do j=1,i; f=f*j; enddo s=s+d(i,1)*p/f enddo print*,'Interpolation of(',xx,')=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Interpolation of (1.15)= 2.398955

13

Problem: Write a program to find the interpolated value for x = 1.15 , using Newton forward method, for these tabulated data. X

1

1.2

1.4

1.6

1.8

2.0

f

2.317

2.425

2.522

2.609

2.689

2.762

dimension f(0:20),diff(0:20) xo=1 x=1.15 h=0.2 n=5 data (f(i), i=0,5)/2.317,2.425,2.522,2.609,2.689,2.762/ r=(x-xo)/h do i=0,n-1 diff(i)=f(i+1)-f(i) enddo coeff=r yx=f(0)+coeff*diff(0) do i=2,n coeff=coeff*(r-i+1)/(i) do j=0,n-i diff(j)=diff(j+1)-diff(j) enddo yx=yx+coeff*diff(0) enddo print*,yx end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Interpolation of (1.15)= 2.398955

14

Problem: Write a program to find the interpolated value for x = 0.55 , using Newton forward method, for these tabulated data. X

0.5

0.7

0.9

1.1

1.3

1.5

Y

0.47943

0.64422

0.78333

0.89121

0.96356

0.99749

dimension y(20),diff(20) x0=0.5 x=0.55 h=0.2 n=6 m=4 data (y(i), i=1,6) /0.47943,0.64422,0.78333,0.89121,0.96356,0.99749/ j=int((x-x0)/h-0.5) r=(x-x0)/h-j if ((j.ge.0).and.(j.lt.n-m)) goto 10 print*,'Not enough function values' goto 20 10 do i=1,m-1 diff(i)=y(j+i+1)-y(j+i) enddo coeff=r yx=y(1)+coeff*diff(1) do i=2,m-1 coeff=coeff*(r-i+1)/i do j=1,m-i diff(j)=diff(j+1)-diff(j) enddo yx=yx+coeff*diff(1) enddo print*,'Interpolation of(',xx,')=',yx 20 end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Interpolation of (0.55)= 5.227315 E-01

15

Problem: Write a program to find the interpolated value for x = 1.93 , using Newton backward method, for these tabulated data. X

1

1.2

1.4

1.6

1.8

2.0

f

2.317

2.425

2.522

2.609

2.689

2.762

program Newton_Backward_Iterpolation dimension x(20),y(20),d(20,20) n=6 xx=1.93 data (x(i), i=1,6) /1,1.2,1.4,1.6,1.8,2/ data (y(i), i=1,6)/2.317,2.425,2.522,2.609,2.689,2.762/ print*,(x(i), i=1,6) print*,(y(i), i=1,6) s=y(n) do i=1,n-1 do j=1,n-i y(j)=y(j+1)-y(j) d(i,j)=y(j) enddo enddo r=(xx-x(n))/(x(2)-x(1)) do i=1,n-1 p=1 do j=0,i-1 p=p*(r+j) enddo f=1 do j=1,i; f=f*j; enddo s=s+d(i,n-i)*p/f enddo print*,'Interpolation of(',xx,')=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Interpolation of(1.93)=2.737522

16

Problem: Write a program to find the interpolated value for x = 1.47 , using Newton backward method, for these tabulated data. x

0.5

0.7

0.9

1.1

1.3

1.5

y

0.47943

0.64422

0.78333

0.89121

0.96356

0.99749

program Newton_Backward_Iterpolation dimension x(20),y(20),d(20,20) n=6 xx=1.47 data (x(i), i=1,6)/.5,.7,.9,1.1,1.3,1.5/ data (y(i), i=1,6) /0.47943,0.64422,0.78333,0.89121,0.96356,0.99749/ print*,(x(i), i=1,6) print*,(y(i), i=1,6) s=y(n) do i=1,n-1 do j=1,n-i y(j)=y(j+1)-y(j) d(i,j)=y(j) enddo enddo r=(xx-x(n))/(x(2)-x(1)) do i=1,n-1 p=1 do j=0,i-1 p=p*(r+j) enddo f=1 do j=1,i; f=f*j; enddo s=s+d(i,n-i)*p/f enddo print*,'Interpolation of(',xx,')=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Interpolation of(1.47)= 9.9492 E-01

17

Problem: Write a program to find the interpolated value for x = 1.05 , using Stirling Interpolation method, for these tabulated data. X

0.5

0.7

0.9

1.1

1.3

1.5

y

0.48

0.64

0.78

0.89

0.96

0.99

program Stirling_Formula_Iterpolation dimension x(20),y(20),d(20,20) n=6 xx=1.05 data (x(i), i=1,6) /.5,.7,.9,1.1,1.3,1.5/ data (y(i), i=1,6)/.48,.64,.78,.89,.96,.99/ print*,(x(i), i=1,6) print*,(y(i), i=1,6) s=y(4) do i=1,n-1 do j=1,n-i y(j)=y(j+1)-y(j) d(i,j)=y(j) print*,d(i,j) enddo enddo r=(xx-x(4))/(x(2)-x(1)) p1=r*(d(1,3)+d(1,4))/2 p2=r**2/2*d(2,3) p3=r*(r**2-1)/6*(d(3,2)+d(3,3))/2 p4=r**2*(r**2-1)/24*d(4,2) s=s+p1+p2+p3+p4 print*,'Interpolation of(',xx,')=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Interpolation of(1.05)=8.660302 E-01

18

Problem: Write a program to find the interpolated value for x = 1.05 , using Bessel’s Interpolation method, for these tabulated data. X

0.5

0.7

0.9

1.1

1.3

1.5

y

0.48

0.64

0.78

0.89

0.96

0.99

program Bessels_Formula_Iterpolation dimension x(20),y(20),d(20,20) n=6 xx=1.05 data (x(i), i=1,6) /.5,.7,.9,1.1,1.3,1.5/ data (y(i), i=1,6)/.48,.64,.78,.89,.96,.99/ print*,(x(i), i=1,6) print*,(y(i), i=1,6) s=(y(3)+y(4))/2 do i=1,n-1 do j=1,n-i y(j)=y(j+1)-y(j) d(i,j)=y(j) print*,d(i,j) enddo enddo r=(xx-x(3))/(x(2)-x(1)) p1=(r-.5)*d(1,3) p2=r*(r-1)/2*(d(2,2)+d(2,3))/2 p3=r*(r-1)*(r-.5)/6*d(3,2) p4=r*(r-1)*(r+1)*(r-2)/24*(d(4,1)+d(4,2))/2 s=s+p1+p2+p3+p4 print*,'Interpolation of(',xx,')=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Interpolation of(1.05)=8.659448 E-01

19

Problem: Write a program to find the interpolated value for x = 0.7 , using divided difference formula, for these tabulated data. X F(x)

0 -1.25

0.5 -3.5

1 6

1.4 2.32

2 1.5

3.2 1.27

5.5 5.6

program Divided_Differences_Iterpolation dimension x(20),y(20),d(20,20) n=7 xx=.7 data (x(i), i=1,7) /0,.5,1,1.4,2,3.2,5.5/ data (y(i), i=1,7)/-1.25,-3.5,6,2.32,1.5,1.27,5.6/ s=y(1) do i=1,n-1 do j=1,n-i y(j)=(y(j+1)-y(j))/(x(j+i)-x(j)) d(i,j)=y(j) print*,d(i,j) enddo enddo do i=1,n-1 p=1 do j=1,i p=p*(xx-x(j)) enddo s=s+d(i,1)*p enddo print*,'Interpolation of(',xx,')=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Interpolation of (0.7) =2.290929

20

Problem: Write a program to find the interpolated value for x = 1.5 , using divided difference formula, for these tabulated data. X F(x)

1 0

2 1

3 4

4 6

5 10

program Divided_Differences_Iterpolation dimension x(20),y(20),d(20,20) n=5 xx=1.5 data (x(i), i=1,5)/1,2,3,4,5/ data (y(i), i=1,5)/0,1,4,6,10/ s=y(1) do i=1,n-1 do j=1,n-i y(j)=(y(j+1)-y(j))/(x(j+i)-x(j)) d(i,j)=y(j) enddo enddo do i=1,n-1 p=1 do j=1,i p=p*(xx-x(j)) enddo s=s+d(i,1)*p enddo print*,'Interpolation of(',xx,')=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Interpolation of (1.5) = -0.171875

21

Problem: Write a program to find the interpolated value for x = 3 , using Lagrangian Polynomial, from the following data. X F(x)

3.2 22

2.7 17.8

1 14.2

4.8 38.3

program Lagrange_Iterpolation dimension x(20),y(20) n=4 xx=3 data (x(i), i=1,4)/3.2,2.7,1,4.8/ data (y(i), i=1,4)/22,17.8,14.2,38.3/ s=0 do i=1,n p=1 do j=1,n if (i.eq.j) goto 5 p=p*(xx-x(j))/(x(i)-x(j)) 5 enddo s=s+p*y(i) enddo print*,'Interpolation of(',xx,')=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Interpolation of (3) = 20.21196

22

Problem: Write a program to determine the parameters a1 & a 2 so that

f ( x ) = a1 + a 2 x , fits the following data in least squares sense. X y

0 1

0.3 2.7

0.6 4.3

0.9 6

1.2 7.5

1.5 9

1.8 10.6

2.1 12

program List_Square_Fitting dimension x(20),y(20),f(20) n=8 data (x(i), i=1,8)/0,.3,.6,.9,1.2,1.5,1.8,2.1/ data (y(i), i=1,8)/1,2.7,4.3,6,7.5,9,10.6,12/ sx=0; sxx=0;sy=0;sxy=0 do i=1,n sx=sx+x(i) sy=sy+y(i) sxx=sxx+x(i)**2 sxy=sxy+x(i)*y(i) enddo print *,'sx=',sx print *,'sxx=',sxx print *,'sy=',sy print *,'sxy=',sxy d=n*sxx-sx**2 a1=(sxx*sy-sx*sxy)/d a2=(n*sxy-sx*sy)/d print*,'a1=',a1 print*,'a2=',a2 s=0 do i=1,n f(i)=a1+a2*x(i) s=s+(y(i)-f(i))**2 enddo print*,'Standerd Deviation=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

sx=8.4

sxx=12.6

sy=53.1

a1=1.133333 a2=5.242064 Standard deviation = 6.726178 E -02 23

sxy=75.57

Problem: Write a program to determine the parameters A & B so that f ( x ) = A ln( x ) + B , fits the following data in least squares sense. X 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 y 0.27 0.72 1.48 2.66 4.48 7.26 11.43 17.64 26.78 program List_Square_Fitting dimension x(20),y(20),f(20) n=9 data (x(i), i=1,9)/.1,.2,.3,.4,.5,.6,.7,.8,.9/ data (y(i), i=1,9)/.27,.72,1.48,2.66,4.48,7.26,11.43,17.64,26.78/ sx=0; sxx=0;sy=0;sxy=0 do i=1,n sx=sx+log(x(i)) sy=sy+y(i) sxx=sxx+(log(x(i)))**2 sxy=sxy+log(x(i))*y(i) enddo print *,'sx=',sx print *,'sxx=',sxx print *,'sy=',sy print *,'sxy=',sxy d=n*sxx-sx**2 a1=(sxx*sy-sx*sxy)/d a2=(n*sxy-sx*sy)/d a=a2; b=a1 print*,'A=',a print*,'B=',b s=0 do i=1,n f(i)=a*log(x(i))+b s=s+(y(i)-f(i))**2 enddo print*,'Standerd Deviation=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

A=9.74113

B=16.66255

S=260.5141 24

Problem: Write a program to determine the parameters A & C so that

f ( x ) = C e A x , fits the following data in least squares sense. X 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 y 0.27 0.72 1.48 2.66 4.48 7.26 11.43 17.64 26.78 program List_Square_Fitting dimension x(20),y(20),f(20) n=9 data (x(i), i=1,9)/.1,.2,.3,.4,.5,.6,.7,.8,.9/ data (y(i), i=1,9)/.27,.72,1.48,2.66,4.48,7.26,11.43,17.64,26.78/ sx=0; sxx=0;sy=0;sxy=0 do i=1,n sx=sx+x(i) sy=sy+log(y(i)) sxx=sxx+(x(i))**2 sxy=sxy+x(i)*log(y(i)) enddo print *,'sx=',sx print *,'sxx=',sxx print *,'sy=',sy print *,'sxy=',sxy d=n*sxx-sx**2 a1=(sxx*sy-sx*sxy)/d a2=(n*sxy-sx*sy)/d a=a2; c=exp(a1) print*,'A=',a print*,'C=',c s=0 do i=1,n f(i)=c*exp(a*x(i)) s=s+(y(i)-f(i))**2 enddo print*,'Standerd Deviation=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

A= 5.512738 C= .2359106 S=52.53117 25

Problem: Write a program to determine the parameters A & C so that f ( x ) = C x A , fits the following data in least squares sense. X 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 y 0.27 0.72 1.48 2.66 4.48 7.26 11.43 17.64 26.78 program List_Square_Fitting dimension x(20),y(20),f(20) n=9 data (x(i), i=1,9)/.1,.2,.3,.4,.5,.6,.7,.8,.9/ data (y(i), i=1,9)/.27,.72,1.48,2.66,4.48,7.26,11.43,17.64,26.78/ sx=0; sxx=0;sy=0;sxy=0 do i=1,n sx=sx+log(x(i)) sy=sy+log(y(i)) sxx=sxx+(log(x(i)))**2 sxy=sxy+log(x(i))*log(y(i)) enddo print *,'sx=',sx print *,'sxx=',sxx print *,'sy=',sy print *,'sxy=',sxy d=n*sxx-sx**2 a1=(sxx*sy-sx*sxy)/d a2=(n*sxy-sx*sy)/d a=a2; c=exp(a1) print*,'A=',a print*,'C=',c s=0 do i=1,n f(i)=c*x(i)**a s=s+(y(i)-f(i))**2 enddo print*,'Standerd Deviation=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

A= 2.092605 C= 23.42711 S= 75.07242 26

Problem: Write a program to determine the parameters C & D so that D x f ( x ) = C x e , fits the following data in least squares sense. X y

0.1 0.27

0.2 0.72

0.3 1.48

0.4 2.66

0.5 4.48

0.6 7.26

0.7 11.4 3

0.8 17.6 4

0.9 26.7 8

program List_Square_Fitting dimension x(20),y(20),f(20) n=9 data (x(i), i=1,9)/.1,.2,.3,.4,.5,.6,.7,.8,.9/ data (y(i), i=1,9)/.27,.72,1.48,2.66,4.48,7.26,11.43,17.64,26.78/ sx=0; sxx=0;sy=0;sxy=0 do i=1,n sx=sx+x(i) sy=sy+log(y(i)/x(i)) sxx=sxx+x(i)**2 sxy=sxy+x(i)*log(y(i)/x(i)) enddo print *,'sx=',sx print *,'sxx=',sxx print *,'sy=',sy print *,'sxy=',sxy d=n*sxx-sx**2 a1=(sxx*sy-sx*sxy)/d a2=(n*sxy-sx*sy)/d d=a2; c=exp(a1) print*,'D=',d print*,'C=',c s=0 do i=1,n f(i)=c*x(i)*exp(d*x(i)) s=s+(y(i)-f(i))**2 enddo print*,'Standerd Deviation=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

C= 1.993407 D= 3.004763 S= 1.117279 E-01 27

Problem: Write a program to determine the parameters A & B so that A f ( x) = + B , fits the following data in least squares sense. x X y

0.5 3.3

2 2

3.5 1.4

4.1 1.2

5 1

6.2 0.8

7.5 0.64

9.2 0.5

program List_Square_Fitting dimension x(20),y(20),f(20) n=8 data (x(i), i=1,8)/.5,2,3.5,4.1,5,6.2,7.5,9.2/ data (y(i), i=1,8)/3.3,2,1.4,1.2,1,.8,.64,.5/ sx=0; sxx=0;sy=0;sxy=0 do i=1,n sx=sx+1/x(i) sy=sy+y(i) sxx=sxx+(1/x(i))**2 sxy=sxy+y(i)*(1/x(i)) enddo print *,'sx=',sx print *,'sxx=',sxx print *,'sy=',sy print *,'sxy=',sxy d=n*sxx-sx**2 a1=(sxx*sy-sx*sxy)/d a2=(n*sxy-sx*sy)/d a=a2; b=a1 print*,'A=',a print*,'B=',b s=0 do i=1,n f(i)=a/x(i)+b s=s+(y(i)-f(i))**2 enddo print*,'Standerd Deviation=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

A= 1.353131 B= .7405201 S=.7070401 28

Problem: Write a program to determine the parameters D & C so that D f ( x) = , fits the following data in least squares sense. x+C X y

0.5 3.3

2 2

3.5 1.4

4.1 1.2

5 1

6.2 0.8

7.5 0.64

9.2 0.5

program List_Square_Fitting dimension x(20),y(20),f(20) n=8 data (x(i), i=1,8)/.5,2,3.5,4.1,5,6.2,7.5,9.2/ data (y(i), i=1,8)/3.3,2,1.4,1.2,1,.8,.64,.5/ sx=0; sxx=0;sy=0;sxy=0 do i=1,n sx=sx+x(i)*y(i) sy=sy+y(i) sxx=sxx+(x(i)*y(i))**2 sxy=sxy+(x(i)*y(i))*y(i) enddo print *,'sx=',sx print *,'sxx=',sxx print *,'sy=',sy print *,'sxy=',sxy d=n*sxx-sx**2 a1=(sxx*sy-sx*sxy)/d a2=(n*sxy-sx*sy)/d c=-1/a2; d=a1*c print*,'C=',c print*,'D=',d s=0 do i=1,n f(i)=d/(x(i)+c) s=s+(y(i)-f(i))**2 enddo print*,'Standerd Deviation=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

C= 1.369226 D= 6.209051 S=5.723841 E-02 29

Problem: Write a program to determine the parameters A & B so that 1 f ( x) = , fits the following data in least squares sense. A x+B X y

0.5 3.3

2 2

3.5 1.4

4.1 1.2

5 1

6.2 0.8

7.5 0.64

9.2 0.5

program List_Square_Fitting dimension x(20),y(20),f(20) n=8 data (x(i), i=1,8)/.5,2,3.5,4.1,5,6.2,7.5,9.2/ data (y(i), i=1,8)/3.3,2,1.4,1.2,1,.8,.64,.5/ sx=0; sxx=0;sy=0;sxy=0 do i=1,n sx=sx+x(i) sy=sy+1/y(i) sxx=sxx+x(i)**2 sxy=sxy+x(i)*1/y(i) enddo print *,'sx=',sx print *,'sxx=',sxx print *,'sy=',sy print *,'sxy=',sxy d=n*sxx-sx**2 a1=(sxx*sy-sx*sxy)/d a2=(n*sxy-sx*sy)/d a=a2; b=a1 print*,'A=',a print*,'B=',b s=0 do i=1,n f(i)=1/(a*x(i)+b) s=s+(y(i)-f(i))**2 enddo print*,'Standerd Deviation=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

A= 1.953443 E-01 B= 9.250808 E-02 S= 3.864387 30

Problem: Write a program to determine the parameters A & B so that x f ( x) = , fits the following data in least squares sense. A x+B X y

0.5 3.3

2 2

3.5 1.4

4.1 1.2

5 1

6.2 0.8

7.5 0.64

9.2 0.5

program List_Square_Fitting dimension x(20),y(20),f(20) n=8 data (x(i), i=1,8)/.5,2,3.5,4.1,5,6.2,7.5,9.2/ data (y(i), i=1,8)/3.3,2,1.4,1.2,1,.8,.64,.5/ sx=0; sxx=0;sy=0;sxy=0 do i=1,n sx=sx+1/x(i) sy=sy+1/y(i) sxx=sxx+(1/x(i))**2 sxy=sxy+1/(x(i)*y(i)) enddo print *,'sx=',sx print *,'sxx=',sxx print *,'sy=',sy print *,'sxy=',sxy d=n*sxx-sx**2 a1=(sxx*sy-sx*sxy)/d a2=(n*sxy-sx*sy)/d a=a1; b=a2 print*,'A=',a print*,'B=',b s=0 do i=1,n f(i)=x(i)/(a*x(i)+b) s=s+(y(i)-f(i))**2 enddo print*,'Standerd Deviation=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

A= 1.279117 B= -5.697291 E-01 S= 16.410690 31

Problem: Write a program to determine the parameters A & B so that 1 f ( x) = 2 , fits the following data in least squares sense. A x+B

(

X y

)

0.5 3.3

2 2

3.5 1.4

4.1 1.2

5 1

6.2 0.8

7.5 0.64

9.2 0.5

program List_Square_Fitting dimension x(20),y(20),f(20) n=8 data (x(i), i=1,8)/.5,2,3.5,4.1,5,6.2,7.5,9.2/ data (y(i), i=1,8)/3.3,2,1.4,1.2,1,.8,.64,.5/ sx=0; sxx=0;sy=0;sxy=0 do i=1,n sx=sx+x(i) sy=sy+(y(i))**(-.5) sxx=sxx+x(i)**2 sxy=sxy+x(i)*(y(i))**(-.5) enddo print *,'sx=',sx print *,'sxx=',sxx print *,'sy=',sy print *,'sxy=',sxy d=n*sxx-sx**2 a1=(sxx*sy-sx*sxy)/d a2=(n*sxy-sx*sy)/d a=a2; b=a1 print*,'A=',a print*,'B=',b s=0 do i=1,n f(i)=1/(a*x(i)+b)**2 s=s+(y(i)-f(i))**2 enddo print*,'Standerd Deviation=',s end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

A= 9.919496 E-02 B= 5.035566 E-01 S= 2.275830 E-03 32

Problem: Write a Fortran program to find the value of the integral 1 dt ∫0 (t 2 + 1) (3t 2 + 4) using Trapezoidal rule with n = 6 program Trapezoidal_Rule f(t)=1/((t**2+1)*(3*t**2+4))**.5 a=0 b=1 n=6 h=(b-a)/n s=0 do i=1,n-1 x=a+h*i s=s+f(x) enddo t=h*(f(a)/2+s+f(b)/2) print*,'Integration by Trapezoidal Rule=',t end

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Integration by Trapezoidal rule = 4.016085 E-01

33

Problem: Write a Fortran program to find the value of the integral π

2

dx ∫0 sin2 x + 14 cos2 x

using Trapezoidal rule with n = 10

dimension y(20) f(x)=1/((sin(x))**2+1./4*(cos(x))**2) pi=22./7 a=0; b=pi/2; n=10 h=(b-a)/n do i=1,n+1 x=a+h*(i-1) y(i)=f(x) print*,y(i) enddo t=0 do i=1,n t=t+h*(y(i)+y(i+1))/2 enddo print*,'Integration by Trapezoidal rule =',t end

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Integration by Trapezoidal rule = 3.142227

34

Problem: Write a Fortran program to find the value of the integral 1 dt ∫0 (t 2 + 1) (3t 2 + 4) using Simpson’s 1/3 rule with n = 6 ! program Simpson's 1/3 Rule f(t)=1/((t**2+1)*(3*t**2+4))**.5 a=0 b=1 n=6 h=(b-a)/n s1=0;s2=0 do i=1,n-1,2 x=a+h*i s1=s1+f(x) enddo do i=2,n-1,2 x=a+h*i s2=s2+f(x) enddo s=(h/3)*(f(a)+4*s1+2*s2+f(b)) print*,'Integration by Simpsons 1/3 Rule=',s end

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Integration by Simpson's 1/3 rule= 4.021834 E-01

35

Problem: Write a Fortran program to find the value of the integral 1 dt ∫0 (t 2 + 1) (3t 2 + 4) using Simpson’s 1/3 rule with n = 6 ! program Simpson's 1/3 Rule dimension y(20) f(t)=1/((t**2+1)*(3*t**2+4))**.5 a=0 b=1 n=6 h=(b-a)/n s=0 do i=1,n+1 x=a+h*(i-1) y(i)=f(x) enddo do i=1,n,2 s=s+(h/3)*(y(i)+4*y(i+1)+y(i+2)) enddo print*,'Integration by Simpsons 1/3 Rule=',s end

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Integration by Simpson's 1/3 rule= 4.021834 E-01

36

Problem: Write a Fortran program to find the value of the integral 1 dt ∫0 (t 2 + 1) (3t 2 + 4) using Simpson’s 3/8 rule with n = 6 ! program Simpson's 3/8 Rule dimension y(0:10) f(t)=1/((t**2+1)*(3*t**2+4))**.5 a=0 b=1 n=6 h=(b-a)/n s=0 do i=0,n x=a+h*i y(i)=f(x) enddo do i=1,(n/2-1) s=s+y(3*i-3)+3*(y(3*i-2)+y(3*i-1))+y(3*i) enddo s=(3*h/8)*s print*,'Integration by Simpsons 3/8 Rule=',s end

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Integration by Simpson’s 3/8 rule= 4.021832 E-01

37

Problem: Write a Fortran program to find the value of the integral 1 dt ∫0 (t 2 + 1) (3t 2 + 4) using Simpson’s 3/8 rule with n = 6 ! program Simpson's 3/8 Rule dimension y(0:10) f(t)=1/((t**2+1)*(3*t**2+4))**.5 a=0 b=1 n=6 h=(b-a)/n s=0 do i=1,n+1 x=a+h*(i-1) y(i)=f(x) enddo do i=1,n,3 s=s+(3*h/8)*(y(i)+3*y(i+1)+3*y(i+2)+y(i+3)) enddo print*,'Integration by Simpsons 3/8 Rule=',s end

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Integration by Simpson’s 3/8 rule= 4.021832 E-01

38

Problem: Write a Fortran program to find the value of the integral 0.8 t 2 e dt 0

∫

using Trapezium rule and Romberg integration

with error of order

h8 .

dimension y(20),trap(20),romb(20,20) f(t)=exp(t**2) a=0; b=0.8; m=4 do k=1,m n=2**k h=(b-a)/n do i=1,n+1 x=a+h*(i-1) y(i)=f(x) enddo t=0 do i=1,n t=t+h*(y(i)+y(i+1))/2 enddo trap(k)=t print*,t enddo print*,'---------------------------' fact=2**2 do k=1,m-1 fact=fact**k do j=1,m-k trap(j)=(fact*trap(j+1)-trap(j))/(fact-1) romb(k,j)=trap(j) print*,romb(k,j) enddo print*,'---------------------------' enddo end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

39

1.048701 1.019178 1.011646 1.009753 -------------------------------------1.009338 1.009135 1.009122 -------------------------------------1.009122 1.009121 -------------------------------------1.009121

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

40

Problem: Write a Fortran program to find the value of the integral 0.8 t 2 e dt 0

∫

using Trapezium rule and Romberg integration

with error of order

h8 .

dimension trap(20,20) f(t)=exp(t**2) a=0; b=0.8; n=2; tol=1e-8 h=(b-a)/n sum=(f(a)+f(b))/2 do i=1,n-1 sum=sum+f(a+i*h) enddo trap(1,1)=h*sum print*,trap(1,1) j=1 5 do k=0,n-1 sum=sum+f(a+k*h+h/2) enddo n=2*n j=j+1 factor=1 h=h/2 trap(1,j)=h*sum print*,trap(1,j) do i=1,j-1 factor=2**2*factor trap(i+1,j)=(factor*trap(i,j)-trap(i,j-1))/(factor-1) enddo if (abs(trap(j,j)-trap(j-1,j)).lt.tol) goto 15 goto 5 15 print*,'+++++++++++++++++++++++++++' do k=1,j do i=1,k print*,trap(i,k) enddo print*,'---------------------------' enddo end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 41

1.048701 1.019178 1.011646 1.009753 +++++++++++++++++++++ 1.048701 -------------------------------------1.019178 1.009338 -------------------------------------1.011646 1.009135 1.009122 -------------------------------------1.009753 1.009122 1.009121 1.009121 --------------------------------------

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

42

Problem: Write a Fortran program to find the value of the integral 0.8 t 2 e dt 0

∫

using Simpson’s 1/3 rule and Romberg integration

with error of order

h8 .

dimension y(20),simps(20),romb(20,20) f(t)=exp(t**2) a=0; b=0.8; m=4 do k=1,m n=2**k h=(b-a)/n do i=1,n+1 x=a+h*(i-1) y(i)=f(x) enddo s=0 do i=1,n,2 s=s+(h/3)*(y(i)+4*y(i+1)+y(i+2)) enddo simps(k)=s print*,s enddo print*,'---------------------------' fact=2**2 do k=1,m-1 fact=fact**k do j=1,m-k simps(j)=(fact*simps(j+1)-simps(j))/(fact-1) romb(k,j)=simps(j) print*,romb(k,j) enddo print*,'---------------------------' enddo end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

43

1.012070 1.009338 1.009135 1.009122 -------------------------------------1.008427 1.009067 1.009117 -------------------------------------1.009110 1.009121 -------------------------------------1.009121

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

44

Problem: Write a Fortran program to find first, second, third, and fourth at x = 1 and h = 0.01 derivation of the function f ( x) = ln( x 3 )

! Derivation using Different Formula f(x) = log(x**3) x=1 h=0.01 ! df: Forward, db: Backward, dc: Central df=(f(x+h)-f(x))/h db=(f(x)-f(x-h))/h dc=(f(x+h)-f(x-h))/(2*h) print*,'Forward derivation is', df print*,'Backward derivation is', db print*,'Central derivation is', dc print*,'-----------------------------------------' ! d3p: Three Point, d5p: Five Point d3p=(f(x+h)-f(x-h))/(2*h) d5p=(-1* f(x+2*h)+8*f(x+h)-8*f(x-h)+f(x-2*h))/(12*h) print*,'Three Point derivation is', d3p print*,'Five Point derivation is', d5p print*,'-----------------------------------------' ! d2: Second derivation, d3: Third Derivation, d4: Fourth Derivation" d2=(f(x+h)-2*f(x)+f(x-h))/h**2 d3=(f(x+2*h)-2*f(x+h)+2*f(x-h)-f(x-2*h))/(2*h**3) d4=(f(x+2*h)-4*f(x+h)+6*f(x)-4*f(x-h)+f(x-2*h))/h**4 print*,'Second derivation is', d2 print*,'Third derivation is ', d3 print*,'Fourth derivation is ', d4 end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

45

Forward derivation is 2.985099 Backward derivation is 3.015098 Central derivation is 3.000097 ------------------------------------------------------------------Three Point derivation is 3.000097 Five Point derivation is 2.999997 ------------------------------------------------------------------Second derivation is -3.000144 Third derivation is 6.001784 Fourth derivation is -18.005940

46

Problem: Write a Fortran program to find first derivation from the following data using derivation of Newton Forward polynomial. x 1 2 3 4 5 6 7 y = x3 − 2 x 2 + 2 x − 5 y -4 -1 10 35 80 151 254 ! Derivative from Newton Forward dimension x(10),y(10),d(10,10) n=7 data (x(i), i=1,7)/1,2,3,4,5,6,7/ data (y(i), i=1,7)/-4,-1,10,35,80,151,254/ !print*,(x(i), i=1,7) !print*,(y(i), i=1,7) do i=1,n-1 do j=1,n-i y(j)=y(j+1)-y(j) d(i,j)=y(j) enddo enddo s=0 do i=1,n-1 s=s+d(i,1)/i*(-1)**(i+1) enddo fd=s/(x(2)-x(1)) print*,'First derivation is', fd end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The First Derivation= 1.000000

47

Problem: Write a Fortran program to find first derivation from the following data at x = 0.1 using derivation of Lagrange polynomial. x 0.1 0.2 0.3 0.4 y 0.01 0.04 0.09 0.16

! Derivation using Lagrange Formula dimension x(40),y(40) n=4 xx=0.1 data (x(i), i=1,4)/.1,.2,.3,.4/ data (y(i), i=1,4)/.01,.04,.09,.16/ d=0 do i=1,n s=0 do j=1,n if (j.eq.i) goto 10 p=1 do k=1,n if ((k.eq.j).or.(k.eq.i)) goto 20 p=p*(xx-x(k))/(x(i)-x(k)) 20 enddo s=s+p/(x(i)-x(j)) 10 enddo d=d+s*y(i) enddo print*,'The First Derivation=',d end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The First Derivation= 2.000000E-01

48

Problem: Write a Fortran program to find first derivation of the function f ( x) = ln(5x 2 ) + 8x 3 − 8 at x = 1.1 using derivation of Lagrange polynomial. ! Derivation using Lagrange Formula f(x)=log(5*x**2)+8*x**3-8 n=4 xo=1.1 h=.1 d=0 do i=1,n s=0 do j=1,n if (j.eq.i) goto 10 p=1 do k=1,n if ((k.eq.j).or.(k.eq.i)) goto 20 p=p*(1-k)/(i-k) 20 enddo s=s+p/((i-j)*h) 10 enddo d=d+s*f(xo+h*(i-1)) enddo print*,'The First Derivation=',d end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The First Derivation= 30.856790

49

Problem: Write a Fortran program to find first derivation of the function f ( x) = ln(5x 2 ) + 8x 3 − 8 at x = 1.1 using derivation of Lagrange polynomial.

! Derivation using Lagrange Formula f(x)=log(5*x**2)+8*x**3-8 n=4 xo=1.1 h=.1 d=0 do i=1,n s=0 do j=1,n if (j.eq.i) goto 10 p=1 do k=1,n if ((k.eq.j).or.(k.eq.i)) goto 20 p=p*(-h*(k-1))/(h*(i-1)-h*(k-1)) 20 enddo s=s+p/(h*(i-1)-h*(j-1)) 10 enddo d=d+s*f(xo+h*(i-1)) enddo print*,'The First Derivation=',d end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The First Derivation= 30.856790

50

Problem: Write a Fortran program to find second derivation of the function: f ( x) = ln(5x 2 ) + 8x 3 − 8 at x = 1.1 using derivation of Lagrange polynomial.

! Second Derivation using Lagrange Formula f(x)=log(5*x**2)+8*x**3-8 n=4 xo=1.1 h=.1 d=0 do i=1,n s=0 p=1 do j=1,n if (j.eq.i) goto 10 s=s+(-h*(j-1)) p=p*(i-j)*h 10 enddo d=d+2*s/p*f(xo+h*(i-1)) enddo print*,'The Second Derivation=',d end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The Second Derivation= 51.200130

51

Problem: Write a program in Fortran to solve the differential equation y′ = e −2 x − 2 y using Euler method over [0,2] , with y(0) = 0.1, let h = 0.1 (i.e. n = 20 ).

! Euler's Method dimension x(0:30), y(0:30) a=0 b=2 n=20 y(0)=0.1 x(0)=a h=(b-a)/n do i=0,n-1 x(i+1)=x(i)+h y(i+1)=y(i)+h*(exp(-2*x(i))-2*y(i)) enddo do i=0,n print*,x(i),y(i) enddo end

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

52

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2

1.000000 E-01 1.800000 E-01 2.258731 E-01 2.477305 E-01 2.530655 E-01 2.473853 E-01 2.346962 E-01 2.178764 E-01 1.989608 E-01 1.793583 E-01 1.600165 E-01 1.415467 E-01 1.243177 E-01 1.085260 E-01 9.424812E-02 8.147950 E-01 7.016230 E-01 6.020606 E-02 5.150217 E-02 4.393411E-02 3.738436 E-02

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ −2 x − 2 y is: Exact solution for y′ = e y=

1 −2 e 10

x

− x e −2

x

⇒

y(2) = 0.03846284

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

53

Problem: Write a program in Fortran to solve the differential equation y′ = e −2 x − 2 y using Runge-Kutta method over [0,2] , with y(0) = 0.1, let n = 20 .

! Runge Kutta Method real k1,k2,k3,k4 dimension x(0:25),y(0:25) a=0 b=2 n=20 x(0)=a y(0)=0.1 h=(b-a)/n do i=0,n-1 xo=x(i) yo=y(i) k1=h*(exp(-2*xo)-2*yo) xo=x(i)+h/2 yo=y(i)+k1/2 k2=h*(exp(-2*xo)-2*yo) xo=x(i)+h/2 yo=y(i)+k2/2 k3=h*(exp(-2*xo)-2*yo) xo=x(i)+h yo=y(i)+k3 k4=h*(exp(-2*xo)-2*yo) y(i+1)=y(i)+1./6*(k1+2*k2+2*k3+k4) x(i+1)=x(i)+h enddo do i=0,n print*,x(i),y(i) enddo end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 54

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2

1.000000 E-01 1.637440 E-01 2.010928 E-01 2.195209 E-01 2.246607 E-01 2.207241 E-01 2.108327 E-01 1.972747 E-01 1.817045 E-01 1.652969 E-01 1.488672 E-01 1.329626 E-01 1.179324 E-01 1.039823 E-01 9.121463 E-02 7.965902 E-02 6.929560 E-02 6.007185 E-02 5.191512 E-02 4.474165 E-02 3.846299 E-02

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Exact solution for y′ = e −2 x − 2 y is: y=

1 −2 e 10

x

− x e −2

x

⇒

y(2) = 0.03846284

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

55

Problem: Write a program in Fortran to solve the differential equation y' = − x y using Taylor series method for values xo = 0.0 , x1 = 0.1, x2 = 0.2 , x3 = 0.3 , x4 = 0.4 , and x5 = 0.5 , with the initial condition y(0) = 1.

! Taylor Series Method dimension x(0:25),y(0:25) x(0)=0 y(0)=1 h=.1 n=5 do i=0,n-1 d1=-1*x(i)*y(i) d2=(x(i)**2-1)*y(i) d3=(-1*x(i)**3+3*x(i))*y(i) d4=(x(i)**4-6*x(i)**2+3)*y(i) y(i+1)=y(i)+h*d1+h**2/2*d2+h**3/6*d3+h**4/24 x(i+1)=x(i)+h enddo do i=0,n print*,x(i),y(i) enddo end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 0 0.1 0.2 0.3 0.4 0.5

1.000000 9.950042 E-01 9.801826 E-01 9.559749 E-01 9.230893 E-01 8.824676 E-01

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ −

x2 2

y(0.5) = e −0.125 = 0.08824969 y=e ⇒ Exact solution is: ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 56