6 Εξίσωση θερμότητας: Πεπερασμένα στοιχεία

6.3 Μέθοδος των Crank–Nicolson

Στη συνέχεια, θα θεωρήσουμε ένα πλήρως διακριτό σχήμα, διακριτοποιώντας και ως προς τον χρόνο το ημιδιακριτό σχήμα (6.4) με την μέθοδο των Crank–Nicolson. Θέτουμε και πάλι k=T/M και tj=jk, j=0,,M, μια διαμέριση του [0,T].

Θεωρούμε λοιπόν τις προσεγγίσεις UnVh, n=1,,M, τέτοιες ώστε

(Un-Un-1k,χ)+((Un-1/2),χ)=0,χVh,U0=gh, (6.30)

όπου ghVh μια προσέγγιση της g και Un-1/2=12(Un+Un-1). Η (6.30) γράφεται και ως εξής

(Un,χ)+k2((Un),χ)=(Un-1,χ)-k2((Un-1),χ),χVhU0=gh. (6.31)

Με ανάλογο τρόπο όπως και για την πεπλεγμένη μέθοδο του Euler στην προηγούμενη παράγραφο, η (6.31) είναι ισοδύναμη με ένα γραμμικό σύστημα της μορφής

(+k2𝒮)αn=(-k2𝒮)αn-1,

όπου και 𝒮 είναι οι πίνακες της (6.6) και αj διανύσματα με συνιστώσες αj=(α1j,,αNj)T, j=0,,N. Εύκολα βλέπουμε ότι ο πίνακας +k2𝒮 είναι συμμετρικός και θετικά ορισμένος και, άρα, αντιστρέφεται. Επομένως,

αn=(+k2𝒮)-1(-k2𝒮)αn-1. (6.32)

Συνεπώς, επειδή γνωρίζουμε το διάνυσμα α0, χρησιμοποιώντας την (6.32) μπορούμε να προσδιορίσουμε αναδρομικά τα αn και κατ’επέκταση την προσέγγιση Un της u(,tn).

Στη συνέχεια, θα δείξουμε την ευστάθεια της μεθόδου (6.30).

Θεώρημα 6.5.

Έστω ότι οι UnVh, n=0,,M, ικανοποιούν την (6.30). Τότε

max0nMUnU0. (6.33)
Απόδειξη.

Επιλέγουμε χ=Un-1/2 στην (6.30). Οπότε

(Un-Un-1k,Un-1/2)+(Un-1/2)2=0. (6.34)

Επομένως, έχουμε ότι

(Un-Un-1,Un+Un-1)0,

από όπου εύκολα βλέπουμε ότι

Un2Un-12U02,

δηλαδή τη ζητούμενη ανισότητα. ∎

Θεώρημα 6.6.

Έστω uC4([0,L]×[0,T]) η λύση της (6.1) με gC2[0,L] και UnVh, n=0,,M, ικανοποιούν την (6.30) με U0=Rhg. Τότε, υπάρχει σταθερά C, ανεξάρτητη των k και h, τέτοια ώστε

max0nMUn-u(,tn)C(k2+h2). (6.35)
Απόδειξη.

Θέτουμε και πάλι ρn=W(,tn)-u(,tn) και ϑn=Un-W(,tn), οπότε Un-u(,tn)=ϑn+ρn. Λόγω της (6.25) αρκεί να εκτιμήσουμε την ϑn. Μπορούμε να δούμε ότι το ϑn ικανοποιεί την ακόλουθη σχέση

(ϑn-ϑn-1k,χ)+((ϑn-1/2),χ)=-(ωn,χ),χVh, (6.36)

όπου ωn=ω1n+ω2n+ω3n και

ω1n =ρn-ρn-1k,ω2n=un-un-1k-ut(,tn-1/2),
ω3n =uxx(,tn-1/2)-uxxn-1/2.

Πράγματι έχουμε

(ϑn-ϑn-1k,χ)+((ϑn-1/2),χ)
 =-(Rhun-Rhun-1k,χ)-((Rhun-1/2),χ)
 =-(ρn-ρn-1k,χ)-(un-un-1k,χ)-(uxn-1/2,χ)
 =-(ρn-ρn-1k,χ)-(un-un-1k,χ)+(uxxn-1/2,χ)
 =-(ρn-ρn-1k,χ)-(un-un-1k,χ)+(ut(,tn-1/2),χ)
-(ut(,tn-1/2),χ)+(uxxn-1/2,χ)
 =-(ρn-ρn-1k,χ)-(un-un-1k,χ)+(ut(,tn-1/2),χ)
-(uxx(,tn-1/2),χ)+(uxxn-1/2,χ)
 =-(ω1n+ω2n+ω3n,χ).

Από την (6.27) έχουμε ότι ω1nCh2. Επίσης, λόγω του αναπτύγματος Taylor έχουμε

ω2n =un-un-1k-ut(,tn-1/2)
=(2k)-1tn-1tn-1/2(t-tn-1)2uttt(t)dt+(2k)-1tn-1/2tn(t-tn)2uttt(t)dt,

από όπου προκύπτει ότι ω2nCk2. Τέλος, χρησιμοποιώντας και πάλι το ανάπτυγμα Taylor, έχουμε

2ω3n =2(uxx(,tn-1/2)-(uxx)n-1/2)
=-tn-1tn-1/2(t-tn-1)uxxtt(t)dt+tn-1/2tn(t-tn)uxxtt(t)dt.

Συνεπώς ω3nCk2. Συνδυάζοντας τις παραπάνω εκτιμήσεις, λαμβάνουμε

ωnC(k2+h2). (6.37)

Επιλέγοντας τώρα χ=ϑn-1/2 στην (6.36), παίρνουμε

ϑn2-ϑn-12+2k(ϑn-1/2)2 =-k(ωn,ϑn+ϑn-1)
kωn(ϑn+ϑn-1).

Οπότε

(ϑn-ϑn-1)(ϑn+ϑn-1)kωn(ϑn+ϑn-1),

από την οποία παίρνουμε

(ϑnϑn-1+kωn.

Συνεπώς, έχουμε

ϑnϑ0+kj=1nωjϑ0+Mkmax0jMωj.

Επομένως, λόγω των (6.25), (6.37) και του γεγονότος ότι ϑ(0)=0, έχουμε ότι ισχύει η επιθυμητή σχέση (6.35). ∎