Quit[];

(* Esercizio 5.2.11 *)
f[t_,y_]:=r*b*(1-y)
eq0=p'[t]==f[t,p[t]];

α=0.01;b0=0.02;d0=0.015;r0=0.1;
tA=0.0; tB=2000.;
n=40;
h=(tB-tA)/n;
t[0]=tA;
t[n]=tB;
Do[t[i+1]=t[i]+h,{i,0,n-1}]
w[0]=α
Do[w[i+1]=w[i]+h*(f[t[i],w[i]]/.{b->b0,d->d0,r->r0}),{i,0,n-1}]
sol=Table[{t[i],w[i]},{i,0,n}];
soltrue=DSolve[{eq0,p[tA]==α},p[t],t][[1]]

trange={tA,tB};
yrange={0,1};
pl1=ListPlot[sol,PlotRange->{trange,yrange}];
(*pl2=Plot[y[t]/.solnum,{t,tA,tB},PlotRange->{trange,yrange}];*)
pl3=Plot[p[t]/.soltrue/.{b->b0,d->d0,r->r0},{t,tA,tB},PlotRange->{trange,yrange},PlotStyle->RGBColor[1,0,0]]
Show[pl1,pl3]

(* Esercizio 5.4.15 con RK4 *)
tA=0.0; tB=0.4;
n=100;
h=(tB-tA)/n;
t[0]=tA;
t[n]=tB;
Do[t[i+1]=t[i]+h,{i,0,n-1}]
n1=1000;
n2=2000.;
n3=3000.;
α=n3;
w[0]=α;
k=6.22*10^(-19);
f[t_,y_]:=k*(n1-y/2)^2*(n2-y/2)^2*(n3-3*y/4)^3
Do[k1=f[t[i],w[i]];
k2=f[t[i]+0.5*h,w[i]+0.5*h*k1];
k3=f[t[i]+0.5*h,w[i]+0.5*h*k2];
k4=f[t[i]+h,w[i]+h*k3];
w[i+1]=w[i]+h*(k1+2*k2+2*k3+k4)/6.,{i,0,n-1}];
sol=Table[{t[i],w[i]},{i,0,n}];
solnum=NDSolve[{y'[t]==f[t,y[t]],y[tA]==α},y[t],{t,tA,tB}][[1]];
(*soltrue=DSolve[{y'[t]==f[t,y[t]],y[tA]==α},y[t],t][[1]]*)
trange={tA,tB};
yrange={0,4000.};
pl1=ListPlot[sol,PlotRange->{trange,yrange},PlotStyle->{PointSize[0.015],RGBColor[1,0,0]}];
pl3=Plot[y[t]/.solnum,{t,tA,tB},PlotRange->{trange,yrange},PlotStyle->Thickness[0.005]];
(*pl3=Plot[y[t]/.soltrue,{t,tA,tB},PlotRange->{trange,yrange},PlotStyle->RGBColor[1,0,0]];*)
Show[pl1,pl3]

(* Esercizio 6.2.16 *)
Clear[a,α];
n=3;
a={{2,1,3},{4,6,8},{6,α,10}};
b={1,5,5};
MatrixForm[a]
MatrixForm[b]

a0=a/.α->6;
MatrixForm[a0]
Do[s[i]=Max[Table[Abs[a0[[i]][[j]]]//N,{j,1,n}]],{i,1,n}];
Do[r[i]=Abs[a0[[i]][[1]]]/s[i]//N,{i,1,n}];
Do[Print["i = ",i, "    s[i] = ",s[i], "   r[i] = ", r[i]],{i,1,n}]
aa={Join[a0[[1]],{b[[1]]}],Join[a0[[2]],{b[[2]]}],Join[a0[[3]],{b[[3]]}]};
MatrixForm[aa]

a1={aa[[1]],aa[[2]]-(aa[[2]][[1]]/aa[[1]][[1]])*aa[[1]],aa[[3]]-(aa[[3]][[1]]/aa[[1]][[1]])*aa[[1]]};
MatrixForm[a1]
Do[s[i]=Max[Table[Abs[a1[[i]][[j]]]//N,{j,2,n}]],{i,2,n}];
Do[r[i]=Abs[a1[[i]][[2]]]/s[i]//N,{i,2,n}];
Do[Print["i = ",i, "    s[i] = ",s[i], "   r[i] = ", r[i]],{i,2,n}]

a2={a1[[1]],a1[[2]],a1[[3]]-(a1[[3]][[2]]/a1[[2]][[2]])*a1[[2]]};
MatrixForm[a2]

(* Esercizio 5.2.12 *)
Quit[];

g[t_]=Exp[-0.06*Pi*t]*Sin[2.0*t-Pi];

c=0.3; r=1.4; L=1.7; 

rhs[t_]=c*g''[t]+g'[t]/r+g[t]/L;
eq=f'[t]==rhs[t];

sol=DSolve[{eq,f[0]==0},f[t],t][[1]]

pl1=Plot[f[t]/.sol,{t,0,10},PlotStyle->RGBColor[0,0,0]]

n=11;
tf=10.0;
dt=tf/(n-1);
tt[1]=0.0;
Do[tt[i+1]=tt[i]+dt,{i,1,n-1}];
f0=0.0;
w[1]=f0;
Do[w[i+1]=w[i]+dt*rhs[tt[i]],{i,1,n-1}]
feul=Table[{tt[i],w[i]},{i,1,n}];

pl2=ListPlot[feul,PlotStyle->{RGBColor[1,0,0],PointSize[0.02]}]

Show[pl1,pl2]

(* CON RK2 *)
Clear[w];
w[1]=f0;
Do[k1=rhs[tt[i]];k2=rhs[tt[i]+0.5*dt];w[i+1]=w[i]+dt*k2,{i,1,n-1}]
frk2=Table[{tt[i],w[i]},{i,1,n}];

pl3=ListPlot[frk2,PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}]

(* CON RK4 *)
Clear[w];
w[1]=f0;
Do[k1=rhs[tt[i]];k2=rhs[tt[i]+0.5*dt];k3=rhs[tt[i]+0.5*dt];k4=rhs[tt[i]+dt];w[i+1]=w[i]+dt*(k1+2.0*k2+2.0*k3+k4)/6.0,{i,1,n-1}]
frk4=Table[{tt[i],w[i]},{i,1,n}];

pl4=ListPlot[frk4,PlotStyle->{RGBColor[0,1,0],PointSize[0.02]}]

Show[pl1,pl2,pl3,pl4]

(* MODELLO MALTHUSIANO *)
Clear[f,y]
b=0.1;d=0.05;
tf=150;
f[t_,y_]=(b-d)*y
eq=y'[t]==f[t,y[t]];
sol=DSolve[{eq,y[0]==1},y[t],t][[1]]
pl1=Plot[y[t]/.sol,{t,0,tf},PlotRange->{0,2000},PlotStyle->{Thickness[0.01],RGBColor[0,0,0]}]

n=21;
dt=tf/(n-1);
tt[1]=0.0;
Do[tt[i+1]=tt[i]+dt,{i,1,n-1}];
f0=1.0;
w[1]=f0;
Do[w[i+1]=w[i]+dt*f[tt[i],w[i]],{i,1,n-1}]
feul=Table[{tt[i],w[i]},{i,1,n}];

pl2=ListPlot[feul,PlotStyle->{RGBColor[1,0,0],PointSize[0.02]}]

Show[pl1,pl2]

(* CON RK2 *)
Clear[w];
w[1]=f0;
Do[k1=f[tt[i],w[i]];k2=f[tt[i]+0.5*dt,w[i]+0.5*dt*k1];w[i+1]=w[i]+dt*k2,{i,1,n-1}]
frk2=Table[{tt[i],w[i]},{i,1,n}];

pl3=ListPlot[frk2,PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}]

Show[pl1,pl2,pl3]

(* CON RK4 *)
Clear[w];
w[1]=f0;
Do[k1=f[tt[i],w[i]];k2=f[tt[i]+0.5*dt,w[i]+0.5*dt*k1];k3=f[tt[i]+0.5*dt,w[i]+0.5*dt*k2];k4=f[tt[i]+dt,w[i]+dt*k3];w[i+1]=w[i]+dt*(k1+2.0*k2+2.0*k3+k4)/6.0,{i,1,n-1}]
frk4=Table[{tt[i],w[i]},{i,1,n}];

pl4=ListPlot[frk4,PlotStyle->{RGBColor[0,1,0],PointSize[0.02]}]

Show[pl1,pl2,pl3,pl4]

(* MODELLO LOGISTICO *)
Clear[f,y]
a=0.1;b=0.01;
tf=150;
f[t_,y_]=a*y-b*y^2
eq=y'[t]==f[t,y[t]];
sol=DSolve[{eq,y[0]==0.1},y[t],t][[1]]
pl1=Plot[y[t]/.sol,{t,0,tf},PlotRange->{0,12},PlotStyle->{Thickness[0.01],RGBColor[0,0,0]}]

n=21;
dt=tf/(n-1);
tt[1]=0.0;
Do[tt[i+1]=tt[i]+dt,{i,1,n-1}];
f0=0.1;
w[1]=f0;
Do[w[i+1]=w[i]+dt*f[tt[i],w[i]],{i,1,n-1}]
feul=Table[{tt[i],w[i]},{i,1,n}];

pl2=ListPlot[feul,PlotStyle->{RGBColor[1,0,0],PointSize[0.02]}]

Show[pl1,pl2]

(* CON RK2 *)
Clear[w];
w[1]=f0;
Do[k1=f[tt[i],w[i]];k2=f[tt[i]+0.5*dt,w[i]+0.5*dt*k1];w[i+1]=w[i]+dt*k2,{i,1,n-1}]
frk2=Table[{tt[i],w[i]},{i,1,n}];

pl3=ListPlot[frk2,PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}]

Show[pl1,pl2,pl3]

(* CON RK4 *)
Clear[w];
w[1]=f0;
Do[k1=f[tt[i],w[i]];k2=f[tt[i]+0.5*dt,w[i]+0.5*dt*k1];k3=f[tt[i]+0.5*dt,w[i]+0.5*dt*k2];k4=f[tt[i]+dt,w[i]+dt*k3];w[i+1]=w[i]+dt*(k1+2.0*k2+2.0*k3+k4)/6.0,{i,1,n-1}]
frk4=Table[{tt[i],w[i]},{i,1,n}];

pl4=ListPlot[frk4,PlotStyle->{RGBColor[0,1,0],PointSize[0.02]}]

Show[pl1,pl2,pl3,pl4]