Listor / Teknisk fysik / Miniprojekt Lutorp

Miniprojekt Lutorp

Inledning

På ett konsultuppdrag från Gaussby kommun har vi blivit anlitade att skapa ett program som simulerar trycket i det nya vattenledningsverket de planerar att bygga i det nya bostadsområdet Lutorp. Med trycket menas skillnaden mellan ledningarnas knytpunkter och den omgivna atmosfären. I våra beräkningar sätts därför en skala där det atmosfäriska trycket definieras till 0.

Det personliga syftet med projektet är att få en djupare förståelse för de berörda momenten i den kursen Beräkningsvetenskap 1. Moment som kommer behandlas är grundläggande programmeringsstrukturer i Matlab, problemlösningsmetodik samt lösning av linjära ekvationssystem med LU-uppdelning med pivotering.

Problemlösning

Vi hade fått till uppgift av Gaussby kommun att beräkna trycket i vattenledningarna i det planerade bostadsområdet Lutorp. Till att lösa uppgiften valde vi att använda programmet MATLAB, dels eftersom det klarar av komplicerade matematiska beräkningarna och dels eftersom beräkningsprocedurerna kan automatiseras genom skriptning.

Del 1:

Först skrev vi ett enkelt program som beräknade trycket i ett litet vattenledningssystem med fyra knutpunkter. Sambandet mellan trycken i de olika knutpunkterna fås av ett ekvationssystem där vänsterleden representeras av en 4x4-matris där varje rad representerar vänsterledet i en ekvation, och varje kolumn är koeffecienterna framför ett visst tryckvärde. Element 1,1 i matrisen är alltså koeffecienten framför p1 i den första ekvationen. Högerledet och variablerna som representerar trycket i de olika punkterna är var sin kolonnmatris. Ekvationssystem löses sedan genom att använda MATLAB:s backslash-operator, som använder sig av sk Gausseliminering för att lösa ekvationssystemet.

Del 2:

En enkel meny programmerades så att användaren kan välja vilket av de tre alternativen av vattennät den vill utföra beräkningar på. När det valet gjorts uppmanas användaren mata in ett värde på vattenreservoarens tryck. En matrisekvation ställs upp och programmet använder LU-faktorisering för att lösa den:

Vi valde LU-faktorisering eftersom det är det effektivaste sättet vi vet för att få datorn att utföra matrisekvationer. Trycket i vattennätets olika knytpunkter skrivs ut och sparas sedan i en fil som användaren döper. Därefter frågar programmet om en resultatet ska visas i form av en graf. Om användaren inte önskar någon grafisk lösning avslutas programmet.

Del 3:

För att uppfylla kraven för del 3 krävdes bara en del modifikation av programmet i del 2. Vår meny utökades med två extra alternativ i form av de nya varianterna på vattennät med två reservoarer. Om någon av dessa två alternativ valts utförs beräkningar på både de gamla och nya vattennäten så att de sedan kan jämföras med varandra. För att kunna utföra beräkningar de valda matriserna flera gånger utan att starta om programmet satte vi in en while-loop:

Efter att användaren matat in önskat tryck och resultatet redovisats kommer differensen mellan de olika vattennäten visas om något av de två senare alternativen valts i menyn .

Sedan sparas resultatet i en fil som användaren döper och en graf visas vid önskan. Därefter frågar programmet om användaren vill göra beräkningar på ett nytt tryckvärde. Om svaret är ja går programmet tillbaka där while-loopen startar beräkningsproceduren görs om, annars avslutas programmet.

Resultat

Körexempel del 1:

A =

-3.7000 0.5000 0.5000 0.7000

0.5000 -1.1600 0 0.5000

0.5000 0 -1.1600 0.5000

0.7000 0.5000 0.5000 -2.0200

b =

-20

g =

8.1172

5.9893

5.7779

Trycket uttryckt i bar:

8.1172

5.9893

5.7779

