1. Starta matlab, starta programmet/gui'n PP med
>>PP. Experimentera! Sök svar på frågorna.
2. Kopiera koden i MyPoissonSolver.m. Studera koden
omsorgsfullt för att förstå hur den är uppbyggd. Komplettera i
funktionen ElementContributions som skall returnera bidragen till
styvhetsmatris och lastvektor från element el. Testa koden med data
motsvarande några problem med känd lösning för att försäkra dig om att
allt fungerar som det ska. Aktuella data specificeras i funktionerna
EqData och BdryData sist i filen MyPoissonSolver.m, eller i separata
filer (med andra namn typ MyEqData och MyBdryData, och motsvarande
anrop MyPoissonSolver(p,t,e,MyEqData,MyBdryData, om inte funktionerna
EqData och BdryData stryks ur MyPoissonSolver.m filen). Glöm sen inte
att försöka hitta några käcka tillämpningar!
3. Utgå gärna från koden i MyHeatEqSolver.m, som
väsentligen bygger på koden MyPoissonSolver för beräkning av styvhets-
och mass matriser och last vektor, vilka nu även kan vara
tidsberoende, och som endast kräver att man specificerar själva
tidsstegningsmetodiken. Än en gång gäller det att få koden
funktionsduglig, testa den, och att tillämpa på problem av
intresse.
Alternativ kan du välja den mera primitiva metoden att hårdkoda in
matriser och lastvektorer utgående från kodskalet under länken "tips" nedan.
Programmet skall förstås vara lätt att modifiera vad gäller t.ex.
steglängderna h och k i rums respektive tidsled och
olika typer av randvillkor. Skall även vara möjligt att studera
tidsutvecklingen av lösningen genom
att lösningen plottas efter varje nytt tidssteg.
Förslag till programstruktur:
Definiera tidssteg och rumssteg/antal element, samt sluttid.
k=.01; %tidssteg
n=20; %antal delintervall/element i rumsled
h=1/n; %rumssteg
T=5; %sluttid
Definerar styvhetsmatris och massmatris.
clear s,m,g %radera tidigare styv och massmatris och vektor g
for i=1:n-1
s(i,i)= %diagonalelement
s(i,i+1)= %superdiagonal
s(i+1,i)= %subdiagonal
m(i,i)= %diagonalelement
m(i,i+1)= %superdiagonal
m(i+1,i)= %subdiagonal
end
s(n,n)= %avvikande diagonalelement i detta fall
m(n,n)= %avvikande diagonalelement i detta fall
Definierar begynnelsevärden för temperaturen u.
x=(1:n)/n; %vektor med x koordinater
u=x;
Initierar grafik och plottar u0:
figure(1)
clf
set(gcf,'renderer','zbuffer') %väljer ritverktyg
handle=line([0 x],[0 u]); %plottar initialdata och ger kurvan en
identitet
set(gca,'YLim',[-1 1])
Startar tidsstegning.
for i=1:T/k %krävs T/k tidssteg
unew= %ger tempen efter i tidssteg uttryckt i tempen u efter i-1 tidssteg
u=unew'; %uppdaterar u (unew på kolonnform)
set(handle,'YData',[0 u]) %uppdaterar plotkurvans "y-data"
drawnow
end
I uppgift b) definieras begynnelsetempen till 0 och tillkommer
randdata u'(1,t)=g(t). Kodas t.ex. genom att definiera funktionen g(t)
som en text i programmets början:
g_av_t='sin(10*t)';
som sedan evalueras i resp tidssteg efter att aktuell tid t
definierats:
t=(i-1)*k;
g(n)=eval(g_av_t)*k;
vilket ger en n-vektor med sista komponent approximativt lika med
integralen av g(t) över det aktuella tidssteget. Detta "kod-skal" kan
kopieras från ~pdef2/PDEF2_98_lab_tips.
7. Modifiera programmet för värmeledningsekvationen.
Last modified: Wed May 5 10:39:29 MET DST 1999