(* * * *  *  * * * *)
(*  Dawid Ejdeholm  *)
(* * *  *  * *  *  *)

(* * * *  *  *  *)
(* Seminare 1   *)
(* * * *  *  *  *)

(* * * *  *  *  *)
(* Uppgift  1.3 *)
(* * * *  *  *  *)
Remove["Global`*"];
h=Table[N[10^-i],{i,1,5}];
For[i=1,i<6,i++,
	y[x_] = (Sqrt[x]);
	ℓ[x_] = k*x+m;
	(Subscript[x, 1])=1;
	(Subscript[y, 1])=1;
	(Subscript[x, 2])=1-h[[i]];
	(Subscript[y, 2])=(Sqrt[1-h[[i]]]);
	km= Solve[{(Subscript[y, 1])==k*(Subscript[x, 1])+m,(Subscript[y, 2])==k*(Subscript[x, 2])+m},{k,m}];
	nl[x_] = ℓ[x]/.km[[1]];
	δ[x_] = y[x]-nl[x];
	dδ =Solve[D[δ[x],x]==0,x]; (* x* *)
	val = δ[x]/.dδ[[1]];
	Print["h: ",h[[i]], " - x*: ",dδ[[1]][[1]] ," - δ*: ",val]
]
Remove["Global`*"];
y[x_] = (Sqrt[x]);
ℓ[x_] = k*x+m;
(Subscript[x, 1])=1;
(Subscript[y, 1])=1;
(Subscript[x, 2])=1-h;
(Subscript[y, 2])=(Sqrt[1-h]);
km= Solve[{(Subscript[y, 1])==k*(Subscript[x, 1])+m,(Subscript[y, 2])==k*(Subscript[x, 2])+m},{k,m}];
nl[x_] = ℓ[x]/.km[[1]];
δ[x_] = y[x]-nl[x];
dδ =Solve[D[δ[x],x]==0,x];
val = δ[x]/.dδ[[1]];

taylor =Normal[Series[val,{h,0,2}]]

(* * * *  *  *  *)
(* Uppgift 2.3  *)
(* * * *  *  *  *)
Remove["Global`*"]
data =Transpose[(({{1, 2, 3, 4, 5}, {1, 3, 2, 4, 2}}))];
deg = 4;
p4[x_] =C[0]+
	C[1](x-data[[1]][[1]])+
	C[2](x-data[[1]][[1]])(x-data[[2]][[1]])+
	C[3](x-data[[1]][[1]])(x-data[[2]][[1]])(x-data[[3]][[1]])+
	C[4](x-data[[1]][[1]])(x-data[[2]][[1]])(x-data[[3]][[1]])(x-data[[4]][[1]]);

const = Table[C[i], {i, 0,deg}];
sol=Solve[Table[p4[data[[i]][[1]]]==data[[i]][[2]],{i,1,deg+1}],const];

p[x_]=Fit[data,{1,x,x^2,x^3,x^4},x];

Print["Our Vandermond formulation: ",Simplify[p4[x]/.sol[[1]]]]
Print["Fit function:               ",p[x]]
Print["Are they the same? ",Expand[p4[x]/.sol[[1]]]==Expand[p[x]]]

pl=Plot[p4[x]/.sol[[1]],{x,data[[1]][[1]],data[[deg+1]][[1]]},PlotLegends->{"p4(x)"}];
lp=ListPlot[data,PlotStyle->{Red,PointSize[Large]},PlotLegends->{"data"}];
Show[{pl, lp}]

Print["Compare with the result given by Fit. Estimate y(4.2) and y'(4.2). Compare with y(4.2) from 2.2"]
Print["(From 2.2) Interpolation[data, 4.2]=",Interpolation[data, 4.2] ](*from 2.2*)
Print["Estimated y(4.2) with Fit. p[4.2]=",p[4.2]] (*Estimated y(4.2) *)
Print["Estimated y(4.2) with our function. p4[4.2]=",p4[4.2] /.sol[[1]]](* Estmited y(4.2) *)
p4d[x_] = D[p4[x],x]/.sol[[1]]; (* cool derivate function *)
Print["Estimated y'(4.2) with our p4d[4.2]=",p4d[4.2]]