Enligt resultatet i del 1 är vattentrycket betydligt lägre i de punkter som ligger långt från reservoaren jämfört med de som ligger nära, vilket är mycket logiskt.

Körexempel del 2:

Uppgift2

Bestäm vattentryck 10

Trycket uttryckt i bar:

8.1172

5.9893

5.7779

Välj filnamnvattennät4

x1 =

1 2 3 4

Vill du se en graf?(ja/nej)

Figur 1: Graf av resultatet från del 2.

Körexempel del 3:

Uppgift3

Bestäm vattentryck 10

Vattennät 1, trycket uttryckt i bar:

7.5342

5.0115

4.0926

3.0834

2.0510

1.6750

1.2619

0.8394

0.6855

0.5165

0.3435

0.2806

0.2114

0.1406

0.1148

0.0865

0.0576

0.0470

0.0354

0.0236

0.0193

0.0146

0.0097

0.0080

0.0061

0.0041

0.0035

0.0029

0.0021

0.0020

Vattennät 2, trycket uttryckt i bar:

7.5348

5.0126

4.0945

3.0864

2.0546

1.6802

1.2697

0.8484

0.6986

0.5356

0.3657

0.3128

0.2583

0.1948

0.1936

0.2013

0.1900

0.2395

0.3159

0.3472

0.4896

0.6999

0.8005

1.1572

1.6807

1.9363

2.8115

4.0946

4.7231

6.8630

Elapsed time is 0.000970 seconds.

Differensen mellan de olika vattennäten uttryckt i bar:

0.0006

0.0011

0.0018

0.0029

0.0035

0.0052

0.0078

0.0090

0.0131

0.0192

0.0222

0.0322

0.0470

0.0542

0.0788

0.1148

0.1324

0.1925

0.2805

0.3236

0.4703

0.6854

0.7908

1.1492

1.6746

1.9321

2.8079

4.0918

4.7210

6.8610

Välj filnamn vatten40_2

Vill du se en graf?(ja/nej) ja

Vill du prova ett annat vattentryck?(ja/nej) nej

Du är nu klar

Figur 2: Graf av resultatet från del 3. Den blå kurvan representerar vattennätet med en reservoar och den orangea kurvan nätet med två reservoarer.

Enligt resultatet i del 3 minskar vatten trycket drastiskt ju längre bort från reservoaren man går. I punkterna längst ifrån är trycket i princip obefintligt. Införandet av en extra vattenreservoar löser inte problemet, det blir visserligen bättre, men punkterna 10-30 ligger för långt bort. Om vattennätet ska innehålla så här många förbindningspunkter är vårt råd att införa ännu en reservoar mellan de nuvarande.

Diskussion

Vi är nöjda med utförandet av vårt projekt och tycker att resultatet är rimligt. LU-faktorisering som lösningsmetod är effektiv och gör att vi sparar en hel del tid när vi upprepar beräkningarna som i del 3. Det finns såklart förbättringsmöjligheter i vårt arbete, antalet variabler skulle kunna minskas och strukturen på programmet skulle kunna göras betydligt bättre. Eftersom det här är vårt första arbete i Matlab och våra kunskaper inom ämnet ännu är låga tror vi det är nåt som kommer med tiden. När vi fått lite mer rutin och lärt oss mer om hur program kan utformas så enkelt som möjligt kommer både resultatet och vår arbetsgång bli betydligt bättre.

Appendix

Del 1:

%Skapar matris A

A=[-3.70,0.50,0.50,0.70;0.50,-1.16,0,0.50;0.50,0,-1.16,0.50;0.70,0.50,0.50,-2.02]

%Skapar vektorn b

b=[-20;0;0;0]

%Löser med Matlabs backslash-kommando

g=A\b

disp('Trycket uttryckt i bar: ');

disp(g)

Del 2:

%Här skapas en meny med de fyra alternativen för vattennät

%n - antalet rader

