Algebra Symboliczna
Wykład VIII

Andrzej Odrzywolek
Instytut Fizyki, Zakład Teorii Względnosci i Astrofizyki
17-12-2008, środa, 13:15

Dane kontaktowe

dr Andrzej Odrzywołek
pokój 447, IV pietro
E-mail: odrzywolek@th.if.uj.edu.pl
Wykład: środy 13.15-15.00, sala 128
Cwiczenia: piątki 10.30-12.00
Konsultacje: środy 11-13, czwartki 12-13
WWW: http://ribes.if.uj.edu.pl/alsymb/

UWAGA! Notebook zawiera wcześniejsze obliczenia. W celu jego samodzielnego przeliczenia wybierz z Menu Cell->Delete AllOutput.

Analiza danych

,,Danymi'' określamy zbiór liczb, którymi zamierzamy manipulować w pewien sposób. Ich pochodzenie nie jest istotne, ale typowo są to wyniki ,,pomiarów'' w realnym świecie lub wyniki skomplikowanych obliczeń numerycznych.

Mathematica nie nadaje się do szybkiego przetwarzania gigantycznych strumieni danych, ale pozwala na wykonywanie na nich nieomal dowolnych operacji matematycznych. Jest idealna do złożonej analizy ,,średniej wielkości'' danych i tworzenia algorytmów ich analizy

Możliwości importu/eksportu danych są bardzo rozbudowane, od plików tekstowych i graficznych  aż do HDF5

Przykład: fitowanie liniowe i nieliniowe prostych formuł

Wczytywanie i przekształcanie danych

Pobieramy zbiór z danymi. Może być zlokalizowany na dysku lub w sieci. Oto przykładowe dane:

Import["https://braeburn.pha.jhu.edu/~ariess/R06/Davis07_R07_WV07.dat"]
(* Import["C:/AO/WykladMATHEMATICA 2008/Davis07_R07_WV07.dat"] *)

Jeżeli danych jest bardzo dużo, możemy wyświetlić ich część (np: 20 linii) :

Short[%,20]

Jak widać plik zawiera interesujace nas liczby i nagłówek omawiający dane. Jest to sytuacja typowa; aby dostać same liczby odrzucamy pierwsze 19 linii instrukcją Drop. "Table" oznacza dene w postaci wierszy i kolumn tablice. Inne typowe konwersje to "TSV" (dane oddzielone tabulatorami, np. wygenerowane przez printf("%d\t%lf\n",n, x); ) "CSV" (przecinkiem), "Lines" (linie w postaci tekstowej ) itd.

sn= Drop[Import["https://braeburn.pha.jhu.edu/~ariess/R06/Davis07_R07_WV07.dat","Table"] ,19] 
(* Drop[Import["C:/AO/WykladMATHEMATICA 2008/Davis07_R07_WV07.dat","Table"] ,19]*)

Ilość kolumn i wierszy można odczytać poleceniem Dimensions. Zwraca ono listę wymiarów, pierwszy to liczba wierszy:

num=Dimensions[sn] //First

Istnieje bardzo dużo sposobów operowania takimi tablicami. Chcemy narysować wykres trzeciej kolumny jako funkcji drugiej. Mój ulubiony sposób to przepisanie jawnie elementów do nowej tablicy za pomocą Table:

zm=Table[{sn[[ii,2]],sn[[ii,3]]},{ii,1,num}];

Inny sposób to manipulowanie wierszami i kolumnami. Nie trzeba ,,znać'' wymiarów tablicy:

zmINV={Transpose[sn][[3]],Transpose[sn][[2]]}//Transpose; (* osie zostały zamienione ! *)

zmINV[[1]]

Pierwszą czynnością jest narysowanie wykresu. Pokazuje to przede wszystkim czy dane zostały odczytane poprawnie:

ListPlot[zm]

Wykresy w skali liniowej i logarytmicznej

Aby zorientować się w charakterze zależności funkcyjnej (o ile nie jest znana z góry) możemy narysować 4 rodzaje wykresów: liniowy, logarytmiczny,   podwójnie logarytmiczny, półlogarytmiczny.

funkcje={1/x, x, x^3, Exp[x], Exp[x]+100Exp[-x], x^5 Exp[-x], Log[x]+10}

f=Table[Tooltip[funkcje[[ii]],funkcje[[ii]]],{ii,1,7}]

wLin=Plot[f,{x,0.001,10},PlotRange->{0,10},PlotLabel->"Liniowy"]