(* * * *  *  *  *)
(* Uppgift 2.5 *)
(* * * *  *  *  *)
Remove["Global`*"]
xv={12.8,17.6, 21.5,24.0, 27.1, 30.2, 32.5};
fv={2.4,0.6, 4.0, 6.0, 4.7, 1.8,-0.6};
data=Transpose[{xv,fv}];

n=Length[xv];
{dfl,dfr}={0,0};

deg=3; (* s(x) spline-funktion av grad deg *)
k = (fv[[1]]-4.2)/(xv[[1]]-7.4);
m = 4.2-k*7.4;

l[x_] = k*x + m;
s1[x_]=Sum[a[i](x-xv[[n-6]])^i,{i,0,deg}];
s2[x_]=Sum[b[i](x-xv[[n-5]])^i,{i,0,deg}];
s3[x_]=Sum[c[i](x-xv[[n-4]])^i,{i,0,deg}];
s4[x_]=Sum[d[i](x-xv[[n-3]])^i,{i,0,deg}];
s5[x_]=Sum[e[i](x-xv[[n-2]])^i,{i,0,deg}];
s6[x_]=Sum[f[i](x-xv[[n-1]])^i,{i,0,deg}];


const=Join[
Table[a[i],{i,0,deg}],Table[b[i],{i,0,deg}],
Table[c[i],{i,0,deg}],Table[d[i],{i,0,deg}],
Table[e[i],{i,0,deg}],Table[f[i],{i,0,deg}]];

sol=Solve[{
s1[xv[[1]]] == fv[[1]],
s1[xv[[2]]]== s2[xv[[2]]] == fv[[2]],
s2[xv[[3]]] == s3[xv[[3]]] ==  fv[[3]],
s3[xv[[4]]] == s4[xv[[4]]] ==  fv[[4]],
s4[xv[[5]]] == s5[xv[[5]]] ==  fv[[5]],
s5[xv[[6]]] == s6[xv[[6]]] ==  fv[[6]],
s6[xv[[7]]] == fv[[7]],

s1'[xv[[2]]] == s2'[xv[[2]]],
s2'[xv[[3]]] == s3'[xv[[3]]],
s3'[xv[[4]]] == s4'[xv[[4]]],
s4'[xv[[5]]] == s5'[xv[[5]]],
s5'[xv[[6]]] == s6'[xv[[6]]],

s1''[xv[[2]]]==s2''[xv[[2]]],
s2''[xv[[3]]]==s3''[xv[[3]]],
s3''[xv[[4]]]==s4''[xv[[4]]],
s4''[xv[[5]]]==s5''[xv[[5]]],
s5''[xv[[6]]]==s6''[xv[[6]]],

s1'[xv[[1]]] == k,
s6''[xv[[7]]] == 0
},const
];

s[x_] = l[x]UnitBox[((x-7.4)/(xv[[1]]-7.4))-(1/2)]+
  s1[x]UnitBox[((x-xv[[1]])/(xv[[2]]-xv[[1]]))-(1/2)] + s2[x]UnitBox[((x-xv[[2]])/(xv[[3]]-xv[[2]]))-(1/2)] + 
s3[x]UnitBox[((x-xv[[3]])/(xv[[4]]-xv[[3]]))-(1/2)] + s4[x]UnitBox[((x-xv[[4]])/(xv[[5]]-xv[[4]]))-(1/2)] +
s5[x]UnitBox[((x-xv[[5]])/(xv[[6]]-xv[[5]]))-(1/2)] +
s6[x]UnitBox[((x-xv[[6]])/(xv[[7]]-xv[[6]]))-(1/2)];
pl = Plot[s[x]/.sol[[1]], {x, 7.4, xv[[7]]}];
lp=ListPlot[data,PlotStyle->{Red,PointSize[Large]}];
Show[{pl, lp}]