meny=menu('Välj vattennät','Vattennät 4x4','Vattennät 40x40','Vattennät 400x400','Avsluta');

switch meny

case 1

Matris1='lutorp4.mat';

n=4;

case 2

Matris1='lutorp40.mat';

n=40;

case 3

Matris1='lutorp400.mat';

n=400;

case 4

error('du har inte valt något vattennät')

end

%mat-filerna laddas och programmet LU-faktoriserar

load(Matris1);

p=input('Bestäm vattentryck ');

b=zeros(n,1);

b(1)= (-1)*k*L*p;

[L,U,P]=lu(A);

d=L\(P*b);

g=U\d;

disp('Trycket uttryckt i bar: ');

disp(g);

x=[1:1:n]; %x representerar x-axeln i den graf som ska visas

filnamn=input('Välj filnamn','s');

save(filnamn,'g');

x1=[1:1:n]

%frågar användaren om den vill se en graf,

%om ja--->grafskiss med plot, annars ---> programmet avslutas

svar1=input('Vill du se en graf?(ja/nej)','s');

if strcmp(svar1,'ja');

plot(x1,g);

title('Vattennätets tryck i bar');

ylabel('Bar');

xlabel('Sammanbindningspunkter');

end

Del 3:

%Här skapas en meny med de fem alternativen för vattennät

meny=menu('Välj vattennät','Vattennät 4x4','Vattennät 40x40','Vattennät 400x400','Vattennät 40x40_2','Vattennät 400x400_2','Avsluta');

switch meny

case 1

Matris1='lutorp4.mat';

n=4;

r=0;

case 2

Matris1='lutorp40.mat';

n=40;

r=0;

case 3

Matris1='lutorp400.mat';

n=400;

r=0;

case 4

Matris1='lutorp40.mat';

Matris2='LUtorp40_2.mat';

n=40;

r=1;

case 5

Matris1='lutorp400.mat';

Matris2='LUtorp400_2.mat';

n=400;

r=1;

case 6

error('du har inte valt något vettennät')

end

%mat-filerna laddas och programmet LU-faktoriserar

load(Matris1);

[Q,U,P]=lu(A);

%B laddas endast om alternativ fyra eller fem valts i menyn

if r==1

load(Matris2);

[G,H,J]=lu(A);

end

%Här skapas en while-loop som gör att användaren kan köra så många uträkningar den vill

ekvloop=1;

while ekvloop==1

p=input('Bestäm vattentryck ');

b=zeros(n,1);

b(1) = (-1)*k*L*p;

tic

d=Q\(P*b);

g=U\d;

disp('Trycket uttryckt i bar: ');

disp(g)

%Nu är ekvationen med matris A uträknad

%ekvationen med matris B räknas endast om alternativ fyra elr fem valts

if r==1

w=zeros(n,1);

w(1)=L*(-1*p*k);

w(n)=w(1);

a=G\(J*w);

f=H\a;

disp('Trycket uttryckt i bar: ');

disp(f)

end

toc

diff=f-g;

disp('Differensen mellan de olika vattennäten uttryckt i bar: ')

disp(diff);

x1=[1:n];%x representerar x-axeln i den graf som ska visas

filnamn=input('Välj filnamn','s');

save(filnamn,'g','f','diff');

svar1=input('Vill du se en graf?(ja/nej)','s');

%om svar1=ja, skissa grafen

if strcmp(svar1,'ja')

plot(x1,g);

hold on

plot(x1,f);

hold off

title('Vattennätets tryck i bar');

xlabel('sammanbindningspunkter')

ylabel('Bar')

end

%frågar användaren om den vill köra programmet med nytt vattentryck, om

%input=ja ---> kör loopen igen. annars--->programmet avslutas

svar2=input('Vill du prova ett annat vattentryck?(ja/nej)','s');

if strcmp(svar2,'ja')

ekvloop=1;

else

ekvloop=0;

disp('Du är nu klar')

end

Publiceringsdatum: 2014-11-15