wLog=LogPlot[f,{x,0.001,10},PlotRange->{0,10},PlotLabel->"Logarytmiczny"]

wLogLog=LogLogPlot[f,{x,0.001,10},PlotRange->{0,10},PlotLabel->"Podwójnie logartmiczny (LogLog)"]

wLogLinear=LogLinearPlot[f,{x,0.001,10},PlotRange->{0,10},PlotLabel->"Półlogarytmiczny (Semilog, LogLinear)"]

w={wLin,wLog,wLogLog,wLogLinear};

Z poniższego przykładu nalezy zapamiętać, że linia prosta na zwyklym wykresie to funkcja a x+b, na LogPlot funkcje wykładnicze (a Exp[b x])
dają proste, a na wykresie LogLogPlot funkcje potęgowe (a x^b ):

Manipulate[Show[w[[ii]]],{ii,{1,2,3,4}}]

Sprawdzamy pod tym kątem nasze dane:

ListPlot[zm]

ListLogPlot[zm]

ListLogLogPlot[zm]

ListLogLinearPlot[zm]

Prosty fit nieliniowy

Ponieważ na wykresie ListLogLogPlot dane układają się na linii prostej przyjmujemy roboczą hipotezę tymczasową, że m = a z^b. Oczywiście, na pierwszy rzut oka równie dobrze pasuje linia prosta do ostatniego wykresu (m = a Log[z] + b ). Manipulate idealnie nadaje się do wstępnego dopasowania:

Manipulate[Show[ListLogLogPlot[zm],LogLogPlot[a z^b,{z,0.01,10.0}]],{{a,40},10,50},{{b,0.1},0.01,0.2}]

Z Manipulate'ora ustalamy, że a=44.6 i b=0.064. Najlepsze dopasowanie znajduje FindFit, który operuje także na funkcjach nieliniowych w parametrach. Dla przykładu funkcja A Cos[x] + B Log[x] jest liniowa w parametrach A i B (daje kombinację liniową pewnych funkcji) natomiast A Sin[B x]  czy 
a x^b już nie są liniowa w parametrach.

FindFit[zm,a z^b, {{a,44.6},{b,0.064}},z] (* wartości startowe nie są w tym prostym przypadku niezbędne *)

FindFit[zm,a z^b, {a,b},z]

Fitowanie liniowe

Należy pamiętać, że fitowanie funkcji nieliniowej nie daje 100% pewnego rezultatu. Dlatego warto znać alternatywne podejście. Najpierw obliczamy Log danych:

zmLogLog=Log[zm];  (* tak wywołana funkcja, tutaj Log oblicza Log z każdego elementu tablicy (listy) zm! NIE JEST to Log w sensie logarytmu z macierzy ! Funkcje macierzowe mają specjalne nazwy np. MatrixExp itp.*)

A następnie znajdujemy dopasowanie funkcji liniowej, co daje gwarantowany (globalne minimum sumy kwadratów) wynik:

Fit[zmLogLog,{1,X}, X]

Jeżeli interesują nas bardziej szczegółowe informacje (w szczególności błędy znalezionych parametrów), używamy LinearModelFit:

fit=LinearModelFit[zmLogLog,{1,X},X]

Jak widać mamy do dyspozycji szeroki wachlarz informacji na temat procesu fitowania:

info=fit["Properties"]

fit["BestFitParameters"] (* dopasowane parametry *)

fit["ParameterConfidenceIntervals",ConfidenceLevel->0.99] (* przedział ufności na poziomie 99.9% (tzw. 3 σ )*)

fit["ParameterErrors",ConfidenceLevel->0.9](*  błędy *)

fit["ParameterTable"] (* współczynniki i ich błędy *)

Nie szkodzi jednak narysować tej lini prostej, mogło dojśc do pomyłki itp:

Fit[zmLogLog,{1,X}, X]

Show[ListPlot[zmLogLog],Plot[%,{X,-5,1}]]

Zaawansowana analiza nieliniowa

Podobnie można analizować fitowanie nieliniowe; mamy do tego specjalną komende:

fitNieliniowy1=NonlinearModelFit[zm,a z^b, {a,b},z]

fitNieliniowy1["ParameterTable"] (* współczynniki i ich błędy *)

{fitNieliniowy1["AdjustedRSquared"],fit["AdjustedRSquared"]}

A teraz uważny student zauważa, że jak narysujemy funkcję odwrotną (lub wykres LogLinear ) to do danych pasuje inny wcześniej wspomniany wzór: a Log10[z] + b

