(** ESERCIZIO 4.1.2,3 **)

f[x_]:=Sin[x];
x[0]=0.5;x[1]=0.6; x[2]=0.7;h=0.1;
y[0]=0.4794;y[1]=0.5646; y[2]=0.6442;
Plot[f''[x],{x,0,0.7}]
Print["Derivate"]
dy[0]=(y[1]-y[0])/(x[1]-x[0])
dyf[1]=(y[2]-y[1])/(x[2]-x[1])
dyb[1]=(y[1]-y[0])/(x[1]-x[0])
dy[2]=(y[2]-y[1])/(x[2]-x[1])

Print["Errori"]
e[0]=Abs[dy[0]-f'[x[0]]]
ef[1]=Abs[dyf[1]-f'[x[1]]]
eb[1]=Abs[dyb[1]-f'[x[1]]]
e[2]=Abs[dy[2]-f'[x[2]]]
Print["Err. bound"]
b[0]=Abs[h/2]*Abs[f''[x[1]]]
b[1]=Abs[h/2]*Abs[f''[x[2]]]
b[2]=Abs[h/2]*Abs[f''[x[2]]]

f[x_]:=Exp[x]-2x^2+3x-1;
Plot[Abs[f''[x]],{x,0,0.4}]

x[0]=0.0;x[1]=0.2; x[2]=0.4;h=0.2;
y[0]=f[x[0]]
y[1]=f[x[1]]
y[2]=f[x[2]]

Print["Derivate"]
dy[0]=(y[1]-y[0])/(x[1]-x[0])
dyf[1]=(y[2]-y[1])/(x[2]-x[1])
dyb[1]=(y[1]-y[0])/(x[1]-x[0])
dy[2]=(y[2]-y[1])/(x[2]-x[1])
Print["Errori"]
e[0]=Abs[dy[0]-f'[x[0]]]
ef[1]=Abs[dyf[1]-f'[x[1]]]
eb[1]=Abs[dyb[1]-f'[x[1]]]
e[2]=Abs[dy[2]-f'[x[2]]]
Print["Err. bound"]
b[0]=Abs[h/2]*Abs[f''[x[0]]]
b[1]=Abs[h/2]*Abs[f''[x[1]]]
b[2]=Abs[h/2]*Abs[f''[x[1]]]

(** ESERCIZIO 4.4.1 e 2 **)

(* (a) *)
f[x_]:=x*Log[x];
a=1.; b=2.;
true1=Integrate[f[x],x]
true=Integrate[f[x],{x,a,b}]//N

x[0]=a;
n=4;
x[n]=b;
h=(b-a)/n;
Do[x[i+1]=x[i]+h,{i,0,n-2}];
sum=0.0;
Do[sum+=f[x[i]],{i,1,n-1}];
trap=0.5*h*(f[a]+2.0*sum+f[b]);
err=Abs[trap-true];
Print["true = ",true];
Print["trap = ", trap];
Print["error = ", err];

w[0]=1.;w[n]=1.;
Do[w[i]=4.,{i,1,n-1,2}]
Do[w[i]=2.,{i,2,n-2,2}]
sum=0.0;
Do[sum+=w[i]*f[x[i]],{i,0,n}];
simps=h*sum/3.;
err=Abs[simps-true];
Print["true = ",true];
Print["simps = ", simps];
Print["error = ", err];
Plot[f[x],{x,a,b}]

(* (e) *)
Clear[x,f];
f[x_]:=Exp[2x]*Sin[3x];
a=0.; b=2.;
true1=Integrate[f[x],x]
true=Integrate[f[x],{x,a,b}]//N

x[0]=a;
n=8;
x[n]=b;
h=(b-a)/n;
Do[x[i+1]=x[i]+h,{i,0,n-2}];
sum=0.0;
Do[sum+=f[x[i]],{i,1,n-1}];
trap=0.5*h*(f[a]+2.0*sum+f[b]);
err=Abs[trap-true];
Print["true = ",true];
Print["trap = ", trap];
Print["error = ", err];

w[0]=1.;w[n]=1.;
Do[w[i]=4.,{i,1,n-1,2}]
Do[w[i]=2.,{i,2,n-2,2}]
sum=0.0;
Do[sum+=w[i]*f[x[i]],{i,0,n}];
simps=h*sum/3.;
err=Abs[simps-true];
Print["true = ",true];
Print["simps = ", simps];
Print["error = ", err];
Plot[f[x],{x,a,b}]

(** ESERCIZIO 4.4.7 **)

Clear[x,f];
f[x_]:=Exp[2x]*Sin[3x];
a=0.; b=2.;
true1=Integrate[f[x],x];
true=Integrate[f[x],{x,a,b}]//N

(* Stima teorica *)
Simplify[f''[x]]
Plot[f''[x],{x,a,b}]
Simplify[f''''[x]]
Plot[f''''[x],{x,a,b}]

(* Trapezi *)
Clear[h];
tol=0.0001;
errbound=(b-a)/12*h^2*f''[b]
hbound=Sqrt[12.*tol/((b-a)*f''[b])]
nbound=(b-a)/hbound

(* Simpson *)
Clear[h];
tol=0.0001;
errbound=Abs[(b-a)/180*h^4*f''''[b]]
hbound=Sqrt[Sqrt[Abs[180.*tol/((b-a)*f''''[b])]]]
nbound=(b-a)/hbound

(* Calcolo numerico *)
(* Trapezi *)
x[0]=a;
n=800;
x[n]=b;
h=(b-a)/n;
Do[x[i+1]=x[i]+h,{i,0,n-2}];
sum=0.0;
Do[sum+=f[x[i]],{i,1,n-1}];
trap=0.5*h*(f[a]+2.0*sum+f[b]);
err=Abs[trap-true];
Print["true = ",true];
Print["trap = ", trap];
Print["error = ", err];

(* Simpson *)
x[0]=a;
n=34;
x[n]=b;
h=(b-a)/n;
Do[x[i+1]=x[i]+h,{i,0,n-2}];
w[0]=1.;w[n]=1.;
Do[w[i]=4.,{i,1,n-1,2}]
Do[w[i]=2.,{i,2,n-2,2}]
sum=0.0;
Do[sum+=w[i]*f[x[i]],{i,0,n}];
simps=h*sum/3.;
err=Abs[simps-true];
Print["true = ",true];
Print["simps = ", simps];
Print["error = ", err];

(** ESERCIZIO 4.4.8 **)

Clear[x,f];
f[x_]:=x^2*Cos[x];
a=0.; b=Pi;
true1=Integrate[f[x],x];
true=Integrate[f[x],{x,a,b}]//N

(* Stima teorica *)
Simplify[f''[x]]
Plot[f''[x],{x,a,b}]
Simplify[f''''[x]]
Plot[f''''[x],{x,a,b}]

(* Trapezi *)
Clear[h];
tol=0.0001;
errbound=(b-a)/12*h^2*f''[b]
hbound=Sqrt[12.*tol/((b-a)*f''[b])]
nbound=(b-a)/hbound

(* Simpson *)
Clear[h];
tol=0.0001;
errbound=Abs[(b-a)/180*h^4*f''''[b]]
hbound=Sqrt[Sqrt[Abs[180.*tol/((b-a)*f''''[b])]]]
nbound=(b-a)/hbound

(* Calcolo numerico *)
(* Trapezi *)
x[0]=a;
n=400;
x[n]=b;
h=(b-a)/n;
Do[x[i+1]=x[i]+h,{i,0,n-2}];
sum=0.0;
Do[sum+=f[x[i]],{i,1,n-1}];
trap=0.5*h*(f[a]+2.0*sum+f[b]);
err=Abs[trap-true];
Print["true = ",true];
Print["trap = ", trap];
Print["error = ", err];

(* Simpson *)
x[0]=a;
n=20;
x[n]=b;
h=(b-a)/n;
Do[x[i+1]=x[i]+h,{i,0,n-2}];
w[0]=1.;w[n]=1.;
Do[w[i]=4.,{i,1,n-1,2}]
Do[w[i]=2.,{i,2,n-2,2}]
sum=0.0;
Do[sum+=w[i]*f[x[i]],{i,0,n}];
simps=h*sum/3.;
err=Abs[simps-true];
Print["true = ",true];
Print["simps = ", simps];
Print["error = ", err];