Numerične metode v
fiziki: domače naloge
Vaje iz predmeta Numerične metode v fiziki
Igor Grešovnik
Revision 3, February 2015.
Revision 0: January 2008.
Vsebina:
1 Naloga 1a - aritmetična in
geometrijska zaporedja
2 Naloga 1b - numerična integracija
3 Naloga 1c - pospeševanje avtomobila
4 Dodatna naloga - integracija
5.1 Namig za kopiranje iz konzole
7 Naloga 2b – Numerično odvajanje
8 Naloga 2c – Reševanje nelinearnih enačb
ene spremenljivke
9 Naloga 2d – Matrike in reševanje sistema
linearnih enačb
10.1 Namig - zapisovanje rezultatov v datoteko
10.1 Namig – predstavitev vektorjev in matrik s
tabelami
11 Naloga 3 – Numerično odvajanje, integriranje in
diferencialne enačbe
12 Naloga 4 – reševanje sistemov linearnih enačb
14 Naloga 5b – Minimizacija funkcij
Napiši program za računanje delnih vsot aritmetičnih
in geometrijskih zaporedij:
Program naj izpiše vse delne vsote do danega n, izračunane s seštevanjem členov in z neposredno formulo.
Izpiši člene in
delne vsote od 1 do 20 aritmetičnega zaporedja s podatki
a0 = 5
d = 0,5
in geometrijskega zaporedja s podatki
a = 2
r = 0,5
Implementiraj numerično integracijo realne funkcije ene spremenljivke po
Simpsonovi formuli.
Simpsonovo in trapezno metodo uporabi za izračun integrala
.
in opazuj, kako se povečuje natančnost pri obeh
metodah, ko večamo število intervalov.
Namig:
Enostaven program za numerično integriranje je implementiran v projektu 01integration\
v direktoriju csharp/gresovnik na študentskem FTP računu. Do projekta najlažje
pridete tako, da v Visual Studio C# odprete datoteko test_solution.sln iz tega
direktorija, ki si ga najprej presnamete na lokalni disk.
Navedeni
program uporablja metode in iz projekta 00library (definirane so v datoteki integration.cs).
Metode imajo med svojimi argumenti delegata RealFunction, ki je definiran
v datoteki 0definitions.cs projekta 00library.
Reši nalogo o pospeševanju avtomobila na strain:
http://www2.ijs.si/~zidansek/tekst1.html
Pri reševanju uporabi trapezno in Simpsonovo formulo. Funkcij ni potrebno
risati, lahko jih samo tabeliraš in
skopiraš tabele v poročilo.
To je dodatna naloga za študente, ki niso znali rešiti naloge. Ti študentje
naj najprej rešijo to nalogo, nato pa še redne naloge.
Naloga:
Izračunaj določeni integral v mejah od 0 do 1 funkcije
f(x)
= x ex .
Nedoločeni
integral funkcije je
I(x) = x
ex - ex + C,
Kjer je C nedoločena konstanta.
Integral izračunaj po trapezni
metodi v mejah od s 4, 16, 32,
1. Naredite kopijo direktorija
csharp/template_solution in ga poimenujte s ”<priimek>.<ime>
a. v imenu nadomestite <priimek> in <ime> z vašim priimkom in imenom
2. V tem direktoriju preimenujte datoteko
template_solution.sln v ”<priimek>.<ime>.sln”. To bo rešitev, v
kateri boste naredili nalogo.
3. Odprite rešitev v Visual Studiu in v rešitev
dodajte nov projekt – konzolno aplikacijo z imenom »domaca_naloga_01dodatna_<priimek>«
4. Datoteko v projektu, kjer imate metodo Main,
preimenujte v Program_01.cs.
5. Imenski
prostor (angleško »namespace«) v
tej datoteki preimenujte v Naloga01.
6. V tej datoteki vstavite kak enostaven izpis na
konzolo in poženite program, da vidite, če pravilno deluje.
7. V rešitev vključite obstoječ projekt, ki je v
csharp/gresovnik/00library.
8. V projektu za nalogo naredite novo datoteko v C# z imanom MyIntegration.cs, v njej pa definirajte
razred MyIntegration, ki naj vsebuje
statično metodo MyIntegralTrapezoidal. Razred naj bo definiran v imenskem prostoru MyLib! Metoda naj ima isti podpis (torej argumente istega tipa) kot prva metoda IntegralTrapezoidal definirana v 00integration.cs v projektu
00library, ki ste ga prej vključili. Poskrbeti boste morali, da bodo uporabljeni razredi in drugi tipi vidni v
vaši izvorni datoteki! Metodo boste implementirali pozneje.
9. Naredite novo datoteko v C# z imenom definitions.cs! V njej v imenskem
prostoru MyLib
definirajte razred Definitions, ki naj vsebuje statično metodo Function1Dodatna , katere podpis ustreza delegatu RealFunction
(definirana v projektu 00library v datoteki 0definitions.cs). Metodo
implementirajte tako, da vrne vrednost funkcije, ki jo morate integrirati! . Poskrbeti boste morali, da bodo
uporabljeni razredi in drugi tipi vidni v vaši izvorni datoteki!
10. V razredu Definitions implementirajte še statično
metodo Integral_function1Dodatna. Funkcijo implementirajte tako, da vrne nedoločen integral (brez
nedoločene konstante) funkcije, ki jo morate integrirati. Analitično vrednost
določenega integrala za primerjavo lahko potem izračunate tako, da kličete to
funkcijo pri zgornji meji integral in od nje odštejete vrednost te funkcije pri
spodnji meji.
11. V razredu projekta za rešitev naloge, v
katerem je metoda Main, definirajte metodo void Naloga1Dodatna_IntegrationFromLibrary(void). V tej metodi implementirajte integriranje z uporabo metode za integriranje IntegralTrapezoidal , ki je definirana v integration.cs v
projektu 00library. Integriranje naj bo po navodilih iz naloge! Primerjajte
izračunane integrale z analizično rešitvijo! Na začetku metode vstavite izpis,
ki bo na dobro viden način povedal, katero metodo ste klicali (to bo prišlo
prav pozneje, ko boste klicali več podobnih metod). Tudi tu boste morali
poskrbeti za to, da bodo uporabljeni tipi vidni v datoteki, v kateri
programirate. Metodo kličite iz metode Main in poženite program! Če boste imeli
kakšne napake in bo videti, da rezultati niso pravilni, metodo popravite.
a. Za izračun analitične vrednosti uporabite metodo Integral_function1Dodatna, ki vrne nedoločen integral funkcije, ki jo integrirate! V metodi
definirajte dve spremenljivki tipa RealFunction
(to je delegat, ki je definiran v datoteki 0definitions.cs v datoteki
00library). Eni spremenljivki priredite metodo Function1Dodatna , drudi pa metodo Integral_function1Dodatna iz razreda Definitions v datoteki definitions.cs. Pri naslavljanju teh
metod morate navesti ime razreda (ki ga s piko ločite od imena metode), ker sta
metodi statični. Poskrbeti morate tudi za to, da je razred viden.
12. Implementirajte medoto MyIntegralTrapezoidal v razredu MyIntegration.
13. V razredu projekta za rešitev naloge, v
katerem je metoda Main, definirajte metodo void
Naloga1Dodatna_MyIntegration(void). V tej metodi implementirajte integriranje z uporabo svoje svoje metode MyIntegralTrapezoidal, metodo pa naredite podobno kot tisto za uporabo
integriranja iz knjižničnega projekta 00lbrary. Metodo prav tako kličite iz
metode Main. Klic te metode samo dodajte za klicem prejšnje metode in prejšnji
klic zakomentirajte!
14. V razredu projekta za rešitev naloge, v
katerem je metoda Main, definirajte še metodo void Naloga1Dodatna (void). V tej metodi implementirajte integriranje z
uporabo obeh metod za integracijo (metode iz knjižnice in svoje lastne) tako,
da pri vsakem številu podintervalov tudi primerjate rešitev po eni in drugi
metodi ter analitično rešitev. Tudi klic te metode dodajte v metodo Main in tam zakomentirajte prejšnja
dva klica! V metodi Main boste tako vedno lahko odkomentirali klice metod, ki
jih želite izvajati, in zakomentirali ostale.
15. Ko bodo stvari delovale, preoblikujte
metodo Naloga1Dodatna_MyIntegration tako, da jo boste lahko uporabili za
poljubno funkcijo, pri kateri poznate analitičen izraz za nedoločen integral.
Metodo skopirajte in jo preimenujte v Naloga1Dodatna_MyIntegration_Generic. Metodi (ki je bila prej brez argumentov)
dodajte dva argumenta tipa RealFunction. Prvi argument naj bo metoda, ki vrne
vrednost funkcije, ki jo integrirate. Drugi argument naj bo metoda, ki vrne vrednost nedoločenega integrala te
funkcije. To bo zelo enostavno narediti, ker boste samo obe lokalni delegatski
spremenljivki premaknili iz telesa metode med argumente. Spremenljivki boste
morali seveda prirediti pred klicem metode. Tudi to metodo kličite v metodi
Main, pred klicem pa definitajte obe delegatski spremenljivki, ki ju boste
podali kot argumenta, in jima priredite ustrezne vrednosti.
16. Sedaj nalogo uporabite za integracijo druge
funkcije in sicer f(x) = x ex.
Nedoločen integral te funkcije je I(x)
= x + ex + C. Nalogo rešite enostavno tako, da v razredu Definitions definirate metodi
za izračun vrednosti funkcije in vrednosti njenega nedoločenega integrala, ter
v metodi Main kličete metodo Naloga1Dodatna_MyIntegration_Generic, ki ji preko argumentov prenesete
delegate tipa RealFunction, ki sta inicializirana na ti dve novo definirani
metodi. Metodi
poimenujte Function1Dodatna2 in Integral_function1Dodatna2.
Za vsako vajo napiši kratko in razumljivo poročilo v urejevalniku besedil (npr.
MS Word).
Programe in projekte v C# umesti v dogovorjeno direktorijsko strukturo.
Nekaj tekstov in nalog iz numeričnih metod je na naslednjih straneh:
http://www2.ijs.si/~zidansek/num2003.html
Dodatno gradivo najdeš na FTP strežniku z naslednjimi podatki za dostop:
port: 2222
username: student
password: student
Kopiranje v konzoli je možno le, če to nastaviš pri lastnostih konzole.
Klikneš na kono konzole v levem zgornjem kotu okna, izbereš Properties/Options/
in pod Edit Options vse odkljukaš.
Kljub temu je kopiranje lahko nerodno, če je teksta za več vrstic ali so v
konzoli lomljene vrstice. Problem lahko rešiš tudi tako, da poženeš program, ki
izpisuje v konzolo, in mu preusmeriš standardni izhod (ki je prednastavljen na
konzolo) v izbrano datoteko. To narediš tako, da program poženeš na naslednji
način:
program.exe argumenti... > out.txt
Tu je predpostavljeno ime programa program.exe.
Ko prevedeš projekt v C#, katerega rezultat je program, lahko ta program najdeš
nekje znotraj projektnega direktorija, tipično v direktoriju /bin/debug
relativno glede na projektni direktorij.
Z opisanim načinom so lahko problemi, kadar program bere podatke, ki jih
vstavi uporabnik preko konzole. V takšnem primeru je najbolje v programu
izključiti vsa interaktivna branja in potem pognati program s preusmerjenim
standardnim izhodom.
V rešitvi test_solution v
direktoriju csharp je v projektu 00library datoteka ClassComplex.cs. V tej datoteki je definicija razreda Complex, ki predstavlja kompleksna
števila. V komentarjih so pojasnila o relevantnih konceptih programskega jezika.
Razred Complex vsebuje
statično metodo Example(), ki
demonstirra uporabo definiranega razreda. Metodo lahko kličeš v projektu test_console.
Po zgledu razreda Complex
definiraj razred Cvector, ki
predstavlja vektor kompleksnih števil. V razredu naj bodo definirane
operacije seštevanja in odštevanja vektorja, skalarni produkt dveh vektorjev
(upoštevaj konjugirano simetričnost skalarnega produkta), množenje vektorja s
kompleksnim ali realnim skalarjem, ter indeksna operacija za dostop do
komponent vektorja. Indeksi
komponent naj tečejo od 0 naprej. Vektor naj ima lastnost Dim, ki vrne dimenzijo vektorja (vendar se dimenzija ne da tudi nastaviti preko te lastnosti).
Tip naj vsebuje konstruktor, katerega element je dimenzija števila, pri
konstrukciji pa naj bodo vse komponente inicializirane na (0 + 0 i). Za vajo poleg tega definiraj še
konstruktor, ki mu podaš dimenzijo in eno kompleksno število, na katerega se
postavijo vse komponente vektorja. Komponente naj se kopirajo, vsako posebej
torej alociraš z izrazom new Complex(...), kjer v npr. v konstruktorju podaš
realni in kompleksni del kompleksnega števila, katerega kopijo tvoriš, dobiš ju
z lastnostmi Re
in Im.
Nasveti za izdelavo te naloge:
Še enkrat si poglej razred Complex
v datoteki ClassComplex.cs. Pozorno preberi vse komentarje. Dobro je, da
poskusiš nekajkrat izvesti funkcijo Complex.Example()
(tako, da metodo kličeš iz metode Main
v nekem projektu). Funkcijo lahko spreminjaš in opazuješ, kako razred Complex deluje.
To, da se dimenzije vektorja ne da nastaviti drugače, lahko dosežeš tako,
da definiraš lastnost Dim, v kateri
definiraš blok set kot private. Na ta način je blok set
lastnosti Dim viden le znotraj
razreda in lahko lastnosti Dim le
znotraj razreda prirediš vrednost. V tem primeru bi za podporo lastnosti Dim definiral še privatno spremenljivko
tipa int, v kateri je dejansko shranjena dimenzija (lastnosti – angleško “properties” so sintaktični konstrukti, ki delujejo v tem smislu
bolj kot metode in sami po sebi ničesar ne hranijo). V konstrukrorju v tem
primeru lahko prirediš vrednost lastnost Dim
(zaradi konsistentnosti jo tudi moraš nastaviti), saj je konstruktor del
razreda in mu je zato viden tudi blok get
lastnosti Dim, ki je definiran kot private.
Še boljša možnost je, da v Dim
sploh ne definiraš bloka set, blok get pa naj vrne kar trenutno dimenzijo
tabele, v kateri so shranjene komponente. V tem primeru ni treba posebej
skrbeti za konsistentnost lastnosti Dim,
ker ta vedno vrne dejansko število komponent vektorja. Če tabelo komponent imenuješ
_components, dobiš dimenzijo tabele (število elementov, ki jih vsebuje) z
izrazom _components.Length.
Za primer, kako se definirajo lastnosti, poglej lastnosti Complex.Re in Complex.fi (slednja nima definiranega bloka set in se torej ne da neposredno nastaviti).
Komponente vektorja naj bodo shranjene v tabeli kompleksnih števil, ki ni
vidna navzven (definirana bo kot private ali protected). Dostop do posameznih komponent bo tako možen samo preko
indeksnega operatorja.
Tabelo objektov tipa Complex definiraš
na naslednji način:
private Complex[] _components;
Prostor za tabelo alociraš tako (to je obvezno,
preden začneš v tabelo shranjevati elemente):
_components = new
Complex[dim];
Tabelo komponent alociraš v konstruktorju, ki mu podaš dimenzijo vektorja
(in s tem dimenzijo tabele komponent). S to dimenzijo (argument dim v zgornji vrstici tabele) alociraš
tabelo. Ne pozabi, da indeksi tabele tečejo od 0 do dim-1! V konstruktorju naj se vsa števila postavijo na (0 + 0 i). Lahko bi definiral tudi konstruktor,
ki mu podaš tabelo kompleksnih števil.
V razredu potrebuješ še način za nastavljanje in dostopanje do posameznih
komponent. V ta namen definiraš indeksni
operator. Primer definicije indeksnega operatorja si poglej v razredu Complex. Skica definicije je
public Complex
this[int index]
{ get{ … } set { … } }
V bloku get uporabiš besedo return,
da vrneš zahtevano komponento vektorja, v bloku set pa preko preddefinirane spremenljivke value (ki je istega tipa kot indeksni operator, v tem primeru tipe
Complex) dobiš vrednost, ki jo prirediš določeni komponenti.
Pri
prirejanju komponent (blok set v
definiciji indeksnega operatorja) ne priredi le reference, ampak kompleksno
število kopiraj. To narediš tako, da narediš novo kompleksno število z
operatorjem new, kjer v oklepaju
navedeš komponente kompleksnega števila, ki ga kopiraš. To lahko dosežeš tudi z
uporabo konstruktorja, ki ima za argument kompleksno število, npr.
_components[index] = new Complex(value);
Implementiraj metodi za prve in druge numerične odvode funkcije realne
spremenljivke.
Metodi naj imata naslednje parametre:
·
funkcija,
ki jo odvajamo, v oblki delegata RealFunction,
ki je definiran v datoteki 0definitions.cs
v projektu 00library.
·
vrednost
argumenta funkcije, pri kateri iščemo odvod
·
velikost
koraka pri odvajanju
Delegat RealFunction je definiran
na naslednji način:
public delegate
double RealFunction(double x);
Funkcijo uporabi za odvajanje naslednjih funkcij:
,
kjer je r
v formuli za funkcijo f2(x) enakomerno porazdeljeno naključno
število na intervalu [0,1]. Odvode teh funkcij računaj v 20 ekvidistančnih točkah na intervalu . Za
računanje uporabi korake h = 0,1, , , in . Primerjaj
numerično izračunane vrednosti z
analitičnimi.
Za generacijo enakomerno porazdeljenih naključnih števil na intervalu [0,1]
lahko uporabiš metodo NextDouble() razreda Random:
Random
rnd = new Random();
...
double
x = rnd.NextDouble();
Opozorilo:
Funkcija f2(x) je ekvivalentna funkciji f1(x),
ki ji je dodan beli šum z amplitudo pol milijonine vrednosti funkcije, zaradi
tega sta analitični prvi in drugi odvod funkcije f2(x)
enaka kot pri f1(x). Pri dodatku naključnega člena v
definiciji funkcije f2(x) je treba paziti, da se pri
vsakem izračunu funkcije generira novo naključno število. Programski
generatorji naključnih števil so deterministični in generirajo vedno isto
zaporedje števil po tem, ko jih inicializirano na določen način. Če bi pred
vsakim klicem metode NextDouble() s klicem new Random() na novo
generilali objekt, na katerem metodo kličemo, bi vsakič dobili začetno število
tega zaporedja, torej vedno isto število. Zato je ključno, da objekt razreda Random
naredimo z operatorjem new le enkrat in potem vsakič kličemo NextDouble()
na istem objektu. Enostaven način za dosego tega je, da generator naključnih
števil definiramo kot javno statično spremenljivko v nekem javnem razredu in ga
tam tudi inicializiramo:
...
public
class MyClass
{
...
public
static Random rnd = new Random();
...
}
Za numerično odvajanje uporabi
formuli
,
kjer
je h korak numeričnega odvajanja.
Implementiraj
metodi za dva postopka reševanja nelinearnih enačb - z Newtonovo metodo in nato
še ali s sekantno metodo ali z bisekcijo.
Kriterij
za ustavitev iteracije naj bo vezan na absolutno vrednost funkcije. Enačba naj
bo podana v standardni obliki
.
Pri
Newtonovi metodi naj bodo parametri funkcija in
njen odvod (parametra
naj bosta oba tipa RealFunction),
začetna vrednost neodvisne spremenljivke in toleranca, pod katero mora pasti
absolutna vrednost funkcije f
(toleranca mora biti pozitivno število, v metodah je to potrebno preveriti in
vreči izjemo, če ta pogoj ni izpolnjen).
Pri
sekantni metodi naj bodo parametri funkcija , izhodiščni
vrednosti neodvisne spremenljivke in toleranca, pod katero mora pasti absolutna
vrednost funkcije f. Funkcijski
parameter naj bo tipa RealFunction. Metoda
naj preveri, ali sta vrednosti funkcije v izhodiščnih točkah nasprotno
predznačeni, in naj vrže izjemo, če nista.
Pri
bisekciji naj bodo parametri funkcija , izhodiščni
vrednosti neodvisne spremenljivke in toleranca, pod katero mora pasti absolutna
vrednost funkcije f. Metoda naj tudi
tu preveri, ali sta vrednosti funkcije v izhodiščnih točkah nasprotno
predznačeni, in naj vrže izjemo, če nista. Uporabiš lahko različne pristope za
iskanje dveh začetnih točk, v katerih je funkcija različno predznačena. Lahko na primer
implementiraš metodo, ki izpiše vrednosti funkcije v 20 ali 30 ekvidistančnih
točkah na danem intervalu in s pomočjo te metode na roke (s premikanjem krajišč
intervala) poiščeš dve takšni točki. Do dveh primernih točk lahko prideš tudi z
analizo grafa funkcije, ali pa razviješ poseben algoritem za iskanje takšnih
točk.
Delegatski
tip RealFunction je definiran v csharp\gresovnik\00library\0definitions.cs.
Pri sekantni metodi in bisekciji Priporočljivo je, da pri vsaki od metod uvedeš
še dodatni parameter, ki določa največje dovoljeno število iteracij.
Iteracijski postopek naj se ustavi, če je preseženo največje dovoljeno število
iteracij, in sicer ne glede na to, ali je absolutna vrednost funkcije, katere
ničlo iščemo, že padla pod predpisano toleranco ali ne.
Metodi
v prvem delu naloge uporabi pri reševanju naslednjih enačb:
.
V
vseh primerih pri sekantni metodi uporabi izhodiščni vrednosti in , pri Newtonovi metodi pa začetni približek .
V
drugem delu naloge rešuj z Newtonovo metodo razred enačb
za
10 vrednosti t, ki so enakomerno
razdeljene na intervalu . Začetni približek
naj bo vedno .
V
csharp\gresovnik\00library\ClassMatrixGeneric.cs
in
csharp\gresovnik\00library\ClassVectorGeneric.cs
so
definirani razredi, ki predstavljajo matrike in vektorje.
Najprej
dopolni definicije teh razredov za računanje produktov matrik in produktov
matrik z vektorji. Nato implementiraj še reševajnje sistemov linearnih enačb z
Gaussovo eliminacijsko metodo.
Gaussovo
eliminacijsko metodo uporabi pri reševanju sistema enačb
.
Po
rešitvi sistema enačb izračunaj normo , kjer je izračunana rešitev.
Za vsako napiši kratko in razumljivo poročilo v urejevalniku besedil (npr.
MS Word).
Programe in projekte v C# umesti v dogovorjeno direktorijsko strukturo.
Nekaj tekstov in nalog iz numeričnih metod je na naslednjih straneh:
http://www2.ijs.si/~zidansek/num2003.html
Dodatno gradivo najdeš na FTP strežniku z naslednjimi podatki za dostop:
port: 2222
username: student
password: student
V datoteko se zapisujejo rezultati podobno kot v konzolo, za to lahko
uporabiš razred StreamWriter.
Primer:
static
void FileTest()
{
StreamWriter
sw = new StreamWriter("c:\\temp\\test.txt",
false /* append */);
sw.WriteLine("Zapis v
testno datoteko... ");
sw.WriteLine(" Cas: " + DateTime.Now.ToString());
sw.WriteLine(" Pi ~ " + Math.PI.ToString());
sw.WriteLine("Konec
zapisa. ");
sw.Close();
}
Za predstavitev vektorjev pogosto uporabljamo enodimenzionalne tabele števil, za predstavitev matrik pa dvodimenzionalne tabele oziroma tabele tabel števil (“jagged arrays” v angleščini). V jeziku C# se tabele obnašajo kot objekti in imajo zapisano svojo dolžino v lastnosti Length.
Primer:
Console.WriteLine(Environment.NewLine + "Matrix
- vector multiplication example:" + Environment.NewLine);
double[] x, b; // two
1D arrays to represent vectors
double[][] A; // 2D
array to represent a matrix (note double square brackets [][])
int d1 = 2, d2 = 3; // dimensions of a rectangular matrix
//
Allocate vector b:
x = new double[d2];
//
Allocate the matrix:
A = new double[d1][]; //
allocate outer-most array (array of arrays of double)
for (int i = 0; i < d1; ++i)
A[i] = new double[d2]; //
allocate inner-most arrays (arrays of double, representing rows)
//
Assign elements of vector x:
for (int i = 0; i < x.Length; ++i) // note
use of the Length property
x[i] = (double)i;
//
Assign elements of matrix A:
for (int i = 0; i < A.Length; ++i)
for (int j = 0;
j < A[i].Length; ++j) // note use of the Length property
A[i][j] = (double)i +
0.1 * (double)j; // note
double square brackets [][]
//
Perform multiplication of matrix A and vector x, store result in b:
b = new double[A.Length]; //
first, allocate the vector; use number of rows of A as its dimension
for (int i = 0; i < A.Length; ++i)
{
// Each elemennt of b will become dot product of the
correspondig row of A and of x:
double
product = 0; // auxiliary variable to calculate the dot product
for (int j = 0;
j < A[i].Length; ++j)
product += A[i][j] * x[j];
b[i] = product;
}
// Now,
output elements of the matrix and vectors:
Console.WriteLine("Matrix
A (dimension " + A.Length + "*" + A[0].Length + "):");
for (int i = 0; i < A.Length; ++i)
{
Console.Write("
");
for (int j = 0;
j < A[i].Length; ++j)
Console.Write(A[i][j]
+ " ");
Console.WriteLine();
}
Console.WriteLine();
Console.WriteLine("Vector
x (dimension " + x.Length + "):");
for (int i = 0; i < x.Length; ++i)
Console.WriteLine("
" + x[i] + "
");
Console.WriteLine();
Console.WriteLine("Vector
b = A*x (dimension " +
b.Length + "):");
for (int i = 0; i < b.Length; ++i)
Console.WriteLine("
" + b[i] + "
");
Console.WriteLine(Environment.NewLine + "Matrix
- vector multiplication example finished." + Environment.NewLine);
Ta naloga je sestavljena na podoben način kot naloge, ki jih bodo dobili
študentje, ki v predvidenem roku pred vajami ne bodo uspeli oddati domače
naloge s prejšnjih vaj. Namen takšnih nlog je pridobivanje dodatnih praktišnih
izkušenj, ki bodo v pomoč pri reševanju podobnih nalog, kot so bile podane pri
domačih nalogah. Poleg same naloge je orisan tudi način reševanja.
V tej nalogi izvedemo odvajanje, integracijo in reševanje diferencialne
enačbe na danem intervalu.
Rešujemo diferencialni enačbi z naslednjim začetnim pogojem na danem
intervalu:
.
ter
.
Diferencialni enačbi rešujemo po Eulerjevi metodi. Interval, na katerem
rešujemo enačbi, razdelimo na n enako
dolgih podintervalov dolžine h. Po
vrsti izračunamo vrednosti funkcije v krajiščnih točkah podintervalov po
formuli
,
pri čemer je začetna vrednost funkcije podana. V vsaki iteraciji postopka izračunamo iz
vrednosti, ki smo jo
izračunali v prejšnji iteraciji, in iz diferencialne enačbe.
Numerična rešitev diferencialne enačbe bo predstavljena s tabelo n+1 vrednosti funkcije v vozliščih, ki
delijo interval reševanja diferencialne enačbe na n podintervalov. Pri reševanju diferencialne enačbe poleg vrednosti
funkcije v vozliščih neposredno dobimo tudi približke odvodov. To lahko
uporabimo za preizkus, ali smo diferencialno enačbo pravilno reševali. Preizkus
naredimo tako, da iz izračunanih vrednosti funkcije izračunamo numerične odvode
in jih primerjamo z odvodi, ki jih izrčunamo iz predpisa diferencialne enačbe.
Numerične odvode izračunamo po formuli
.
Odvod v točki postavimo na isto vrednost, kot je izračunan
odvod v .
Iz predpisa diferencialne enačbe izračunamo npr. odvode za 1. funkcijo na
nasleden način:
.
Vidimo, da pri računanju odvodov po predpisu uporabimo samo vrednosti v
točki, v kateri računamo odvod. Za preizkus, ali smo metodo pravilno
implementirali, lahko preverimo ujemanje odvodov izračunanih po obeh metodah.
Poleg navedenega lahko izračunamo še določeni integral odvoda funkcije (ki
je stranski produkt pri reševanju diferencialne enačbe) in ga primerjamo z
rešitvijo diferencialne enačbe. Integral odvoda želimo primerjati z rešitvijo
enačbe v vseh vozliščnih točkah delitve, ne le v zgornji meji integracijskega
intervala. Rabili bomo torej metodo za numerično integriranje, ki izračuna
integral v vseh točkah. Za to lahko uporabimo običajno metodo, ki jo pokličemo
za vsak interval posebej. Tak način pa je zelo neučinkovit,
zato raje sprogramiramo metodo, ki sproti akumulira vrednost integrala v
vozliščih po vrstnem redu. Uporabili bomo trapezno pravilo, pri katerem
izračunamo integral F funnkcije f v točki ,
,
po naslednjem iteracijskem predpisu:
.
Enostavno lahko preverimo, da imata obe diferencialni enačbi, ki ju
rešujemo, isto analitično rešitev
.
Analitični odvod rešitve je
.
To bomo uporabili pri primerjavi napak metod.
V okviru te naloge naredi naslednje:
·
Vsako
posebej reši obe diferencialni enačbi, pri tem tabeliraj izračunane verdnosti in odvode.
·
Tabeliraj
analitično rešitev in njen odvod
·
Izračunaj
numerični odvod analitične funkcije in numerični integral njenega analitičnega
odvoda .
·
Izračunaj
numerični integral tabeliranega odvoda funkcije, ki ga dobiš pri reševanju
diferencialne enačbe, ter numerični odvod tabelirane rešitve diferencialne
enačbe za obe funkciji.
Za vrednosti n=30, n=100, n=200 in n=1000
primerjaj:
·
Analitično
funkcijo z numeričnim integralom
·
Analitično
funkcijo z numeričnim odvodom
·
Reštev
diferencialnih enačb z analitično funkcijo
·
Odvod
funkcije, ki ga dobiš pri reševanju diferencialnih enačb, z analitičnim odvodom
·
Numerični
odvod funkcije, ki jo dobiš pri reševanju diferencialnih enačb, z odvodom, ki
ga dobiš pri reševanju diferencialnih enačb.
·
Numerični
integral odvoda, ki ga dobiš pri reševanju diferencialnih enačb, z numerično
rešitvijo diferencialne enačbe
·
Numerični
odvod funkcije, ki jo dobiš pri reševanju diferencialnih enačb, z analitičnim
odvodom
Pri n=30 primerjaj vrednosti v
vseh točkah. Pri ostalih izračunaj
samo mero za napako, ki je definirana kot
.
Tu je približna vrednost funkcje v točki xi.
Pri reševanju takšnih nalog si najprej pripravimo načrt reševanja. Najprej
identificiarmo, kakšne metode bomo rabili za rešitev problema, od najbolj
osnovnih do tistih, ki so odvisne od le-teh. Nato te osnovne gradnike povežemo
v višjenivojske in rešimo roblem po delih.
Za rešitev pričujočega primera bomo rabili naslednje osnovne metode::
·
metoda
za vzorčenje, ki naredi tabelo n+1 točk
·
metoda
za tabeliranje brednosti funkcije, ki iz zabele naredi tabelo vrednosti . Argumet
metode bo delegat, ki predstavlja realno funkcijo ene spremenljivke.
·
Metoda
za integriranje tabele tako, da so v tabeli rezultatov shranjeni vsi novi
rezultati
·
Metoda za numerično odvajanje na podoben način
·
Metoda
za numerično intgiranje na podoben način.
·
Metoda
za reševanje diferencialne enačbe. Vhodna podatka sta tabela abscis in delegat,
ki ima dva realna argumenta (absciso točke in vrednost funkcije v tej točko)
·
Metoda
za primerjavo dveh tabel
·
Mnetoda
za izračun napake iz dveh tabel
Metode za numerično integriranje, odvajanje in reševanje diferencialnih
enačb sprogramiramo tako, da kot vhodni podatek vzamejo tabele vrednosti in ne
funkcijske predpise.
V tej nalogi rešujemo sisteme linearnih enačb oblike
,
kjer je A matrika sistema z n
vrsticami in n stolpci, b vektor desnih strani dolžine n
in x vektor rešitev, ki zadošča
zgornji enačbi.
Sisteme bomo reševali z uporabo numerične knjižnice Math.Net Iridium. Demonstracija uporabe knjižnice za množenje
matrik in reševanje sistemov enačb je na FTP v direktoriju csharp\gresovnik\05LinearAlgebra\ .
Reši sisteme linearnih enačb z n neznankami, kjer je . Matriko sistema A in vektor desnih strani b generiraj naključno z uporabo metode Matrix.Random(). Vsak sistem reši najprej z uporabo razcepa LU, nato pa isti sistem še z uporabo razcepa QR. Za vsak n primerjaj normo ostanka
,
Število neznank povečuj za faktor 2, dokler skupni
čas rešvevanja sistema enačb po obeh metodah ne postane večji od minute. Pri
tem lahko opazuješ, kako je do neke meje čas računanja komaj opazen, nato pa
zelo hitro narašča s številom neznank. To je zaradi tega, ker je vodilni člen
pri številu operacij za obe metodi sorazmeren z .
Naredi tabelo, s katero primerjaš norme ostankov pri obeh metodi za vse n, pri katerih rešuješ sistem enačb. Pri
katerem n (približno) čas računana
preseže minuto?
Ko bo čas presegel minuto, izvedi za občutek reševanje sistema enačb še za
nekaj naslednjih n, ki ustrezajo
potencam števila 2.
Pri delo s knjižnico Math.Net Iridium
za predstavitev desnih strani in rešitve sistema ne uporabiš tipa Vector, ampak kar tip Matrix dimenzije n*1. Za normo uporabiš
Forbeniusovo normo, ki pri matrikah z enim stolpcem ustreza evklidski normi
ustreznega vektorja.
Z metodo po svoji izbiri (ali s streljanjem ali z relaksacijo) reši robni problem
.
Analitična rešitev je
.
Problem reši pri različnih delitvah intervala: n=50 in n = 200.
Opise metod najdeš med drugim na naslednjih straneh:
·
http://homepage.univie.ac.at/franz.vesely/cp_tut/nol2h/new/c4od_s3bv.html
·
http://www.aip.de/groups/soe/local/numres/bookcpdf/c17-1.pdf
(streljanje)
·
http://www.aip.de/groups/soe/local/numres/bookcpdf/c17-3.pdf
(relaksacija)
Z metodo iz knjižnice najdi minimum funkcije
.
Pri vsakem približku, pri katerem metoda izračuna vrednost funkcije, izpiši
vrednosti neodvisnih spremenljivk, vrednost funkcije in oddaljenost od
minimuma. Minimum funkcije je v (1,1).
Za minimizacijo uporabi metdo L-BFGS iz knjižnice alglib.net.
Spletna stran knjižnice je na naslovu
http://www.alglib.net/ . Na spletni strani
dobiš izvorno kodo metode, ki jo vključiš v svoj projekt. Spletna stran vsebuje
še razlago algoritma in navodila za uporabo metode s primeri. Potrebne datoteke
v C# s spletne strani vključiš v svoj projekt.
Koda metode L-BFGS in njena spletna stran z navodili sta vključeni tudi
rešitev test-solution v projektu AlgLib, ki se nahaja v direktoriju
csharp/external/ALGLIB/. Navodila so v datoteki doc/lbfgs.php.htm, odprete jih
lahko tudi neposredno iz Visual Studia (z desnim gumbom kliknete na datoteko in
izberete ”View in Browser”). Metodo L-BFGS lahko uporabite enostavno tako, da v
svojo rešitev vključite projekt AlgLib in se sklicujete nanj. V projekt je
vključena starejša verzija kode in je tam le informativno ter za primer, da bi
bila spletna stran projekta AlgLib nedostopna.