fitNieliniowy2=NonlinearModelFit[zm, a Log10[z]+b, {a,b}, z]

Na podstawie błędów dopasowania fitu nie jesteśmy w stanie okreslić który ze wzorów: m = a z^b czy m= a Log10[z] + b pasuje lepiej:

{fitNieliniowy1["AdjustedRSquared"],fitNieliniowy2["AdjustedRSquared"]}

fitNieliniowy1["ParameterTable"] (* współczynniki i ich błędy *)

fitNieliniowy2["ParameterTable"] (* współczynniki i ich błędy *)

Być może uwzględnienie błędów pomiaru coś wyjaśni?

Analiza i wykresy z uwzględnieniem błędów ,,pomiarów''

Needs["ErrorBarPlots`"]

zmErrors=Table[{{sn[[ii,2]],sn[[ii,3]]},ErrorBar[sn[[ii,4]]]},{ii,1,num}];

wykresErr=ErrorListPlot[zmErrors]

wagi=Table[1/sn[[ii,4]]^2,{ii,1,num}];

fitNieliniowy1w=NonlinearModelFit[zm,a z^b, {a,b},z,Weights->wagi]

fitNieliniowy2w=NonlinearModelFit[zm,a Log10[z]+b, {a,b},z,Weights->wagi]

fitNieliniowy2["ParameterTable"] (* współczynniki i ich błędy *)

fitNieliniowy2w["ParameterTable"]

{fitNieliniowy1["AdjustedRSquared"],fitNieliniowy2["AdjustedRSquared"],fitNieliniowy1w["AdjustedRSquared"],fitNieliniowy2w["AdjustedRSquared"]}

Dobre i złe dopasowanie; ocena wizualna i numeryczna

Generalnie mozna dopasować wszystko, im więcej parametrów tym łatwiej:

fitZLY=NonlinearModelFit[zm,a + b /z +c z^2 ,{a,b,c},z,Weights->wagi]

{fitZLY["AdjustedRSquared"],fitNieliniowy2w["AdjustedRSquared"]}

fitZLY["ParameterTable"]

,,Optyczna'' ocena pozwala odrzucić bardzo zły fit, ale niewiele daje informacji o funkcjach a z^b i a Log10[z]+b:

Show[wykresErr,Plot[{fitZLY["BestFit"],fitNieliniowy2w["BestFit"],fitNieliniowy1w["BestFit"]},{z,0.01,2}]]

Ładne wykresy

Poniższy przykład może przydać się jeżeli zależy nam na ,,eleganckiej'' prezentacji wyników.

p={a,b} /.fitNieliniowy2w["BestFitParameters"]

pErr=fitNieliniowy2w["ParameterErrors"]

opis="μ = a Log10[z] +b \n"<>"a="<>ToString[NumberForm[p[[1]],{6,3}]]<>"±"<>ToString[NumberForm[pErr[[1]],{6,3}]]<>"\n"<>"b="<>ToString[NumberForm[p[[2]],{6,3}]]<>"±"<>ToString[NumberForm[pErr[[2]],{6,3}]]

Show[wykresErr,Plot[{fitNieliniowy2w["BestFit"]},{z,0.01,2},PlotStyle->{Black,Dashed}],Frame->True,FrameLabel->{"z","μ [mag]"},Epilog->Inset[opis,{1.5,36}],Axes->None,ImageSize->300]

Export["C:/Fit.pdf",%]

Tego typu analiza jak powyżej niczego nie zakłada na temat zależności funkcyjnej m(z). Tak naprawdę nawet nie wiemy co to jest m i z. Zwykle jednak dane analizujemy w ramach pewnej teorii lub modelu. To wymaga sporej wiedzy. W opisywanym przykladzie dane to zależność jasności obserwowanej m supernowej typu Ia od jej przesunięcia ku czerwieni z. W najprostszym modelu zależnośc ta opisywana jest prawem Hubble'a :

dL==c z/H

W powyższym wzorze dL to odległość jasnościowa, c=299792458 m/s, H - stała Hubble's, z - przesumięcie ku czerwieni. Odległość jasnościową można wyrazić poprzez m jako:

m==5 Log10[dL]+25

%/. dL->c z/H

Solve[%,m]

c=299792.458

FindFit[zm,((5 (5 Log[10]+Log[((c z)/H)]))/(Log[10])),{H},z]

