Rjesavanje trodijagonalnog sustava LU - dekompozicijom

(* b je gornja, a c donja dijagonala *)
LUtri[n_,a_,b_,c_,y_]:=Module[{d,e,f,z},
z=Table[0,{i,n}]; s=Table[0,{i,n}]; e=Table[0,{i,n}]; f=Table[0,{i,n}];d=Table[0,{i,n}]; 
(* Izracunavanje *)
e[[1]]=a[[1]];f[[1]]=b[[1]];d[[1]]=c[[1]]/e[[1]]; Do[e[[i]]=a[[i]]-d[[i-1]] f[[i-1]];f[[i]]=b[[i]]; d[[i]]=c[[i]]/e[[i]],{i,2,n-1}];
e[[n]]=a[[n]]-d[[n-1]] f[[n-1]];

(* Rjesavanje sustava LUx=y *)
z[[1]]=y[[1]]; 
Do[z[[i]]=y[[i]]-d[[i-1]]z[[i-1]],{i,2,n}];
s[[n]]=z[[n]]/e[[n]];
Do[s[[n-i-1]]=(z[[n-i-1]]-f[[n-i-1]] s[[n-i]] )/e[[n-i-1]],{i,0,n-2}];
s//N]

Interpolacija preko sustava jednadzbi

(* Definiranje podataka *)
f[x_]:=1-Abs[x-1]; a=0; b=2; n=10;
(*f[x_]:=Abs[Sin[x]]; a=0; b=2Pi;*)
 x=Table[(b-a)i/n,{i,0,n}];                  y=f[x];
(* x=Table[(a+b)/2+N[.5(b-a)Cos[(n-k)Pi/n]],{k,0,n}];   y=f[x]; *)    (* nultocke Ceb. pol *)
pod=Table[{x[[i]],y[[i]]},{i,n+1}];
Print[MatrixForm[Transpose[N[pod,2]]]]
                 (* Matrica *)              
mat=Table[0,{i,0,n},{j,0,n}];
Do[
Do[
If[j<2,mat[[i,1]]=1,mat[[i,j]]=N[x[[i]]^(j-1)]]
,{i,n+1}]
       ,{j,n+1}]; 
		                        (* Interpolacijski polinom *)
koef=LinearSolve[mat,y];
intpol[xx_]:=Sum[koef[[i]]xx^(i-1),{i,n+1}];
           (* Slika funkcije f interp. polinoma p i podataka *)
slf=Plot[f[x],{x,x[[1]],x[[n+1]]}];
slpod=ListPlot[pod,PlotStyle->AbsolutePointSize[4]];
slp=Plot[intpol[xx],{xx,x[[1]],x[[n+1]]}];
slf=Show[slf,slpod];
sl2=Show[slp,slpod,slf];
GraphicsGrid[{{slf,sl2}}]
		
cond=Norm[mat,Infinity]Norm[Inverse[mat],Infinity];
Print["Uvjetovanost =",cond]

Funkcija omega

omega[xx_]:=Product[(xx-x[[i]]),{i,n+1}]
Plot[omega[xx],{xx,x[[1]],x[[n+1]]},PlotRange->{-.01,.01}]

Pogreska

err[f_,p_,m_]:=Max[Table[Abs[f[a+i (b-a)/m]-intpol[a+i (b-a)/m]],{i,(b-a)/m,m}]]//N
err[f,intpol,100]

Lagrangeov oblik interpolacijskog polinoma

intpol[xx_]:=Sum[y[[i]] Apply[Times,Drop[Table[xx-x[[k]],{k,n+1}],{i}]] /Apply[Times,Drop[Table[x[[i]]-x[[k]],{k,n+1}],{i}]],{i,n+1}]

slp=Plot[intpol[xx],{xx,x[[1]],x[[n+1]]}];
sl2=Show[slp,slpod];
GraphicsGrid[{{slf,sl2}}]

err[f,intpol,100]

Newtonov oblik interpolacijskog polinoma

(* Podjeljene razlike *)
n=4; x={-2,-1,0,1};y={-4,-6,-2,2};
a=y;t=a; Do[kk=k-1;
	        Do[ ii=n-i-kk;
	  	   t[[ii]]=(y[[ii+1]]-y[[ii]])/(x[[ii+kk+1]]-x[[ii]]),{i,n-1}];
                                      	a[[k+1]]=t[[1]];y=t,{k,n-1}];a

Vrijednost polinoma u tocki z - Hornerova shema

p=a[[n]];z=1;
Do[ii=n-i;	p=p(z-x[[ii]])+a[[ii]],{i,n-1}];p

Podjeljene razlike

n=5; x={0,.25,.5,.75,1};y=Sin[Pi x];z=1;
a=y;t=a; Do[kk=k-1;
	        Do[ ii=n-i-kk;
	  	   t[[ii]]=(y[[ii+1]]-y[[ii]])/(x[[ii+kk+1]]-x[[ii]]),{i,n-1}];
                                      	a[[k+1]]=t[[1]];y=t,{k,n-1}];a

Vrijednost polinoma u tocki z - Hornerova shema

z=.25;p=a[[n]];
Do[ii=n-i;	p=p(z-x[[ii]])+a[[ii]],{i,n-1}];p

Vrijednost polinoma P (x) = an x^n + ... + a1 x + a0 
u tocki z - Hornerova shema

n=4;a={-1,0,5,-2}; z=1;
p=a[[n]];
Do[p=p z+a[[n-i]],{i,n-1}];p

Interpolacijski spline

Do[
spline=Interpolation[pod,InterpolationOrder->k];
sl[k]=Plot[spline[xx],{xx,x[[1]],x[[n+1]]}];
s[k]=Show[sl[k],slpod,PlotLabel->k];,
	{k,3}];
GraphicsGrid[{{slf,s[1]},{s[2],s[3]}}]

Spline interpolacija  - PROGRAM prema skripti

(* Definiranje podataka preko funkcije *)
(* m=10;SeedRandom[17];
f[t_]:=Sin[x]+1;
slf=Plot[f[t],{t,0,m}, DisplayFunction->Identity];
x=Table[i,{i,m}];
y=Table[f[i]+Random[NormalDistribution[0,.05]],{i,m}]; *)

(* Definiranje podataka direktno *)
n=5;
x={0,.25,.5,.75,1};
y={1,2,1,0,1};
h=Table[x[[i+1]]-x[[i]],{i,n-1}];
d=Table[(y[[i+1]]-y[[i]])/h[[i]],{i,n-1}];
pod=Table[{x[[i]],y[[i]]},{i,n}];
Print["podaci: ",TableForm[pod]]
slpod=ListPlot[pod,PlotStyle->{AbsolutePointSize[5]},
		PlotRange->{0,Max[y]}];

(* Izracunavanje *)
aa=Table[2(h[[i]]+h[[i+1]]),{i,n-2}];
bb=Table[h[[i]],{i,2,n-2}]; cc=bb;
yy=Table[6(d[[i+1]]-d[[i]]),{i,n-2}];
LUtri[n-2,aa,bb,cc,yy];
s=Prepend[s,0]; s=Append[s,0];

(* Definiranje kubicnih funkcija *)
 b=Table[d[[i]]-(s[[i+1]]-s[[i]])  h[[i]]/ 6,{i,n-1}];
c[xx_,i_]:=(y[[i-1]]-s[[i-1]]h[[i-1]]^2/6)+b[[i-1]](xx-x[[i-1]])+s[[i-1]](x[[i]]-xx)^3/(6h[[i-1]])+s[[i]](xx-x[[i-1]])^3/(6h[[i-1]]);
Do[
Print["C(t) = ",Expand[c[t,i]] ],{i,2,n}]
Spline[t_]:=If[t<x[[2]],c[t,2],If[t<x[[3]],c[t,3],If[t<x[[4]],c[t,4],c[t,5]]]];
sls=Plot[Spline[t],{t,x[[1]],x[[n]]}];
Show[sls,slpod]

(* Definiranje podataka preko funkcije *)
 n=10;SeedRandom[17];
f[t_]:=Sin[t]+1;
slf=Plot[f[t],{t,0,2Pi}];
x=Table[i 2Pi/n,{i,n}];
y=Table[f[x[[i]]]+Random[NormalDistribution[0,.05]],{i,n}]; 
pod=Table[{x[[i]],y[[i]]},{i,n}];
Print["podaci: ",TableForm[pod]]
slpod=ListPlot[pod,PlotStyle->{AbsolutePointSize[5]},
		PlotRange->{0,Max[y]}];
Show[slf,slpod]