(* * * *  *  *  *)
(* Uppgift 3.1 *)
(* A, B, C      *)
(* * * *  *  *  *)
Remove["Global`*"]
xv={1,2,3.5,5.8};
yv={1.1,1.8,4,5.5};

data=Transpose[Append[{xv},yv]];
𝔸={#^2,#,1}&/@data[[All,1]];
p ={a,b,c};
fx={x^2,x,1}.p;
𝔸T=Transpose[𝔸];
sol = Solve[𝔸T.𝔸.p == 𝔸T.yv,p];
answer=fx/.sol;
Print["Min lösning:      ",answer[[1]]]
Print["Mathematicas Fit: ",Fit[data,{1,x,x^2},x]]

pl=Plot[answer,{x,xv[[1]],xv[[4]]}];
lp=ListPlot[data,PlotStyle->{Red,PointSize[Large]}];
Show[lp,pl]

(* * * *  *  *  *)
(* Uppgift 3.4 *)
(* * * *  *  *  *)
Remove["Global`*"]
xv={1,2.5,4,6,8};
yv={1,1.5,2,2.5,3};
data=Transpose[Append[{xv},yv]];
y = (x/((a*x+b)));
z=x/y;
zv =xv/yv;
𝔸={#,1}&/@data[[All,1]];
p={a,b};
𝔸T =Transpose[𝔸];
sol = Solve[𝔸T.𝔸.p == 𝔸T.zv,{a,b}];
Print["Min lösning:      ",FullSimplify[sol[[1]]]];
pf = FindFit[data,(x/(a*x+b)),{a,b},x];
Print["Mathematicas Fit: ",pf]
answer = y/.sol;
pl=Plot[answer,{x,xv[[1]],xv[[5]]},PlotLegends->{"Min"}];
pfp = Plot[y/.pf,{x,xv[[1]],xv[[5]]} ,PlotStyle->{Green,Dashed}, PlotLegends->{"Fit"}];
lp=ListPlot[data,PlotStyle->{Red,PointSize[Large]}];
Show[{lp,pl, pfp}]

(* * * *  *  *  *)
(* Uppgift 2.6 *)
(* * * *  *  *  *)
Print["*******"]
Print["* 2.6 *"]
Print["*******"]
yn=Table[s[i],{i,7.4,32.5}]/.sol[[1]];
xn=Table[i,{i,7.4,32.5}];
data2=Transpose[Append[{xn},yn]];
lp2=ListPlot[data2,PlotStyle->{Red,PointSize[Large]}];
Show[lp2]

(* * * *  *  *  *)
(* Uppgift 4.5 *)
(* * * *  *  *  *)
Remove["Global`*"]
xv={1.4,1.5,1.7};
fv={0.98545,0.99749,0.99167};
data=Transpose[Append[{xv},fv]];
c=Table[C[i],{i,0,2}];
p=Table[x^i,{i,0,2}].c;
a=Transpose[Table[Table[xv[[i]]^j,{i,1,3}],{j,0,2}]];
sol=Solve[a.c==Table[fv[[i]],{i,1,3}],c];
Print["(Subscript[p, 2])(x) = ",pp=p/.sol[[1]]," | (Subscript[p, 2])'(",π/2,") = " , D[pp,x]/.x->π/2]
Print["M Fit = ",Fit[{data[[1]],data[[2]],data[[3]]},{1,x,x^2},x]]

(* * * *  *  *  *)
(* Uppgift 4.7 *)
(* * * *  *  *  *)
Remove["Global`*"]
f[x_]=Sin[x];
(* f(x) =sin(x). Estimate f'(pi/2) med back, forward och center för h=10^-k, k=1,2,3.. tills datorn dör. Jämför f'[pi/2]=cos(pi/2)=0*)
forward[x_,h_]:=((f[x+h]-f[x])/h);
backward[x_,h_]:=((f[x]-f[x-h])/h);
center[x_,h_]:=((f[x+h]-f[x-h])/(2h));
k=1;
Print["------------"]
While[forward[Pi/2,N[10^-k]]≠0,Print["Forward[Pi/2, h=",N[10^-k],"]=",forward[Pi/2,N[10^-k]]];k++];

Print["Forward[Pi/2, h=",N[10^-k],"]=",forward[Pi/2,N[10^-k]]]
Print["Med k=",k," får vi 0"]
Print["------------"]
k=1;
While[backward[Pi/2,N[10^-k]]≠0,Print["Backward[Pi/2, h=",N[10^-k],"]=",backward[Pi/2,N[10^-k]]];k++];
Print["Backward[Pi/2, h=",N[10^-k],"]=",backward[Pi/2,N[10^-k]]]
Print["Med k=",k," får vi 0"]
Print["------------"]
k=1;
While[center[Pi/2,N[10^-k]]≠0;k++];Print["Center[Pi/2, h=",N[10^-k],"]=",center[Pi/2,N[10^-k]]]
Print["Med k=",k," får vi 0"]

(* * * *  *  *  *)
(* Seminare 2   *)
(* * * *  *  *  *)

(* * * *  *  *  *)
(* Uppgift 5.7  *)
(* * * *  *  *  *)
Remove["Global`*"]
f[x_]=Log[Cos[x]]; 
{a,b}={0,1};
n = 50;
h=(b-a)/n;
s = f[a];
For[i = 1, i≤n-1,i=i+1,
xi =a+h*i;
s=s+2*f[xi];
]
s=s+f[b];
N[h*s/2,4]
N[NIntegrate[Log[Cos[x]],{x,0,1}],4]

(* * * *  *  *  *)
(* Uppgift 6.2  *)
(* * * *  *  *  *)
Remove["Global`*"]
f[x_]:=x^3+x-1;
Plot[f[x],{x,0.5,0.8}]
{a,b}={0,1};
epsilon = 0.001;
While[b-a>epsilon,
c=0.5*(a+b);
Print[{a,c,b}];
If[f[a]*f[c]<0,b=c,a=c]
]
NSolve[f[x]==0,x][[3]]

(* * * *  *  *  *)
(* Uppgift 6.5  *)
(* * * *  *  *  *)
Remove["Global`*"]
f[x_]:=x*Sin[x]-4*Cos[x]+Exp[x];
newton[x_] := x - f[x]/f'[x];
epsilon = 0.0001;
diff = 1;
xold = -4;
While[diff>epsilon,
xnew=newton[xold];
diff=Abs[xnew-xold];
xold=xnew;
]
Print[{"x "-> N[xnew]}];
FindRoot[f[x]==0,{x,-4}]

(* * * *  *  *  *)
(* Uppgift 6.9  *)
(* * * *  *  *  *)
Remove["Global`*"]
a=2;
f[x_]=x^2-a;
g[x_]={0.5(x+(a/x)),0.5(x+(a/x))};
list={{"k","g1", "g2"}};
kmax=3;
xold={-1,1};
For[k=1,k≤kmax,++k,
For[i=1,i≤2,i++,
xnew[i]=g[xold[[i]]][[i]];
xold[[i]]=xnew[i];
];
AppendTo[list,Flatten[{k,Table[g[xnew[j]][[j]],{j,1,2}]}]]
]
Grid[list,Frame->All]

(* * * *  *  *  *)
(* Uppgift 6.10  *)
(* * * *  *  *  *)
(* Newton heat/cooling *)
Remove["Global`*"]
T1[t_]=100(1-Exp[-0.2*t]);
T2[t_]=40*Exp[-0.01*t];
f[x_]=T1[x]-T2[x];
newtons[x_]:=x-((f[x])/(f'[x]));
g[x_]=newtons[x];
kmax=10;
xold=1;
For[k=1,k≤kmax,++k,
xnew=g[xold];
xold=xnew;
];

Print["t* = ", xnew]
Print["T1(t*) = ",  T1[xnew]]
Print["T2(t*) = ",  T2[xnew]]
Print ["FindRoot(...) = ",FindRoot[T1[x]-T2[x]==0,{x,0}]]