NonlinearModelFit[zm, ((5 (5 Log[10]+Log[((c z)/H)]))/(Log[10])),{H},z]["ParameterTable"]

Manipulator poniżej pokazuje dopasowanie prawa Hubble'a z różnymi wartościami stałej Hubble'a H. Wartość współcześnie przyjmowana do 72 km/s/Mpc z dosyć dużym błędem rzędu kilka km/s/Mpc, ale na pewno nie 50.8±0.7 jak wynika z fitu powyżej. Prawo Hubble's dL=cz/H stosuje się tylko zla z<<1, czyli v<<c.

Manipulate[Show[wykresErr,Plot[{((5 (5 Log[10]+Log[5895.947075851495` z]))/(Log[10])),((5 (5 Log[10]+Log[((c z)/H)]))/(Log[10]))},{z,0,2},PlotStyle->{Black,Dashed}]],{{H,72},30,100}]

Musimy więc dopasować liniowe prawo Hubble'a dla podzbioru danych spełniającego z<<1. Spróbujemy narysować  wykres Hbestfit(zMAX). Musimy w tym celu wybrać dane dla których z<zMAX. Dla wygody tworze tabele zawierającą odległość dL zamiast ,,wielkości gwiazdowej'' m i sortuję po z:

Solve[m==5 Log10[dL]+25,dL]

zdL=Table[{zm[[ii,1]],(10^(-5+(m/5)))/.m->zm[[ii,2]]},{ii,1,num}]//Sort;

z0=zdL[[1]] //First(* najniższe z *)

z1=zdL //Last //First(* najwyższe z *)

Select[zdL,#[[2]]>10000&]

FindFit[Select[zdL,#[[1]]<0.05&],299792.458 z/H,{H},z]

Manipulate[Show[ListPlot[Select[zdL,#[[1]]<zMAX&]],Plot[Evaluate[(299792.458 z/H) /. FindFit[Select[zdL,#[[1]]<zMAX&],299792.458 z/H,{H},z]],{z,0,2},PlotStyle->{Red,Dashed}],Plot[299792.458 z/Hubble,{z,0,2},PlotStyle->Black],PlotLabel->"H (najlepszy fit, czerwony)="<>ToString[H/. FindFit[Select[zdL,#[[1]]<zMAX&],299792.458 z/H,{H},z]],Frame->True, FrameLabel->{"z","dL [Mpc]"}],{{zMAX,1},z0,z1,0.02},{{Hubble,50.8},40,90}]

ListPlot[Table[{zMAX, H /.FindFit[Select[zdL,#[[1]]<zMAX&],c z/H,{H},z]},{zMAX,z0,z1,0.02}] ,Joined->True,InterpolationOrder->0]

Fitowanie skomplikowanych modeli danych w postaci całek, równań różniczkowych itp.

Powyżej widać, że dopiero dopasowanie stałej Hubble'a dla z<<1 pokazuje, że dla z>1 mamy istotne odstępstwo od prawa liniowego. Nie znając zakresu stosowalności tego prawa, doszlibyśmy do wniosku że stała Hubble'a wynosi 40-50 km/s/Mpc. W rzeczywistości prawo to jest nieliniowe, ale nie wynika to bynajmniej z samych danych które analizujemy. Potrzebujemy dodatkowej informacji wynikającej z teorii, w tym wypadku Ogólnej Teorii Względności. Oto wzór na odległość jasnościową we Wszechświecie z zawartością masy Ωm, stałą kosmologiczną Ωλ i stałą Hubble'a H:

(* Integrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] *)

dL1[z_,Ωm_?NumberQ,Ωλ_?NumberQ,H0_?NumberQ]:=Piecewise[{{299792.458/H0(1+z)/Sqrt[1-Ωm-Ωλ]Sinh[Sqrt[1-Ωm-Ωλ] NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ],Ωm+Ωλ<1},{299792.458/H0(1+z)NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ,Ωm+Ωλ==1},{299792.458/H0(1+z)/Sqrt[Ωm+Ωλ-1]Sin[Sqrt[Ωm+Ωλ-1] NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ],Ωm+Ωλ>1}}]

dL1[0.9,0.3, 0.9, 70]

Tutaj jest podana inna forma funkcji powyżej, która zgodnie z manualem pozwala na jej użycie w funkcji FindFit. Teraz funkcja jako argumenty posiada wyłącznie parametry fitowane, a zwraca nie liczbę, ale funkcję (interpolującą). Ten sposób działa również dla innych ,,trudnych'' funkcji np. generowanych przez NDSolve, FindRoot itp.

dL2[Ωm_?NumberQ,Ωλ_?NumberQ,H0_?NumberQ]:=Interpolation[Table[{z,Piecewise[{{299792.458/H0(1+z)/Sqrt[1-Ωm-Ωλ]Sinh[Sqrt[1-Ωm-Ωλ] NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ],Ωm+Ωλ<1},{299792.458/H0(1+z)NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ,Ωm+Ωλ==1},{299792.458/H0(1+z)/Sqrt[Ωm+Ωλ-1]Sin[Sqrt[Ωm+Ωλ-1] NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ],Ωm+Ωλ>1}}]},{z,0,2,0.25}]]

dL2[0.3,0.9,70][0.99]

Definicja poniżej wykorzystuje mechanizm zapamiętywania raz obliczonych wartości funkcji. Przeczytaj w manualu stronę ,,FunctionsThatRememberValuesTheyHaveFound''. Jeżeli ,,normalna'' funkcja jest zdefiniowana np. jako f[t_]:=NIntegrate[x^x Exp[-x^2],{x,-Infinity,t}] to funkcja ,,pamiętająca'' wygląda tak: f[t_]:=f[t]=NIntegrate[x^x Exp[-x^2],{x,-Infinity,t}] . Wartości już obliczone są zapisywane w pamięci a potem tylko odczytywane, co dla funkcji powyżej prowadzi przy fitowaniu do przyspieszenia około 100x !!!

dL3[Ωm_?NumberQ,Ωλ_?NumberQ,H0_?NumberQ]:=(dL3[Ωm,Ωλ,H0]=Interpolation[Table[{z,Piecewise[{{299792.458/H0(1+z)/Sqrt[1-Ωm-Ωλ]Sinh[Sqrt[1-Ωm-Ωλ] NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ],Ωm+Ωλ<1},{299792.458/H0(1+z)NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ,Ωm+Ωλ==1},{299792.458/H0(1+z)/Sqrt[Ωm+Ωλ-1]Sin[Sqrt[Ωm+Ωλ-1] NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ],Ωm+Ωλ>1}}]},{z,0,3,0.25}]])

Alternatywnie, możemy wygenerować tablicę wartości funkcji dL i stworzyć funkcję interpolującą. Trwa to dosyć długo (4 wymiary), ale otrzymaną funkcję można zapisać na dysku w celu późniejszego natychmiastowego użycia. UWAGA! Linie są wykomentowane, bo obliczenia zajmują bardzo dużo czasu i pamięci.

DistributeDefinitions[dL3]

(* dL4=Flatten[ParallelTable[{z,OmegaM, OmegaΛ, H, dL3[OmegaM,OmegaΛ,H][z]//Re},{z,0.001,2.5,0.1},{OmegaM,0.01,2.0,0.1},{OmegaΛ,-1,2,0.1},{H,40,90,5}],3]//Interpolation *)

(* Export["dL.m",dL4] *)

dL4=Import["https://ribes.if.uj.edu.pl/alsymb/Wyklad/dL.m"]

dL1[1.1,0.3,0.8,55]

dL2[0.3,0.8,55][1.1]

dL3[0.3,0.8,55][1.1]

dL4[1.1, 0.3,0.8,55]

Obniżamy żądaną dokładność całkowania do 1%, podobnie dla fitowania. Bez tego nie doczekałem się wyniku po około 15 minutach. UWAGA: używając ,,funkcji pamiętajacej'' (dL3) dostajemy wynik bez ograniczania dokładności po około 2 sekundach (!). Jeżeli ktoś posiada szybki komputer, może pominąć 2 linie poniżej.

SetOptions[NIntegrate,AccuracyGoal->2]

SetOptions[FindFit,AccuracyGoal->2]

dL1[0.3,0.7,72]

dL2[0.3,0.7,72][1.1]

Na początek spróbujemy manipulatora:

Manipulate[Show[ListPlot[zdL],Plot[dL3[OmegaM,OmegaΛ,Hubble][z],{z,0,2}]],{Hubble,40,90},{OmegaM,0,2},{OmegaΛ,-2,2}]

UWAGA! poniższe obliczenia mogą być bardzo czaso- lub pamięciochłonne, w zależności od funkcji (dL1, dL2, dL3 lub dL4). Na ,,słabszych'' PC, najlepiej użyć  dL3 w połączeniu z ograniczeniem żądanej dokładności. Linijki definiujące pozostałe funkcje (powyżej) najlepiej pominąć w celu zaoszczędzenia pamięci.

(* FindFit[zdL,dL[OmegaM,OmegaΛ,Hubble][z],{OmegaM,OmegaΛ,Hubble},z] //Timing (* na słabszym komputerze to może trwać bardzo długo *)

(* {176.891`,{OmegaM->0.348213546549425`,OmegaΛ->1.4964962156776132`,Hubble->75.81960138266479`}} naprawdę to wyliczyło ! *)

NonlinearModelFit[zdL,dL3[OmegaM,OmegaΛ,Hubble][z],{{OmegaM,0.3},{OmegaΛ,0.8},{Hubble,50}},z] //Timing

kosmo=%[[2]]

kosmo["ParameterTable"]

A teraz zakładamy, że znana jest stała Hubble'a H=72 km/s/Mpc:

NonlinearModelFit[zdL,dL3[OmegaM,OmegaΛ,72][z],{{OmegaM,0.3},{OmegaΛ,0.7}},z] //Timing

kosmoH=%[[2]]

kosmoH["ParameterTable"]

A teraz zakładamy, że przestrzeń jest płaska, czyli Ωm+Ωλ==1 (dL4 zżera >1GB pamięci, lepsza jest funkcja dL3 ):

NonlinearModelFit[zdL,{dL3[OmegaM,1-OmegaM,H][z]},{{OmegaM,0.3},{H,72}},z] //Timing

kosmoPlaski=%[[2]]

kosmoPlaski["ParameterTable"]

Funkcje pamiętające

Istnieje w MATHEMATICE mechanizm definiowania funkcji, które ,,zapamiętują'' raz obliczone rezultaty.

Rozważmy 3 definicje tej samej funkcji:

f1[X_]=NIntegrate[Abs[x]^x Exp[-x^2],{x,0,X}] (* przypisanie natychmiastowe = *)

f2[X_]:=NIntegrate[Abs[x]^x Exp[-x^2],{x,0,X}] (* klasyczna definicja funkcji *)

f3[X_]:=f3[X]=NIntegrate[Abs[x]^x Exp[-x^2],{x,0,X}] (* funkcja ,,pamiętająca'' *)

{Plot[f1[t],{t,-3,3}],Plot[f1[t],{t,-3,3}], Plot[f1[t],{t,-3,3}]  } //Timing

{Plot[f2[t],{t,-3,3}],Plot[f2[t],{t,-3,3}], Plot[f2[t],{t,-3,3}]  } //Timing

{Plot[f3[t],{t,-3,3}],Plot[f3[t],{t,-3,3}], Plot[f3[t],{t,-3,3}]  } //Timing

Jak widać gołym okiem, ponowne wywołanie funkcji f3 daje wynik nieomal natychmiastowo:

Plot[f3[t],{t,-4,4}] //Timing

Narysowanie jednak w innym przedziale, nawet zawierającym wcześniejszy, nie daje przyspieszenia:

Plot[f3[t], {t, -1, 1}] // Timing

Analiza przebiegu ,,czasowego'': chłodzenie PC

Nieco egzotyczny jak na ten wykład przykład, ale pokazujący uniwersalność MATHEMATICI w analizie danych. Kod poniżej wczytuje  i rysuje wykresy zależności temperatury (oraz prędkości wentylatorów) od czasu zapisane programami SpeedFan i RivaTuner.

Najpierw ustalamy ,,dzisiejszą'' datę (w przykładzie poniżej jest ona ustalona ,,ręcznie'' na 16 grudnia 2008):

Quit[]

(* now=DateList[] *)

now={2008,12,16,20,49,59.9375`9.530273452750455}

today={now[[1]],now[[2]],now[[3]],0,0,0}

(* abstime=AbsoluteTime[today] *)

abstime=3438374400

Tworzymy końcówkę pliku zawierajacą datę:

stamp=If[now[[3]]>9,ToString[today[[1]]]<>""<>ToString[today[[2]]]<>""<>
ToString[today[[3]]],ToString[today[[1]]]<>""<>ToString[today[[2]]]<>"0"<>
ToString[today[[3]]]]

Nagłówek pliku zawiera informacje o zapisanych danych: czas, temperatury procesorów, chipsetu, dysku, karty graficznej oraz prędkości 3 wentylatorów.

Take[Import["https://ribes.if.uj.edu.pl/alsymb/Wyklad/SFLog"<>stamp<>".csv","Table"],2]

Odczytujemy dane odrzucając nagłówek:

antec=Drop[Import["https://ribes.if.uj.edu.pl/alsymb/Wyklad/SFLog"<>stamp<>".csv","Table"],1];

Rozmiar danych:

num=Dimensions[antec][[1]]

SpeedFan zapisuje liczbę sekund które upłyneły od północy każdego dnia, dlatego czas pomiaru to abstime+czas, gdzie abstime to liczba sekund od 1 stycznia 1900.

c2d=Table[{abstime+antec[[ii,1]],antec[[ii,2]]},{ii,1,num}]; (* temperatura procesora *)
c2d1=Table[{abstime+antec[[ii,1]],antec[[ii,3]]},{ii,1,num}]; (* temperatura rdzenia 1 *)
c2d2=Table[{abstime+antec[[ii,1]],antec[[ii,4]]},{ii,1,num}]; (* temperature rdzenia 2 *)

Dane ,,czasowe'' możemy rysować poleceniem DateListPlot. Jest to zależność temperatur w ciągu całego dnia:

figCPU=DateListPlot[{c2d,c2d1,c2d2},Joined->True,PlotStyle->{Black,Red,Blue},AspectRatio->1/2]

Podobnie możemy narysować obroty wentylatora na procesorze:

c2dFAN=Table[{abstime+antec[[ii,1]],antec[[ii,8]]},{ii,1,num}];

figFAN=DateListPlot[c2dFAN,Joined->True,AspectRatio->1/2,PlotRange->All]

Wykres sprawdzający jak działa automatyczna regulacja obrotów w zależności od temperatury:

GraphicsArray[{{figFAN},{figCPU}}]

Poniżej nieco bardziej złozony przykład, który odczytuje log zapisany przez RivaTuner

p35=Table[{abstime+antec[[ii,1]],antec[[ii,5]]},{ii,1,num}];

g80=Table[{abstime+antec[[ii,1]],antec[[ii,7]]},{ii,1,num}];

samsung=Table[{abstime+antec[[ii,1]],antec[[ii,6]]},{ii,1,num}];

sysFAN=Table[{abstime+antec[[ii,1]],antec[[ii,10]]},{ii,1,num}];

t1=abstime+antec[[1,1]]
t2=abstime+antec[[-1,1]] (* czas który obejmuje log *)

Take[Import["https://ribes.if.uj.edu.pl/alsymb/Wyklad/HardwareMonitoring.hml","CSV"],7]

rt=Drop[Import["https://ribes.if.uj.edu.pl/alsymb/Wyklad/HardwareMonitoring.hml","CSV"],8];

num2=Dimensions[rt]//First

δii=Round[num2/1000] (* nie potrzeba więcej niż 1000 punktów, bo monitor nie ma takiej rozdzielczości ! *)

DateList[{rt[[-1,2]],{"Day","-","Month","-","Year"," ","Hour",":","Minute",":","Second"}}] (* format zapisu czasu i daty jest tu dosyć skomplikowany typu "15-12-2008 15:44:48" *)

DateString[%]

gpuFAN=Table[{DateList[{rt[[ii,2]],{"Day","-","Month","-","Year"," ","Hour",":","Minute",":","Second"}}],rt[[ii,9]]},{ii,1,num2,δii}];

gpuAmbient=Table[{DateList[{rt[[ii,2]],{"Day","-","Month","-","Year"," ","Hour",":","Minute",":","Second"}}],rt[[ii,8]]},{ii,1,num2,δii}];

g80=Table[{DateList[{rt[[ii,2]],{"Day","-","Month","-","Year"," ","Hour",":","Minute",":","Second"}}],rt[[ii,6]]},{ii,1,num2,δii}];

Oto efekt końcowy:

DateListPlot[{Tooltip[sysFAN,"SYS FAN"],Tooltip[c2dFAN,"CPU FAN"],Tooltip[gpuFAN,"GPU Fan"]},Joined->True,PlotStyle->{Red,Blue,Black},AspectRatio->1/3,PlotRange->{{t1,t2},{0,3000}}];

DateListPlot[{Tooltip[g80,"GeForce"],Tooltip[gpuAmbient,"GPU ambient"],Tooltip[p35,"Chipset"],Tooltip[c2d,"CPU"],Tooltip[c2d1,"Core 1"],Tooltip[samsung,"Samsung"]},Joined->True,PlotRange->{{t1,t2},{30,90}},PlotStyle->{Blue,Blue,Black,Red,Red,Green,Green},AspectRatio->1/3];

GraphicsArray[{{%},{%%}}]

Przykład analizy ,,orbity'' eliptycznej (Keplerowskiej)

Naszym zadaniem jest znalezienie parametrów elipsy na podstawie danych punktów opisujących jej fragment. Jest to ,,obcięta'' wersja zagadnienia polegajacego na wyznaczeniu orbity poprzez znane pozycje w kilkunastu punktach. Kepler jako pierwszy rozwiązał to zadanie poprawnie podając, że punkty te leżą na krzywej opisanej wzorem r=p/(1+e Cos[ϕ]).

Generujemy ,,sztuczne dane'':

Quit[]

X0=RandomReal[{-5,5}] (* ognisko elipsy, X *)

Y0=RandomReal[{-5,5}] (* ognisko elipsy Y *)

mimosrod=RandomReal[] (* mimośród elipsy *)

Δϕ=RandomReal[{0,2 Pi}] (* kąt obrotu elipsy *)

ϕMAX=RandomReal[{0,2 Pi}] (* obszar w którym ,,znamy'' pozycje *)

Niestety, MATHEMATICA bezpośrednio nie fituje funkcji o wartościach wektorowych. Jednym ze sposobów jest wyliczenie ,,ręcznie'' sumy kwadratów różnic i minimalizowanie tak zdefiniowanej funkcji poprzez np. NMinimize. Niżej podaję bardziej eleganckie rozwiązanie,działajace w 2 wymiarach, czyli użycie liczb zespolonych. Elipsę traktuję jako funkcję zespoloną parametru rzeczywistego:

dane1=Table[{ϕ,Cos[ϕ]/(1+mimosrod Cos[ϕ+1.2])+RandomReal[]/10+X0+I (Sin[ϕ]/(1+mimosrod Cos[ϕ+1.2])+RandomReal[]/10+Y0)},{ϕ,0,ϕMAX,0.1}]

{Re[Transpose[dane1][[2]]],Im[Transpose[dane1][[2]]]}//Transpose

punkty1=ListPlot[%]

Zawsze polecam Manipulate jako wstęp do fitowania:

Manipulate[Show[punkty1,ParametricPlot[{Cos[ϕ] p/(1+e Cos[ϕ+ϕ0])+x0,Sin[ϕ] p/(1+e Cos[ϕ+ϕ0])+y0},{ϕ,0,2Pi}],PlotRange->{{-10,10},{-10,10}},AspectRatio->1],{p,1,2},{e,0,1},{ϕ0,0, 2 Pi},{x0,-5,5},{y0,-5,5}]

Aby było możliwe fitowanie funkcji o wartościach zespolonych musimy przedefiniować ,,normalizację'':

UWAGA! Wartości parametrów startowych poniżej należy ustalić na podstawie przybliżonego dopasowania z Manipulate powyżej! Bez tego FindFit prawdopodobnie nie da dobrego wyniku, albo wyświetli komunikat o błędzie!

fit1=FindFit[dane1,Cos[ ϕ] p/(1+e Cos[ϕ+ϕ0])+I (Sin[ ϕ] p/(1+e Cos[ϕ+ϕ0]))+x0+I y0,{{p,0.5,1.5},{e,0.1,0.99},{ϕ0,0.1,2.0},{x0,1,3},{y0,-2,2}},ϕ,NormFunction->(Norm[Abs[#1^2]]&)]

Sprawdzamy dopasowanie:

{Cos[ϕ] p/(1+e Cos[ϕ+ϕ0])+x0,Sin[ϕ] p/(1+e Cos[ϕ+ϕ0])+y0} /. fit1

dopasowanie=ParametricPlot[%,{ϕ,0,2 Pi}]

Show[punkty1,dopasowanie,PlotRange->All]

Fitowanie rozwiązania numerycznego równania różniczkowego

Jest to dobrze wytłumaczone w manualu do FindFit, sekcja Applications/Differential Equations  (2). Jest to bardzo podobne do przykładu z fitowaniem modelu kosmologicznego powyżej.

Zawsze można używać funkcji interpolujących, aczkolwiek generowanie wartości może trwać długo, a sama funkcja zajmuje sporo miejsca w pamięci. Musimy też znać z góry spodziewany zakres parametrów, np. z Manipulate. Funkcję jednak możemy zapisać na dysku do późniejszego użytku.

Optymalnym rozwiązaniem jest użycie funkcji pamiętającej f[x_]:=f[x]=.....definicja funkcji ....