Numeriska metoder f¨
or differentialekvationer
Grundl¨
aggande begrepp
Den enklaste formen f¨or en differentialekvation av f¨orsta ordningen ¨ar
y 0 = h(x).
(1)
En s˚
adan ekvation kan l¨osas direkt. Om
H(x) =
Z
h(x)dx [H 0 (x) = h(x)]
¨ar en primitiv till h(x) s˚
a ¨ar ju
y(x) = H(x) + C
den alm¨anna l¨osningen till (1). Konstanten C best¨ams av n˚
agot begynnelsevillkor.
En linj¨ar differentialekvation av f¨orsta ordningen ¨ar
L(y) ≡ y 0 + g(x)y = h(x).
(2)
H¨ar ¨ar g och h givna kontinuerliga funktioner i ett ¨oppet intervall p˚
a reela axeln x. L kallas
en linj¨
ar differentialoperator (av f¨orsta ordningen) eftersom
L(y1 + y2 ) = (y1 + y2 )0 + g(x)(y1 + y2 ) = y10 + y20 + g(x)y1 + g(x)y2 = L(y1 ) + L(y2 ).
L(αy) = (αy)0 + g(x)(αy) = αy 0 + αg(x)y = αL(y)
och
L(αy1 + βy2 ) = αL(y1 ) + βL(y2 ).
(3)
Om y1 och y2 l¨oser de tv˚
a ekvationerna
L(y) = h1 (x)
respektive
L(y) = h2 (x)
s˚
a l¨oser y1 + y2 ekvationen
L(y) = h1 (x) + h2 (x)
och αy1 l¨oser ekvationen
L(y) = αh1 (x).
Detta kallas superpositionprincipen.
En differentialekvation av f¨orsta ordningen s¨ages vara separabel eller ha separabla variabler
om den kan skrivas p˚
a formen
g(y)y 0 = h(x).
1
(4)
En s˚
adan ekvation kan l¨osas direkt. Betrakta t ex begynnelsev¨
ardesproblemet
(
√
= −k y,
y(0) = h,
dy
dt
y = y(t), t > 0,
Differentialekvationen kan skrivas som
1 dy
= −k
√
y dt
(y > 0).
(5)
Den har separabla variabler. Vi har
d √
1
(2 y) = √ :
dy
y
1 dy
d q
(2 y(t)) = √
,
dt
y dt
(6)
vilket betyder att differentialekvationen kan skrivas som
d √
(2 y) = −k
dt
Detta ¨ar ekvivalent med
√
2 y=
(7)
Z
(−k)dt = −kt + C
√
Begynnelsev¨ardet y(0) = h ger att C = 2 h och allts˚
a
√
y=
√
−k
t+ h
2
(8)
I detta fall kan vi enkla l¨osa y genom kvadrering:
y=
!2
√
−k
t+ h
2
.
(9)
En linj¨ar differentialekvation av andra ordningen ¨ar
M (y) ≡ y 00 + a(x)y 0 + b(x)y = h(x).
(10)
H¨ar ¨ar a, b och h givna kontinuerliga funktioner. M kallas en linj¨ar differentialoperator (av
andra ordningen) eftersom M satisfierar (3).
Ekvationen
y 00 + a(x)y 0 + b(x)y = 0
(11)
kallas den till (10) h¨orande homogena ekvationen. (10) kallas inhomogena ekvationen
L˚
at yp vara en given l¨osning till (10). D˚
a ¨ar funktionen y l¨osning till (10) om och endast
om y ¨ar av formen
y = yh + yp ,
d¨ar funktionen yh ¨ar en l¨osning till motsvarande homogena ekvationen (11). L¨osningen yp kallas
partikul¨
arl¨osning.
Betrakta ekvationen
y 00 + ay 0 + by = 0
2
(12)
med konstanta (komplexa) koefficienterna a och b. L¨osningen till homogena ekvationen (12) ¨ar
av formen
y = C1 er1 x + C2 er2 x ,
r1 6= r2 ,
(13)
r1 = r2 = r,
(14)
eller
y = (C1 + C2 x)erx ,
d¨ar r1 och r2 ¨ar nollst¨allena till motsvarande karakteristiska polynomet
r2 + ar + b.
(15)
y 00 − 4y 0 + 3y = 0
(16)
r2 − 4r + 3
(17)
Exempel 1
Ekvationen
har karakteristiska polynomet
med nollst¨allena r1 = 1 och r2 = 3. Den fullst¨andiga l¨osningen till homogena ekvationen (16)
¨ar
y = C1 ex + C2 e3x .
(18)
y 00 + y = 0
(19)
r2 + 1
(20)
Ekvationen
har karakteristiska polynomet
med komplexa nollst¨allena r1 = i och r2 = −i. Den fullst¨andiga l¨osningen till (19) ¨ar
y = C1 eix + C2 e−ix = C1 cos x + iC1 sin x + C2 cos x − iC2 sin x = C˜1 cos x + C˜2 sin x.
(21)
Randv¨ardesproblemet f¨or linj¨ara differentialekvationen av andra ordningen y 00 + q(x)y = 0
skrivas som
(
y 00 + q(x)y = 0, y = y(x), a < x < b,
y(a) = y0 , y(b) = y1 ,
(22)
d¨ar q(x) ¨ar en given kontinuerlig funktion.
F¨or randv¨ardesproblem med icke-konstant koefficienten q m˚
aste man i allm¨anhet ber¨akna
en approximativ l¨osning.
3
Ibland kan man l¨osa (22) genom att hitta en partikul¨arl¨osning till ekvationen y 00 + q(x)y = 0
f¨or vissa funktioner q(x). Tydligen, ett par {y(x), q(x)} s˚
adant att
{y(x) = g(x),
q(x) = −
g 00 (x)
[g(x) 6= 0, a ≤ x ≤ b]}
g(x)
(23)
l¨oser y 00 + q(x)y = 0 och randv¨ardesproblemet (22) om g(x) satisfierar randvilkoren g(a) =
y0 , g(b) = y1 .
Exempel 2
L˚
at y(x) = P2 (x) = p0 x2 + p1 x + p2 vara ett polynom av andra ordningen. D˚
a l¨oser paret
{y(x) = p0 x2 + p1 x + p2 ,
q(x) = −
p0
x2
2p0
[P2 (x) 6= 0, a ≤ x ≤ b]}
+ p1 x + p2
(24)
ekvationen y 00 + q(x)y = 0. Satisfiera randvillkoren
p0 a2 + p1 a + p2 = y0 ,
p0 b2 + p1 b + p2 = y1
Det ¨ar ett underbest¨amt ekvationssystem eftersom a, b, y0 och y1 ¨ar givna tal (indata) och
antalet ekvationer (tv˚
a) ¨ar mindre ¨an antalet obekanta p0 , p1 , p2 . D˚
a f˚
ar man



p0 =
p1 =


p2 =
var,
− p0 (b + a),
y0 − p1 a − p0 a2 ,
y1 −y0
b−a
Om y0 = y1 och a = 0, d˚
a f˚
ar man



p0 = var,
p = −p0 b,
 1

p2 = y0 ,
och polynomet (l¨osningen) y(x) = P2 (x) = p0 x2 − p0 bx + y0 = p0 x(x − b) + y0 .
Randv¨
ardesproblemet f¨
or v¨
armekvationen
Antag att vi har en stav av l¨angden L med konstant v¨armeledningsf¨orm˚
aga. Antag vidare
att staven ¨ar isolerad i l¨angriktingen och att ¨andarna h˚
als vid konstant temperatur 0. L˚
at
u = u(x, t) beteckna stavens temperatur i punkten x, 0 < x < L, och tidpunkten t ≥ 0. D˚
a
satisfierar u = u(x, t) randv¨ardesproblemet f¨or en-dimensionella v¨
armekvationen
(
ut = c2 uxx , u = u(x, t), t > 0, 0 < x < L,
u(x, 0) = f (x), u(0, t) = 0,
u(L, t) = 0, t ≥ 0.
(25)
d¨ar c = const och f (x) ¨ar begynnelsetemperatur; de f¨orsta och andra partiella derivatorna
ut =
∂u
,
∂t
uxx =
∂ 2u
.
∂x2
Vi s¨oker en l¨osning till randv¨ardesproblemet (25) i formen
u(x, t) = F (x)G(t);
4
dvs vi g¨or det med hj¨alp av Fouriermetod och variabelseparation.
Ins¨attning av ansatsen
∂u
∂ 2u
= F G0 ,
= F 00 G
∂t
∂x2
i v¨armekvationen ut = c2 uxx ger ekvationen
F G0 = c2 F 00 G,
som har separabla variabler x och t och skrivas p˚
a formen
G0
F 00
=
= const = −p2 ,
c2 G
F
och vi f˚
ar differentialekvationer f¨or G(t) och F (x)
G0 + c2 p2 G = 0,
F 00 + p2 F = 0.
Satisfiera randvillkoren:
u(0, t) = F (0)G(t) = 0,
u(L, t) = F (L)G(t) = 0, t ≥ 0.
Det g¨aller
F (0) = 0,
F (L) = 0.
Den fullst¨andiga l¨osningen till ordin¨ara differentialekvationen av andra ordningen f¨or F ¨ar
F = A cos px + B sin px.
D˚
a kan man anv¨anda randvillkoren i (25) och skriva ekvationer f¨or att f˚
a A, B och p
F (0) = 0 :
A = 0;
F (L) = 0 :
B sin pL = 0
sin pL = 0 (B 6= 0),
nπ
pL = nπ, p = pn =
(n = 1, 2, . . .).
L
nπ
F = Fn = sin pn x = sin
x (n = 1, 2, . . .).
L
F¨or G f˚
ar vi en ordin¨ar differentialekvation av f¨orsta ordningen
G0 + λ2n G = 0,
λn =
cnπ
.
L
Den fullst¨andiga l¨osningen blir
2
G(t) = Gn (t) = Bn e−λn t
(n = 1, 2, . . .).
Nu kan man skriva en l¨osning u = un = un (x, t) till ekvationen
∂u
∂ 2u
= c2 2 , 0 < x < L
∂t
∂x
5
som satisfierar randvillkoren
u(0, t) = 0,
u(L, t) = 0, t ≥ 0,
2
un (x, t) = Fn (x)Gn (t) = Bn e−λn t sin
nπ
x (n = 1, 2, . . .).
L
un (x, t) kallas egenfunktioner och talen
cnπ
L
λn =
kallas egenv¨arden.
Nu kan man l¨osa randv¨ardesproblemet (25) genom att skriva
u(x, t) =
∞
X
un (x, t) =
n=1
∞
X
2
Bn e−λn t sin
n=1
nπ
x.
L
F¨or att f˚
a obekanta koefficienterna Bn , satisfierar man begynnelsevillkoren
u(x, 0) =
∞
X
Bn sin
n=1
nπ
x = f (x).
L
D˚
a
2ZL
nπ
f (x) sin
xdx,
L 0
L
blir Fourierskoefficienterna till funktionen f (x).
Bn =
n = 1, 2, . . . ,
Exempel 3 Sinus-begynnelsetemperatur
L¨os randv¨ardesproblemet (25) med sinus-begynnelsetemperatur
f (x) = 100 sin
π
x.
80
L¨
osning. Satisfiera begynnelsevillkoren
u(x, 0) =
∞
X
Bn sin
n=1
nπ
π
x = f (x) = 100 sin x.
80
80
λn = λ1 =
cπ
.
80
D˚
a blir l¨osningen
u(x, t) =
∞
X
2
un (x, t) = B1 e−λ1 t sin
n=1
π
π
2
x = 100e−λ1 t sin x.
80
80
Exempel 4 Triangel-begynnelsetemperatur
L¨os randv¨ardesproblemet (25) med triangel-begynnelsetemperatur
f (x) =
(
x
if 0 < x < L/2,
L − x if L/2 < x < L.
6
L¨
osning. Satisfiera begynnelsevillkoren
u(x, 0) =
∞
X
Bn sin
n=1
nπ
x = f (x).
L
Fourierkoefficienterna till en udda-periodisk f (x) [f (−x) = −f (x)] ¨ar
!
Z L
2ZL
nπ
2 Z L/2
nπ
nπ
xdx =
xdx +
(L − x) sin
xdx =
Bn =
f (x) sin
x sin
L 0
L
L 0
L
L
L/2
2
2 Z L/2
nπ L/2
nπ
− x cos
x +
cos
xdx−
nπ
L 0
nπ 0
L
2 ZL
nπ
2
nπ L
− (L − x) cos
x −
cos
xdx =
nπ
L L/2 nπ L/2
L
−
L
nπ
2L
nπ
L
nπ
2L
nπ
cos
+ 2 2 sin
+
cos
+ 2 2 sin
=
πn
2
π n
2
πn
2
π n
2
(26)
4L
nπ
sin
=
n2 π 2
2
4L
4 (−1)l+1
l+1
(−1)
=
=
(2l − 1)2 π 2
π (2l − 1)2
(
4L
n2 π 2
− n4L
2 π2
if n = 1, 5, 9, . . .,
if n = 3, 7, 11, . . ..
Vi f˚
ar egenv¨ardena
λn =
cnπ
,
L
och l¨osningen
u(x, t) =
∞
X
un (x, t) =
n=1
cπ 2
3cπ 2
π
3π
4L
1
sin xe−( L ) t − sin xe−( L ) t + . . . .
2
π
L
9
L
Numerisk l¨
osning till en-dimensionella v¨
armekvationen i MATLAB
Vi ska skapa MATLABsatser som l¨oser randv¨ardesproblemet (25). Skriv en-dimensionella
v¨armekvationen som
K
∂u
∂ 2u
=
,
2
∂x
∂t
K = c2 , t > 0, 0 < x < L.
u = u(x, t),
7
F¨or att approximera funktionen u(x, t), x ∈ [0, L], t ∈ [0, T ], och v¨armekvationen och l¨osa (25)
numeriskt genom att ers¨atta den med en differensapproximation delar man in intervallet [0, L] i
n (mindre) intervall och intervallet [0, T ] i m intervall genom att ber¨akna funktionen f¨or (n − 1)
L
och (m − 1) stycken x- och t-v¨arden xi = ih och t = jk likformigt f¨ordelade med avst˚
and h =
n
T
och k = :
m
xi = ih,
h=
L
,
n
i = 0, 1, 2 . . . , n − 1, n,
(27)
x0 = 0 < x1 < x2 < . . . < xn−1 < xn = L.
(28)
tj = jk,
j = 0, 1, 2 . . . , m − 1, m,
k=
T
.
m
(29)
Nu kan vi approximera v¨armekvationen med differenserna
ui−1,j+1 − 2ui,j+1 + ui+1,j+1
K
h2
−ui,j + ui,j+1
=
.
k
(30)
d¨ar ui,j (i = 0, 1, 2, . . . , n − 1 j = 0, 1, 2, . . . , m − 1) ¨ar obekanta [approximationen till u(xi , tj )].
D˚
a skrivas (30) som
(1 + 2α)ui,j+1 − α(ui+,j+1 + ui−1,j+1 ) = ui,j ,
i = 0, 1, . . . , n,
j = 0, 1, . . . ,
α = Kk/h2 .(31)
(31) skrivas som ett ekvationssystemet med symmetriska trediagonala matrisen A















b −α 0 . . . 0
−α b −α . . . 0
0 −α b . . . 0
.
.
. ... .
.
.
. ... .
.
.
. ... .
0
0
0 . . . −α
0
0
0 ... b















u1,j+1
u2,j+1
u3,j+1
.
.
.
un−2,j+1
un−1,j+1
















=














u1,j+1 + αu0
u2,j
u3,j
.
.
.
un−2,j
un−1,j + αun















(32)
d¨ar
u0 = u(0, t) = lowb = 0,
un = u(L, t) = hib = 0.
L¨osningen till systemet (32) ger u1 , u2 , . . ., un−1 i tidpunkten tj+1 genom u1 , u2 , . . ., un−1 i
tidpunkten tj (t0 = 0, j = 0, 1, 2, . . .). Man skapar begynnelsetemperatur f (x) som en (n + 1)dimensionell radvektor
f = [f (xi )] = [f (x0 ) = lowb = 0, f (x1 ), . . . , f (xn−1 ), f (xn ) = hib = 0].
MATLABfunktionen
heat
utf¨or den h¨ar metoden som kallas implicit difference solution.
function u = heat(nx,hx,nt,ht,init,lowb,hib,K)
%nx, hx are the number and size of the x-grid sections, index i (nx ¿ 2)
%nt, ht are the number and size of the t-grid sections, index j (nt ¿ 2)
%init is a row vector of nx + 1 initial values of function f(x)
8
%lowb, hib are boundaries at low and high values of x at the ends of x interval
%solution u is determined as an (nt + 1 x nx + 1) matrix
alpha = K∗ht/(hx∗hx);
a = zeros(nx - 1, nx - 1); u = zeros(nt + 1, nx + 1);
u(:,1) = lowb∗ones(nt + 1, 1);
%the first column of solution matrix u consisting of nt+1 elts. is set
u(:,nx + 1) = hib∗ones(nt + 1, 1);
%the last column of solution matrix u consisting of nt+1 elts. is set
u(1,:) = init;
%the first row of solution matrix u consisting of nx+1 elts. is set equal to initial row data
vector init
%start formation of tri-diagonal matrix a of the equation system
a(1,1) = 1 + 2∗alpha; a(1,2) = -alpha;
for i = 2:nx-2
a(i,i) = 1 + 2∗alpha; a(i,i-1) = -alpha; a(i,i+1) = -alpha;
end
a(nx-1,nx-2) = -alpha; a(nx-1,nx-1) = 1 + 2∗alpha;
%tri-diagonal matrix a is formed
%start formation of right-hand side column vector bb of the equation system
bb(1,1) = init(2) + init(1)∗alpha;
for i = 2:nx-2, bb(i,1) = init(i+1); end
bb(nx-1,1) = init(nx) + init(nx+1)∗alpha;
%right-hand side column vector bb of the equation system is formed
[L,U] = lu(a);
%the equation system is prepared for solution
for j = 2:nt+1
y=L\bb; x=U\y;
%the equation system is solved at the current time step in index j
u(j,2:nx) = x’; bb = x;
bb(1,1) = bb(1,1) + lowb∗alpha;
bb(nx-1,1) = bb(nx-1,1) + hib∗alpha;
%the right-hand side of the equation system is prepared for solution at the next step in j
end
%the solution is determined as an (nt + 1 x nx + 1) matrix u
Exempel 5
Ber¨akna temperatur i en tegelsten v¨agg med K = 5 · 10−7 m2 /s, L = 0.3 m och begynnelsetemperatur 100 deg C. Antag att v¨aggs¨andarna h˚
als vid konstant temperatur 20 deg C. Vi vill rita
v¨aggstemperatur under 22000 sekunder (366.67 min), tidintervallet k ¨ar lika med 440 sec (7.33
min); d˚
a blir det 50 t-tidintervall.
Vi ska anv¨anda 15 x-intervall. D˚
a, hx = 0.02, ht = 440, och
α = Kk/h2 = 5 · 10−7 · 440 · (0.3/15)−1 ≈ 0.55.
Anv¨and MATLABfunktionen heat och skriv MATLABsatser som utf¨or
l¨osa motsvarande randv¨ardesproblemet (25).
9
heat f¨or att
K=5e-7; hx=0.02; nx=15; hx=0.02; ht=440; nt=50;
init = 100∗ones(1,nx + 1);
lowb = 20; hib = 20;
%the initial data are set
u = heat(nx,hx,nt,nt,init,lowb,hib,K);
%the solution is determined as (nt + 1 x nx + 1) = (51 x 16) matrix u
surfl(u)
axis([0 16 0 50 0 120])
view([-217 30])
xlabel(’x - node nos.’); ylabel(’time - node nos.’);
zlabel(’temperature’);
Problem 1. Kalla v¨
aggs¨
andarna
Ber¨akna temperatur i en tegelsten v¨agg med K = 5 · 10−7 m2 /s, L = 0.3 m och begynnelsetemperatur 100 deg C. Anv¨and MATLABfunktionen heat och skriv MATLABsatser som utf¨or
heat f¨or att l¨osa motsvarande randv¨ardesproblemet (25). Antag att v¨aggs¨andarna h˚
als vid
konstant temperatur 0 deg C.
Problem 2. K¨
allasmodell
Anv¨and MATLABfunktionen heat och skriv MATLABsatser som utf¨or heat f¨or att l¨osa
motsvarande randv¨ardesproblemet (25). som ber¨aknar temperatur i en tegelsten v¨agg med
K = 5 · 10−7 m2 /s, L = 0.3 m och begynnelsetemperatur 100 deg C; valj Tns = Tsource 6= 100,
d¨ar 1 < ns < nx.
Randv¨
ardesproblemet f¨
or den tv˚
a-dimensionella hyperboliska ekvationen
Antag att vi har en rektangul¨ar platta av l¨angden a och bredden b. Antag vidare att
plattens sidor h˚
als fast (or¨orlig). L˚
at u = u(x, y, t) beteckna h¨ojden av plattens sm˚
a vertikala
vibrationer i punkten (x, y), 0 < x < a, 0 < y < b, och tidpunkten t ≥ 0. D˚
a satisfierar
u = u(x, y, t) randv¨ardesproblemet f¨or tv˚
a-dimensionella hyperboliska ekvationen








H¨ar,







utt = c2 ∆u, u = u(x, y, t), t > 0, 0 < x < a, 0 < y < b,
u(x, y, 0) = ϕ(x, y)
ut (x, y, 0) = ψ(x, y)
u(0, y, t) = 0,
u(a, y, t) = 0,
u(x, 0, t) = 0,
u(x, b, t) = 0, t ≥ 0.
utt =
∂ 2u
,
∂t2
∆u =
(33)
∂ 2u ∂ 2u
+
.
∂x2 ∂y 2
och ϕ(x, y) och ψ(x, y) ¨ar givna funktioner som definierar begynnelsepositionen och begynnelsehastigheten i punkten (x, y).
Vi s¨oker en l¨osning till randv¨ardesproblemet (33) i formen
u(x, y, t) = v(x, y)T (t);
10
dvs vi g¨or det med hj¨alp av Fouriermetod och variabelseparation. Ins¨attning av ansatsen
u(x, y, t) = v(x, y)T (t)
ger ekvationen f¨or T
T 00 + a2 λT = 0
,
och ekvationen samst randv¨ardesproblemet f¨or v(x, y)



vxx + vyy + λv = 0, 0 < x < a, 0 < y < b
v(0, y) = 0,
v(a, y) = 0


v(x, 0) = 0,
v(x, b) = 0
(34)
Vi ska ocks˚
a s¨oka en l¨osning till randv¨ardesproblemet (34) i formen med separabla variabler
v(x, y) = X(x)Y (y).
Ins¨attning av den h¨ar ansatsen i (34) ger ekvationer och randv¨ardesproblem f¨or X(x) och Y (y)
(
X 00 + νX = 0, 0 < x < a
X(0) = 0,
X(a) = 0
(35)
(
Y 00 + µY = 0, 0 < y < b
Y (0) = 0,
Y (b) = 0
(36)
H¨ar µ och ν ¨ar parameter (tal) s˚
adana att µ + ν = λ.
Randv¨ardesproblemen (35) och (36) kallas Sturm–Liouville eigenv¨
ardesproblem.
Vi har
nπ
nπ
X(x) = Xn (x) = sin
x,
νn = ( )2 ;
a
a
mπ 2
mπ
y,
µm = (
).
Y (y) = Ym (y) = sin
b
b
λ=(
nπ 2
mπ 2
) +(
)
a
b
vn,m = An,m sin
nπ
mπ
x sin
y.
a
b
Koefficienterna An,m f˚
ar man genom att anv¨anda ortogonalitetsvillkor
Z aZ b
0
0
2
dxdy = A2n,m
vn,m
Z a
sin2
0
D˚
a
An,m =
och
vn,m (x, y) =
s
Z b
nπ
mπ
xdx
sin2
ydy = 1.
a
b
0
s
4
,
ab
4
nπ
mπ
sin
x sin
y.
ab
a
b
11
L¨osningen till (33) ¨ar
u(x, y, t) =
∞ X
∞
X
q
q
(Cn,m cos λn,m at + Dn,m sin λn,m at)vn,m (x, y),
m=1 n=1
d¨ar
Cn,m =
Z aZ b
0
=
s0
4 Z aZ b
nπ
mπ
ϕ(x, y) sin
x sin
ydxdy,
ab 0 0
a
b
1
Dn,m = q
ϕ(x, y)vn,m (x, y)dxdy =
a2 λn,m
s
4 Z aZ b
nπ
mπ
ψ(x, y) sin
x sin
ydxdy.
ab 0 0
a
b
Randv¨
ardesproblemet f¨
or tv˚
a-dimensionella v¨
armekvationen
Antag att vi har en rektangul¨ar platta av l¨angden a och bredden b med konstant v¨arme
ledningsf¨orm˚
aga. Antag vidare att platten ¨ar isolerad och att sidorna h˚
als vid konstant temperatur 0. L˚
at u = u(x, y, t) beteckna plattens temperatur i punkten (x, y), 0 < x < a, 0 < y < b,
och tidpunkten t ≥ 0. D˚
a satisfierar u = u(x, y, t) randv¨ardesproblemet f¨or tv˚
a-dimensionella
v¨
armekvationen





H¨ar,
ut = c2 ∆u, u = u(x, y, t), t > 0, 0 < x < a, 0 < y < b,
u(x, y, 0) = ϕ(x, y),

u(0, y, t) = 0,
u(a, y, t) = 0,



u(x, 0, t) = 0,
u(x, b, t) = 0, t ≥ 0.
ut =
∂u
,
∂t
∆u =
(37)
∂ 2u ∂ 2u
+
,
∂x2 ∂y 2
och ϕ(x, y) ¨ar en given funktion som definierar begynnelsetemperatur i punkten (x, y).
Problem 3
L¨os randv¨ardesproblemet (25) f¨or en-dimensionella v¨armekvationen med givet c och begynnelsetemperatur
f (x) = sin 0.1πx.
L¨
osning. Vi har
f (x) = sin 0.1πx = sin
π
x.
10
D˚
a L = 10. Satisfiera begynnelsevillkoren
u(x, 0) =
∞
X
n=1
Bn sin
nπ
x = f (x) = sin 0.1πx.
L
λn = λ1 =
12
cπ
.
10
Nu kan man f˚
a l¨osningen
u(x, t) =
∞
X
2
un (x, t) = B1 e−λ1 t sin
n=1
π
π
2
x = e−λ1 t sin x.
10
10
Randv¨
ardesproblemet f¨
or stationera tv˚
a-dimensionella v¨
armekvationen
Antag att vi har en rektangul¨ar platta av l¨angden a och bredden b med konstant v¨armelednings
f¨orm˚
aga. Antag vidare att platten ¨ar isolerad och att sidorna h˚
als vid konstant eller ickekonstant temperatur som beskrivas med fyra givna funktioner f1 (y), f2 (y), 0 < y < b, och f3 (x),
f4 (x), 0 < x < a. L˚
at u = u(x, y, t) beteckna plattens temperatur i punkten (x, y), 0 < x < a,
0 < y < b, och tidpunkten t ≥ 0. Om f¨ors¨oket p˚
ag˚
at l¨ange kommer plattens temperatur att
vara stationar, dvs oberoende av tiden. D˚
a satisfierar u = u(x, y) randv¨ardesproblemet f¨or
tv˚
a-dimensionella Laplace ekvationen i rektangeln med h¨ornen i punkterna (0, 0), (a, 0), (a, b)
och (0, b):



∆u = 0, u = u(x, y), 0 < x < a, 0 < y < b,
u(0, y) = f1 (y),
u(a, y) = f2 (y),


u(x, 0) = f3 (x),
u(x, b) = f4 (x).
(38)
L¨agg m¨arke till att motsvarande homogena randv¨ardesproblemet



∆u = 0, u = u(x, y), 0 < x < a, 0 < y < b,
u(0, y) = 0,
u(a, y) = 0,


u(x, 0) = 0,
u(x, b) = 0
(39)
har triviala l¨osningen.
Numerisk l¨
osning till Laplace ekvationen i en rektangel
Skriv MATLABfunktionen ellipgen som utf¨or numerisk l¨osning till randv¨ardesproblemet
(38) i rektangeln med h¨ornen i punkterna (x0,y0), (xn,y0), (xn,yn), (x0,yn).
function [a,om] = ellipgen(nx,hx,ny,hy,bx0,bxn,by0,byn)
%nx, hx are the number and size of the x-grid sections
%ny, hy are the number and size of the y-grid sections
%bx0, bxn are vector of boundary conditions at x0, xn, with (ny+1) elements each,
%beginning at y0 and corresponding to f1 (y), f2 (y)
%by0, byn are vector of boundary conditions at y0, yn, with (nx+1) elements each
%beginning at x0 and corresponding to f3 (x), f4 (x)
%solution is determined as an (ny + 1 x nx + 1) matrix a
nmax = (nx - 1)∗(ny - 1); r=hy/hx;
a = zeros(ny+1, nx+1); p = zeros(ny + 1, nx + 1);
gg = zeros(ny+1, nx+1); ff = zeros(ny+1, nx+1);
test = 1;
if bx0 == zeros(1, ny + 1), test = test + 1; end
if bxn == zeros(1, ny + 1), test = test + 1; end
if by0 == zeros(1, nx + 1), test = test + 1; end
13
if byn == zeros(1, nx + 1), test = test + 1; end
if test == 5
disp(’WARNING – problem has trivial solution, z = 0.’)
break
end
bx0 = bx0(1, ny+1:-1:1); bxn = bxn(1, ny+1:-1:1);
a(1,:) = byn; a(ny+1,:) = by0;
a(:,1) = bx0’; a(:,nx+1) = bxn’;
for i = 2:ny
for j = 2:nx
nn = (i - 2)∗(nx - 1) + (j - 1);
q(nn,1) = i; q(nn,2) = j; p(i,j) = nn;
end
end
c = zeros(nmax, nmax); e = zeros(nmax, 1); om = zeros(nmax, 1);
g = zeros(nmax, 1);
for i = 2:ny
for j = 2:nx
nn = p(i,j); c(nn,nn)=-(2+2∗r∗r); e(nn) = hy∗hy∗gg(i,j);
g(nn) = g(nn) + hy∗hy∗ff(i,j);
if p(i+1,j)∼=0
np=p(i+1,j); c(nn,np)=1;
else
g(nn) = g(nn) - by0(j);
end
if p(i-1,j)∼=0
np=p(i-1,j); c(nn,np)=1;
else
g(nn) = g(nn) - byn(j);
end
if p(i,j+1)∼=0
np=p(i,j+1); c(nn,np)=r∗r;
else
g(nn) = g(nn) - r∗r∗bxn(i);
end
if p(i,j-1)∼=0
np=p(i,j-1); c(nn,np)=r∗r;
else
g(nn) = g(nn) - r∗r∗bx0(i);
end
end
end
c = c + diag(e); z=c\g;
for nn = 1:nmax
i=q(nn,1); j=q(nn,21);
a(i,j) = z(nn);
end
14
Exempel 6
Ber¨akna temperaturen i en rektangul¨ar platta med konstant v¨armelednings
f¨orm˚
aga. Platten ¨ar isolerad och sidorna h˚
als vid konstant eller icke-konstant temperatur
som beskrivas med fyra givna radvektorer bx0, bxn, by0, byn. I kommandofilen utnyttjar vi
ellipgen med 6 punkter p˚
a varje sida och ritar contour plots:
nx=6; ny=6; hx=0.5; hy = 0.3333;
bx0 = [0 33.33 66.67 100 133.33 166.67 200];
bxn = [0 83.33 166.67 250 333.33 416.67 500];
by0 = [0 0 0 0 0 0 0];
byn = [200 208.33 233.33 275 333.33 408.33 500];
a = ellipgen(nx,hx,ny,hy,bx0,bxn,by0,byn);
aa = flipud(a);
contour(aa)
xlabel(’node numbers in x direction’);
ylabel(’node numbers in y direction’);
Exempel 7. K¨
allasmodell
Nu ritar vi temperaturen i en rektangul¨ar platta n¨ar sidorna h˚
als vid icke-konstant temperatur
som beskrivas med 21-dimensionella radvektorer bx0, bxn, by0, byn.
nx=20; ny=20;
hx=1/20; hy = 1/20;
bx0 = [0 0 0 0.01 0 0 0 0 0 0 0.01 0 0 0 0 0 0 0.01 0 0 0];
bxn = [0 0 0 0.01 0 0 0 0 0 0 0.01 0 0 0 0 0 0 0.01 0 0 0];
by0 = [0 0 0 0.1 0 0 0 0 0 0 0.1 0 0 0 0 0 0 0.01 0 0 0];
byn = [0 0 0 0 0 0 0 0 0 0 0.01 0 0 0 0 0 0 0.01 0 0 0];
a = ellipgen(nx,hx,ny,hy,bx0,bxn,by0,byn);
surfl(a)
axis([1 20 1 20 0 0.2])
xlabel(’x-node nos.’);
ylabel(’y-node nos.’);
zlabel(’Temperature’);
Problem 4
Anv¨and MATLAB funktionen ellipgen f¨or att ber¨akna och rita temperaturen i en rektangul¨ar platta.
Anv¨and MATLABfunktionen heat och skriv MATLABsatser som utf¨or heat f¨or att
l¨osa motsvarande randv¨ardesproblemet (25).
(i) Skriv funktionsfiler som ber¨aknar randvillkoren f1 (y), f2 (y), f3 (x), och f4 (x).
(ii) Skapa x− och y-variabel radvektorer x och y med nx+1 och ny+1 element (anv¨and
MATLABsatsen linspace).
(iii) Skapa x− och y-randvillkor radvektorer bx0, bxn, by0, byn genom att utnyttja
funktionsfiler som ber¨aknar randvillkoren fj och radvektorer x, y.
15
(iv) utf¨or ellipgen och rita temperaturen (anv¨and MATLABsatsen surfl
passande skala).
(v) Valj ett randvillkor s˚
adant att
2
f3 (x) = Ae−Q(x−P ) ,
A, Q > 0,
0<P <a
och valj en
(40)
och fj = 0, j = 1, 2, 4. Det ¨ar en k¨allasmodell d¨ar A ¨ar k¨allasamplitud, Q ¨ar k¨allaskvalitet och
P ¨ar k¨allasposition. Skriv en funktionsfil som ber¨aknar temperaturen f¨or givna A, Q och P .
Rita en graf som visar hur temperaturen u(xic , yic ) i plattens centrala punkt (xic , yic ) ber˚
ar p˚
a
A, Q eller P .
16