author:

  • 'Adrian E. Feiguin' title: 'Computational Physics' ...

Ordinary differential equations

Let’s consider a simple 1st order equation: $$\frac{dy}{dx}=f(x,y)$$

To solve this equation with a computer we need to discretize the differences: we have to convert the differential equation into a “finite differences” equation. The simplest solution is Euler’s method.

Euler’s method

Supouse that at a point $x_0$, the function $f$ has a value $y_0$. We want to find the approximate value of $y$ in a point $x_1$ close to $x_0$, $x_1=x_0+\Delta x$, with $\Delta x$ small. We assume that $f$, the rate of change of $y$, is constant in this interval $\Delta x$. Therefore we find: $$\begin{eqnarray} && dx \approx \Delta x &=&x_1-x_0, \\ && dy \approx \Delta y &=&y_1-y_0,\end{eqnarray}$$ with $y_1=y(x_1)=y(x_0+\Delta x)$. Then we re-write the differential equation in terms of discrete differences as: $$\frac{\Delta y}{\Delta x}=f(x,y)$$ or $$\Delta y = f(x,y)\Delta x$$ and approximate the value of $y_1$ as $$y_1=y_0+f(x_0,y_0)(x_1-x_0)$$ We can generalize this formula to find the value of $y$ at $x_2=x_1+\Delta x$ as $$y_{2}=y_1+f(x_1,y_1)\Delta x,$$ or in the general case: $$y_{n+1}=y_n+f(x_n,y_n)\Delta x$$

This is a good approximation as long as $\Delta x$ is “small”. What is small? Depends on the problem, but it is basically defined by the “rate of change”, or “smoothness” of $f$. $f(x)$ has to behave smoothly and without rapid variations in the interval $\Delta x$.

Notice that Euler’s method is equivalent to a 1st order Taylor expansion about the point $x_0$. The “local error” calculating $x_1$ is then $O(\Delta x^2)$. If we use the method $N$ times to calculate $N$ consecutive points, the propagated “global” error will be $NO(\Delta x^2)\approx O(\Delta x)$. This error decreases linearly with decreasing step, so we need to halve the step size to reduce the error in half. The numerical work for each step consists of a single evaluation of $f$.

Exercise 1.1: Newton’s law of cooling

If the temperature difference between an object and its surroundings is small, the rate of change of the temperature of the object is proportional to the temperature difference: $$\frac{dT}{dt}=-r(T-T_s),$$ where $T$ is the temperature of the body, $T_s$ is the temperature of the environment, and $r$ is a “cooling constant” that depends on the heat transfer mechanism, the contact area with the environment and the thermal properties of the body. The minus sign appears because if $T>T_s$, the temperature must decrease.

Write a program to calculate the temperature of a body at a time $t$, given the cooling constant $r$ and the temperature of the body at time $t=0$. Plot the results for $r=0.1\frac{1}{min}$; $T_s=83^{\circ} C$, $T_0=10^{\circ}$ using different intervals $\Delta t$ and compare with exact (analytical) results.


In [1]:
T0 = 10.   # initial temperature
Ts = 83.   # temp. of the environment
r = 0.1    # cooling rate
dt = 0.05  # time step
tmax = 60. # maximum time
nsteps = int(tmax/dt)  # number of steps

In [13]:
T = T0
for i in range(1,nsteps+1):
    new_T = T - r*(T-Ts)*dt
    T = new_T
    print ('{:20.18f}  {:20.18f}  {:20.18f}'.format(i,i*dt, T))
    # we can also do t = t - r*(t-ts)*dt


1.000000000000000000  0.062500000000000000  10.456250000000000711
2.000000000000000000  0.125000000000000000  10.909648437500001350
3.000000000000000000  0.187500000000000000  11.360213134765626108
4.000000000000000000  0.250000000000000000  11.807961802673341722
5.000000000000000000  0.312500000000000000  12.252912041406633037
6.000000000000000000  0.375000000000000000  12.695081341147840703
7.000000000000000000  0.437500000000000000  13.134487082765666344
8.000000000000000000  0.500000000000000000  13.571146538498380707
9.000000000000000000  0.562500000000000000  14.005076872632765017
10.000000000000000000  0.625000000000000000  14.436295142178810380
11.000000000000000000  0.687500000000000000  14.864818297540193015
12.000000000000000000  0.750000000000000000  15.290663183180566165
13.000000000000000000  0.812500000000000000  15.713846538285688226
14.000000000000000000  0.875000000000000000  16.134384997421403085
15.000000000000000000  0.937500000000000000  16.552295091187520626
16.000000000000000000  1.000000000000000000  16.967593246867597401
17.000000000000000000  1.062500000000000000  17.380295789074676094
18.000000000000000000  1.125000000000000000  17.790418940392960678
19.000000000000000000  1.187500000000000000  18.197978822015503653
20.000000000000000000  1.250000000000000000  18.602991454377907132
21.000000000000000000  1.312500000000000000  19.005472757788044902
22.000000000000000000  1.375000000000000000  19.405438553051869377
23.000000000000000000  1.437500000000000000  19.802904562095296370
24.000000000000000000  1.500000000000000000  20.197886408582199635
25.000000000000000000  1.562500000000000000  20.590399618528561376
26.000000000000000000  1.625000000000000000  20.980459620912757401
27.000000000000000000  1.687500000000000000  21.368081748282051535
28.000000000000000000  1.750000000000000000  21.753281237355288624
29.000000000000000000  1.812500000000000000  22.136073229621818115
30.000000000000000000  1.875000000000000000  22.516472771936680175
31.000000000000000000  1.937500000000000000  22.894494817112075680
32.000000000000000000  2.000000000000000000  23.270154224505123608
33.000000000000000000  2.062500000000000000  23.643465760601966252
34.000000000000000000  2.125000000000000000  24.014444099598204474
35.000000000000000000  2.187500000000000000  24.383103823975716296
36.000000000000000000  2.250000000000000000  24.749459425075869490
37.000000000000000000  2.312500000000000000  25.113525303669145927
38.000000000000000000  2.375000000000000000  25.475315770521213210
39.000000000000000000  2.437500000000000000  25.834845046955454251
40.000000000000000000  2.500000000000000000  26.192127265411983217
41.000000000000000000  2.562500000000000000  26.547176470003158499
42.000000000000000000  2.625000000000000000  26.900006617065638892
43.000000000000000000  2.687500000000000000  27.250631575708979426
44.000000000000000000  2.750000000000000000  27.599065128360798838
45.000000000000000000  2.812500000000000000  27.945320971308543534
46.000000000000000000  2.875000000000000000  28.289412715237865825
47.000000000000000000  2.937500000000000000  28.631353885767630629
48.000000000000000000  3.000000000000000000  28.971157923981582627
49.000000000000000000  3.062500000000000000  29.308838186956698735
50.000000000000000000  3.125000000000000000  29.644407948288218790
51.000000000000000000  3.187500000000000000  29.977880398611418400
52.000000000000000000  3.250000000000000000  30.309268646120095525
53.000000000000000000  3.312500000000000000  30.638585717081845416
54.000000000000000000  3.375000000000000000  30.965844556350084815
55.000000000000000000  3.437500000000000000  31.291058027872896474
56.000000000000000000  3.500000000000000000  31.614238915198690449
57.000000000000000000  3.562500000000000000  31.935399921978699922
58.000000000000000000  3.625000000000000000  32.254553672466336423
59.000000000000000000  3.687500000000000000  32.571712712013422220
60.000000000000000000  3.750000000000000000  32.886889507563338952
61.000000000000000000  3.812500000000000000  33.200096448141067640
62.000000000000000000  3.875000000000000000  33.511345845340187566
63.000000000000000000  3.937500000000000000  33.820649933806812726
64.000000000000000000  4.000000000000000000  34.128020871720522678
65.000000000000000000  4.062500000000000000  34.433470741272266480
66.000000000000000000  4.125000000000000000  34.737011549139317879
67.000000000000000000  4.187500000000000000  35.038655226957196476
68.000000000000000000  4.250000000000000000  35.338413631788711200
69.000000000000000000  4.312500000000000000  35.636298546590033709
70.000000000000000000  4.375000000000000000  35.932321680673844355
71.000000000000000000  4.437500000000000000  36.226494670169635981
72.000000000000000000  4.500000000000000000  36.518829078481076067
73.000000000000000000  4.562500000000000000  36.809336396740569342
74.000000000000000000  4.625000000000000000  37.098028044260942693
75.000000000000000000  4.687500000000000000  37.384915368984309225
76.000000000000000000  4.750000000000000000  37.670009647928154095
77.000000000000000000  4.812500000000000000  37.953322087628606596
78.000000000000000000  4.875000000000000000  38.234863824580926916
79.000000000000000000  4.937500000000000000  38.514645925677292837
80.000000000000000000  5.000000000000000000  38.792679388641808202
81.000000000000000000  5.062500000000000000  39.068975142462797123
82.000000000000000000  5.125000000000000000  39.343544047822405219
83.000000000000000000  5.187500000000000000  39.616396897523515008
84.000000000000000000  5.250000000000000000  39.887544416913989664
85.000000000000000000  5.312500000000000000  40.156997264308280648
86.000000000000000000  5.375000000000000000  40.424766031406356603
87.000000000000000000  5.437500000000000000  40.690861243710067185
88.000000000000000000  5.500000000000000000  40.955293360936877889
89.000000000000000000  5.562500000000000000  41.218072777431025600
90.000000000000000000  5.625000000000000000  41.479209822572080668
91.000000000000000000  5.687500000000000000  41.738714761181007873
92.000000000000000000  5.750000000000000000  41.996597793923626796
93.000000000000000000  5.812500000000000000  42.252869057711606615
94.000000000000000000  5.875000000000000000  42.507538626100910051
95.000000000000000000  5.937500000000000000  42.760616509687778830
96.000000000000000000  6.000000000000000000  43.012112656502232255
97.000000000000000000  6.062500000000000000  43.262036952399093082
98.000000000000000000  6.125000000000000000  43.510399221446597551
99.000000000000000000  6.187500000000000000  43.757209226312554051
100.000000000000000000  6.250000000000000000  44.002476668648100144
101.000000000000000000  6.312500000000000000  44.246211189469050851
102.000000000000000000  6.375000000000000000  44.488422369534866618
103.000000000000000000  6.437500000000000000  44.729119729725276500
104.000000000000000000  6.500000000000000000  44.968312731414492589
105.000000000000000000  6.562500000000000000  45.206010776843150722
106.000000000000000000  6.625000000000000000  45.442223209487877966
107.000000000000000000  6.687500000000000000  45.676959314428579262
108.000000000000000000  6.750000000000000000  45.910228318713400597
109.000000000000000000  6.812500000000000000  46.142039391721439756
110.000000000000000000  6.875000000000000000  46.372401645523183333
111.000000000000000000  6.937500000000000000  46.601324135238662905
112.000000000000000000  7.000000000000000000  46.828815859393422727
113.000000000000000000  7.062500000000000000  47.054885760272213702
114.000000000000000000  7.125000000000000000  47.279542724270513077
115.000000000000000000  7.187500000000000000  47.502795582243820149
116.000000000000000000  7.250000000000000000  47.724653109854799027
117.000000000000000000  7.312500000000000000  47.945124027918204490
118.000000000000000000  7.375000000000000000  48.164217002743718865
119.000000000000000000  7.437500000000000000  48.381940646476571999
120.000000000000000000  7.500000000000000000  48.598303517436093557
121.000000000000000000  7.562500000000000000  48.813314120452119482
122.000000000000000000  7.625000000000000000  49.026980907199295245
123.000000000000000000  7.687500000000000000  49.239312276529297208
124.000000000000000000  7.750000000000000000  49.450316574800986302
125.000000000000000000  7.812500000000000000  49.660002096208479827
126.000000000000000000  7.875000000000000000  49.868377083107176873
127.000000000000000000  7.937500000000000000  50.075449726337758705
128.000000000000000000  8.000000000000000000  50.281228165548149889
129.000000000000000000  8.062500000000000000  50.485720489513475684
130.000000000000000000  8.125000000000000000  50.688934736454015706
131.000000000000000000  8.187500000000000000  50.890878894351175177
132.000000000000000000  8.250000000000000000  51.091560901261480865
133.000000000000000000  8.312500000000000000  51.290988645628594611
134.000000000000000000  8.375000000000000000  51.489169966593415495
135.000000000000000000  8.437500000000000000  51.686112654302206693
136.000000000000000000  8.500000000000000000  51.881824450212818078
137.000000000000000000  8.562500000000000000  52.076313047398990363
138.000000000000000000  8.625000000000000000  52.269586090852747873
139.000000000000000000  8.687500000000000000  52.461651177784915490
140.000000000000000000  8.750000000000000000  52.652515857923759768
141.000000000000000000  8.812500000000000000  52.842187633811732894
142.000000000000000000  8.875000000000000000  53.030673961100411873
143.000000000000000000  8.937500000000000000  53.217982248843533455
144.000000000000000000  9.000000000000000000  53.404119859788259816
145.000000000000000000  9.062500000000000000  53.589094110664582615
146.000000000000000000  9.125000000000000000  53.772912272472929374
147.000000000000000000  9.187500000000000000  53.955581570769972188
148.000000000000000000  9.250000000000000000  54.137109185952660084
149.000000000000000000  9.312500000000000000  54.317502253540453694
150.000000000000000000  9.375000000000000000  54.496767864455826214
151.000000000000000000  9.437500000000000000  54.674913065302973791
152.000000000000000000  9.500000000000000000  54.851944858644827718
153.000000000000000000  9.562500000000000000  55.027870203278297367
154.000000000000000000  9.625000000000000000  55.202696014507807831
155.000000000000000000  9.687500000000000000  55.376429164417132256
156.000000000000000000  9.750000000000000000  55.549076482139525979
157.000000000000000000  9.812500000000000000  55.720644754126155362
158.000000000000000000  9.875000000000000000  55.891140724412863960
159.000000000000000000  9.937500000000000000  56.060571094885283117
160.000000000000000000  10.000000000000000000  56.228942525542251474
161.000000000000000000  10.062500000000000000  56.396261634757614445
162.000000000000000000  10.125000000000000000  56.562534999540382330
163.000000000000000000  10.187500000000000000  56.727769155793254185
164.000000000000000000  10.250000000000000000  56.891970598569542972
165.000000000000000000  10.312500000000000000  57.055145782328480664
166.000000000000000000  10.375000000000000000  57.217301121188924640
167.000000000000000000  10.437500000000000000  57.378442989181493772
168.000000000000000000  10.500000000000000000  57.538577720499112900
169.000000000000000000  10.562500000000000000  57.697711609745994110
170.000000000000000000  10.625000000000000000  57.855850912185083246
171.000000000000000000  10.687500000000000000  58.013001843983929007
172.000000000000000000  10.750000000000000000  58.169170582459031493
173.000000000000000000  10.812500000000000000  58.324363266318663079
174.000000000000000000  10.875000000000000000  58.478585995904168726
175.000000000000000000  10.937500000000000000  58.631844833429767050
176.000000000000000000  11.000000000000000000  58.784145803220830828
177.000000000000000000  11.062500000000000000  58.935494891950703789
178.000000000000000000  11.125000000000000000  59.085898048876011046
179.000000000000000000  11.187500000000000000  59.235361186070534245
180.000000000000000000  11.250000000000000000  59.383890178657594561
181.000000000000000000  11.312500000000000000  59.531490865040986193
182.000000000000000000  11.375000000000000000  59.678169047134481673
183.000000000000000000  11.437500000000000000  59.823930490589887654
184.000000000000000000  11.500000000000000000  59.968780925023700945
185.000000000000000000  11.562500000000000000  60.112726044242300816
186.000000000000000000  11.625000000000000000  60.255771506465784171
187.000000000000000000  11.687500000000000000  60.397922934550372531
188.000000000000000000  11.750000000000000000  60.539185916209433458
189.000000000000000000  11.812500000000000000  60.679566004233123522
190.000000000000000000  11.875000000000000000  60.819068716706667033
191.000000000000000000  11.937500000000000000  60.957699537227249209
192.000000000000000000  12.000000000000000000  61.095463915119580633
193.000000000000000000  12.062500000000000000  61.232367265650083255
194.000000000000000000  12.125000000000000000  61.368414970239768991
195.000000000000000000  12.187500000000000000  61.503612376675768303
196.000000000000000000  12.250000000000000000  61.637964799321544263
197.000000000000000000  12.312500000000000000  61.771477519325785011
198.000000000000000000  12.375000000000000000  61.904155784829995923
199.000000000000000000  12.437500000000000000  62.036004811174805695
200.000000000000000000  12.500000000000000000  62.167029781104965025
201.000000000000000000  12.562500000000000000  62.297235844973059216
202.000000000000000000  12.625000000000000000  62.426628120941977329
203.000000000000000000  12.687500000000000000  62.555211695186088150
204.000000000000000000  12.750000000000000000  62.682991622091172701
205.000000000000000000  12.812500000000000000  62.809972924453106202
206.000000000000000000  12.875000000000000000  62.936160593675275265
207.000000000000000000  12.937500000000000000  63.061559589964801376
208.000000000000000000  13.000000000000000000  63.186174842527520923
209.000000000000000000  13.062500000000000000  63.310011249761721785
210.000000000000000000  13.125000000000000000  63.433073679450707516
211.000000000000000000  13.187500000000000000  63.555366968954139395
212.000000000000000000  13.250000000000000000  63.676895925398177667
213.000000000000000000  13.312500000000000000  63.797665325864436170
214.000000000000000000  13.375000000000000000  63.917679917577785886
215.000000000000000000  13.437500000000000000  64.036944418092929254
216.000000000000000000  13.500000000000000000  64.155463515479851822
217.000000000000000000  13.562500000000000000  64.273241868508108610
218.000000000000000000  13.625000000000000000  64.390284106829938082
219.000000000000000000  13.687500000000000000  64.506594831162246351
220.000000000000000000  13.750000000000000000  64.622178613467482933
221.000000000000000000  13.812500000000000000  64.737039997133308589
222.000000000000000000  13.875000000000000000  64.851183497151225765
223.000000000000000000  13.937500000000000000  64.964613600294029538
224.000000000000000000  14.000000000000000000  65.077334765292192742
225.000000000000000000  14.062500000000000000  65.189351423009114228
226.000000000000000000  14.125000000000000000  65.300667976615301313
227.000000000000000000  14.187500000000000000  65.411288801761457989
228.000000000000000000  14.250000000000000000  65.521218246750450476
229.000000000000000000  14.312500000000000000  65.630460632708263802
230.000000000000000000  14.375000000000000000  65.739020253753835732
231.000000000000000000  14.437500000000000000  65.846901377167867508
232.000000000000000000  14.500000000000000000  65.954108243560568781
233.000000000000000000  14.562500000000000000  66.060645067038308298
234.000000000000000000  14.625000000000000000  66.166516035369312476
235.000000000000000000  14.687500000000000000  66.271725310148255517
236.000000000000000000  14.750000000000000000  66.376277026959826344
237.000000000000000000  14.812500000000000000  66.480175295541329206
238.000000000000000000  14.875000000000000000  66.583424199944190036
239.000000000000000000  14.937500000000000000  66.686027798694539115
240.000000000000000000  15.000000000000000000  66.787990124952699489
241.000000000000000000  15.062500000000000000  66.889315186671751690
242.000000000000000000  15.125000000000000000  66.990006966755046847
243.000000000000000000  15.187500000000000000  67.090069423212824518
244.000000000000000000  15.250000000000000000  67.189506489317750493
245.000000000000000000  15.312500000000000000  67.288322073759516684
246.000000000000000000  15.375000000000000000  67.386520060798517306
247.000000000000000000  15.437500000000000000  67.484104310418530304
248.000000000000000000  15.500000000000000000  67.581078658478418220
249.000000000000000000  15.562500000000000000  67.677446916862933790
250.000000000000000000  15.625000000000000000  67.773212873632544984
251.000000000000000000  15.687500000000000000  67.868380293172336337
252.000000000000000000  15.750000000000000000  67.962952916340015008
253.000000000000000000  15.812500000000000000  68.056934460612893645
254.000000000000000000  15.875000000000000000  68.150328620234063237
255.000000000000000000  15.937500000000000000  68.243139066357599631
256.000000000000000000  16.000000000000000000  68.335369447192860548
257.000000000000000000  16.062500000000000000  68.427023388147901528
258.000000000000000000  16.125000000000000000  68.518104491971982384
259.000000000000000000  16.187500000000000000  68.608616338897164155
260.000000000000000000  16.250000000000000000  68.698562486779053415
261.000000000000000000  16.312500000000000000  68.787946471236679713
262.000000000000000000  16.375000000000000000  68.876771805791449310
263.000000000000000000  16.437500000000000000  68.965041982005246268
264.000000000000000000  16.500000000000000000  69.052760469617709305
265.000000000000000000  16.562500000000000000  69.139930716682599154
266.000000000000000000  16.625000000000000000  69.226556149703327492
267.000000000000000000  16.687500000000000000  69.312640173767675833
268.000000000000000000  16.750000000000000000  69.398186172681633366
269.000000000000000000  16.812500000000000000  69.483197509102367917
270.000000000000000000  16.875000000000000000  69.567677524670472167
271.000000000000000000  16.937500000000000000  69.651629540141286157
272.000000000000000000  17.000000000000000000  69.735056855515409779
273.000000000000000000  17.062500000000000000  69.817962750168433672
274.000000000000000000  17.125000000000000000  69.900350482979874300
275.000000000000000000  17.187500000000000000  69.982223292461256392
276.000000000000000000  17.250000000000000000  70.063584396883371141
277.000000000000000000  17.312500000000000000  70.144436994402852292
278.000000000000000000  17.375000000000000000  70.224784263187828515
279.000000000000000000  17.437500000000000000  70.304629361542907873
280.000000000000000000  17.500000000000000000  70.383975428033267008
281.000000000000000000  17.562500000000000000  70.462825581608058201
282.000000000000000000  17.625000000000000000  70.541182921723006416
283.000000000000000000  17.687500000000000000  70.619050528462238958
284.000000000000000000  17.750000000000000000  70.696431462659347744
285.000000000000000000  17.812500000000000000  70.773328766017726821
286.000000000000000000  17.875000000000000000  70.849745461230114074
287.000000000000000000  17.937500000000000000  70.925684552097422397
288.000000000000000000  18.000000000000000000  71.001149023646817682
289.000000000000000000  18.062500000000000000  71.076141842249029423
290.000000000000000000  18.125000000000000000  71.150665955734979207
291.000000000000000000  18.187500000000000000  71.224724293511641804
292.000000000000000000  18.250000000000000000  71.298319766677195730
293.000000000000000000  18.312500000000000000  71.371455268135463257
294.000000000000000000  18.375000000000000000  71.444133672709611460
295.000000000000000000  18.437500000000000000  71.516357837255171148
296.000000000000000000  18.500000000000000000  71.588130600772331036
297.000000000000000000  18.562500000000000000  71.659454784517507164
298.000000000000000000  18.625000000000000000  71.730333192114272833
299.000000000000000000  18.687500000000000000  71.800768609663563780
300.000000000000000000  18.750000000000000000  71.870763805853172812
301.000000000000000000  18.812500000000000000  71.940321532066590748
302.000000000000000000  18.875000000000000000  72.009444522491179441
303.000000000000000000  18.937500000000000000  72.078135494225605839
304.000000000000000000  19.000000000000000000  72.146397147386693405
305.000000000000000000  19.062500000000000000  72.214232165215520354
306.000000000000000000  19.125000000000000000  72.281643214182921042
307.000000000000000000  19.187500000000000000  72.348632944094276809
308.000000000000000000  19.250000000000000000  72.415203988193681539
309.000000000000000000  19.312500000000000000  72.481358963267467743
310.000000000000000000  19.375000000000000000  72.547100469747050511
311.000000000000000000  19.437500000000000000  72.612431091811131978
312.000000000000000000  19.500000000000000000  72.677353397487308939
313.000000000000000000  19.562500000000000000  72.741869938753012548
314.000000000000000000  19.625000000000000000  72.805983251635808529
315.000000000000000000  19.687500000000000000  72.869695856313086324
316.000000000000000000  19.750000000000000000  72.933010257211122962
317.000000000000000000  19.812500000000000000  72.995928943103550068
318.000000000000000000  19.875000000000000000  73.058454387209152969
319.000000000000000000  19.937500000000000000  73.120589047289101359
320.000000000000000000  20.000000000000000000  73.182335365743540478
321.000000000000000000  20.062500000000000000  73.243695769707642285
322.000000000000000000  20.125000000000000000  73.304672671146974494
323.000000000000000000  20.187500000000000000  73.365268466952301196
324.000000000000000000  20.250000000000000000  73.425485539033843452
325.000000000000000000  20.312500000000000000  73.485326254414886193
326.000000000000000000  20.375000000000000000  73.544792965324788270
327.000000000000000000  20.437500000000000000  73.603888009291509320
328.000000000000000000  20.500000000000000000  73.662613709233440318
329.000000000000000000  20.562500000000000000  73.720972373550736734
330.000000000000000000  20.625000000000000000  73.778966296216040632
331.000000000000000000  20.687500000000000000  73.836597756864691178
332.000000000000000000  20.750000000000000000  73.893869020884281440
333.000000000000000000  20.812500000000000000  73.950782339503760454
334.000000000000000000  20.875000000000000000  74.007339949881867369
335.000000000000000000  20.937500000000000000  74.063544075195110850
336.000000000000000000  21.000000000000000000  74.119396924725137410
337.000000000000000000  21.062500000000000000  74.174900693945602370
338.000000000000000000  21.125000000000000000  74.230057564608443954
339.000000000000000000  21.187500000000000000  74.284869704829645798
340.000000000000000000  21.250000000000000000  74.339339269174459446
341.000000000000000000  21.312500000000000000  74.393468398742115255
342.000000000000000000  21.375000000000000000  74.447259221249979078
343.000000000000000000  21.437500000000000000  74.500713851117168929
344.000000000000000000  21.500000000000000000  74.553834389547688488
345.000000000000000000  21.562500000000000000  74.606622924613020587
346.000000000000000000  21.625000000000000000  74.659081531334194892
347.000000000000000000  21.687500000000000000  74.711212271763358217
348.000000000000000000  21.750000000000000000  74.763017195064833231
349.000000000000000000  21.812500000000000000  74.814498337595679800
350.000000000000000000  21.875000000000000000  74.865657722985702094
351.000000000000000000  21.937500000000000000  74.916497362217043587
352.000000000000000000  22.000000000000000000  74.967019253703185200
353.000000000000000000  22.062500000000000000  75.017225383367545533
354.000000000000000000  22.125000000000000000  75.067117724721498462
355.000000000000000000  22.187500000000000000  75.116698238941992827
356.000000000000000000  22.250000000000000000  75.165968874948603684
357.000000000000000000  22.312500000000000000  75.214931569480171447
358.000000000000000000  22.375000000000000000  75.263588247170915224
359.000000000000000000  22.437500000000000000  75.311940820626091408
360.000000000000000000  22.500000000000000000  75.359991190497183311
361.000000000000000000  22.562500000000000000  75.407741245556579202
362.000000000000000000  22.625000000000000000  75.455192862771852447
363.000000000000000000  22.687500000000000000  75.502347907379530056
364.000000000000000000  22.750000000000000000  75.549208232958406484
365.000000000000000000  22.812500000000000000  75.595775681502416887
366.000000000000000000  22.875000000000000000  75.642052083493027226
367.000000000000000000  22.937500000000000000  75.688039257971198026
368.000000000000000000  23.000000000000000000  75.733739012608879193
369.000000000000000000  23.062500000000000000  75.779153143780078494
370.000000000000000000  23.125000000000000000  75.824283436631446875
371.000000000000000000  23.187500000000000000  75.869131665152494293
372.000000000000000000  23.250000000000000000  75.913699592245293957
373.000000000000000000  23.312500000000000000  75.957988969793760248
374.000000000000000000  23.375000000000000000  76.002001538732542940
375.000000000000000000  23.437500000000000000  76.045739029115466678
376.000000000000000000  23.500000000000000000  76.089203160183501495
377.000000000000000000  23.562500000000000000  76.132395640432349637
378.000000000000000000  23.625000000000000000  76.175318167679648695
379.000000000000000000  23.687500000000000000  76.217972429131648937
380.000000000000000000  23.750000000000000000  76.260360101449577996
381.000000000000000000  23.812500000000000000  76.302482850815522397
382.000000000000000000  23.875000000000000000  76.344342332997925382
383.000000000000000000  23.937500000000000000  76.385940193416686839
384.000000000000000000  24.000000000000000000  76.427278067207836898
385.000000000000000000  24.062500000000000000  76.468357579287783210
386.000000000000000000  24.125000000000000000  76.509180344417231368
387.000000000000000000  24.187500000000000000  76.549747967264622162
388.000000000000000000  24.250000000000000000  76.590062042469213566
389.000000000000000000  24.312500000000000000  76.630124154703779027
390.000000000000000000  24.375000000000000000  76.669935878736879431
391.000000000000000000  24.437500000000000000  76.709498779494779797
392.000000000000000000  24.500000000000000000  76.748814412122939643
393.000000000000000000  24.562500000000000000  76.787884322047176511
394.000000000000000000  24.625000000000000000  76.826710045034374730
395.000000000000000000  24.687500000000000000  76.865293107252909977
396.000000000000000000  24.750000000000000000  76.903635025332576447
397.000000000000000000  24.812500000000000000  76.941737306424244025
398.000000000000000000  24.875000000000000000  76.979601448259089125
399.000000000000000000  24.937500000000000000  77.017228939207470262
400.000000000000000000  25.000000000000000000  77.054621258337419931
401.000000000000000000  25.062500000000000000  77.091779875472809636
402.000000000000000000  25.125000000000000000  77.128706251251102799
403.000000000000000000  25.187500000000000000  77.165401837180780831
404.000000000000000000  25.250000000000000000  77.201868075698399707
405.000000000000000000  25.312500000000000000  77.238106400225291281
406.000000000000000000  25.375000000000000000  77.274118235223880902
407.000000000000000000  25.437500000000000000  77.309904996253735021
408.000000000000000000  25.500000000000000000  77.345468090027154062
409.000000000000000000  25.562500000000000000  77.380808914464481063
410.000000000000000000  25.625000000000000000  77.415928858749083474
411.000000000000000000  25.687500000000000000  77.450829303381908630
412.000000000000000000  25.750000000000000000  77.485511620235769215
413.000000000000000000  25.812500000000000000  77.519977172609301874
414.000000000000000000  25.875000000000000000  77.554227315280499511
415.000000000000000000  25.937500000000000000  77.588263394560001984
416.000000000000000000  26.000000000000000000  77.622086748344003126
417.000000000000000000  26.062500000000000000  77.655698706166859324
418.000000000000000000  26.125000000000000000  77.689100589253314411
419.000000000000000000  26.187500000000000000  77.722293710570482972
420.000000000000000000  26.250000000000000000  77.755279374879421539
421.000000000000000000  26.312500000000000000  77.788058878786429773
422.000000000000000000  26.375000000000000000  77.820633510794010590
423.000000000000000000  26.437500000000000000  77.853004551351546070
424.000000000000000000  26.500000000000000000  77.885173272905603881
425.000000000000000000  26.562500000000000000  77.917140939949945277
426.000000000000000000  26.625000000000000000  77.948908809075263093
427.000000000000000000  26.687500000000000000  77.980478129018536038
428.000000000000000000  26.750000000000000000  78.011850140712169832
429.000000000000000000  26.812500000000000000  78.043026077332712021
430.000000000000000000  26.875000000000000000  78.074007164349382037
431.000000000000000000  26.937500000000000000  78.104794619572203374
432.000000000000000000  27.000000000000000000  78.135389653199879945
433.000000000000000000  27.062500000000000000  78.165793467867374034
434.000000000000000000  27.125000000000000000  78.196007258693200015
435.000000000000000000  27.187500000000000000  78.226032213326362807
436.000000000000000000  27.250000000000000000  78.255869511993068954
437.000000000000000000  27.312500000000000000  78.285520327543110852
438.000000000000000000  27.375000000000000000  78.314985825495966765
439.000000000000000000  27.437500000000000000  78.344267164086616617
440.000000000000000000  27.500000000000000000  78.373365494311073576
441.000000000000000000  27.562500000000000000  78.402281959971631409
442.000000000000000000  27.625000000000000000  78.431017697721813420
443.000000000000000000  27.687500000000000000  78.459573837111051375
444.000000000000000000  27.750000000000000000  78.487951500629108637
445.000000000000000000  27.812500000000000000  78.516151803750176441
446.000000000000000000  27.875000000000000000  78.544175854976742812
447.000000000000000000  27.937500000000000000  78.572024755883134617
448.000000000000000000  28.000000000000000000  78.599699601158860673
449.000000000000000000  28.062500000000000000  78.627201478651613797
450.000000000000000000  28.125000000000000000  78.654531469410045474
451.000000000000000000  28.187500000000000000  78.681690647726227894
452.000000000000000000  28.250000000000000000  78.708680081177945453
453.000000000000000000  28.312500000000000000  78.735500830670588357
454.000000000000000000  28.375000000000000000  78.762153950478904108
455.000000000000000000  28.437500000000000000  78.788640488288407937
456.000000000000000000  28.500000000000000000  78.814961485236608496
457.000000000000000000  28.562500000000000000  78.841117975953878272
458.000000000000000000  28.625000000000000000  78.867110988604167687
459.000000000000000000  28.687500000000000000  78.892941544925392350
460.000000000000000000  28.750000000000000000  78.918610660269607138
461.000000000000000000  28.812500000000000000  78.944119343642924491
462.000000000000000000  28.875000000000000000  78.969468597745162697
463.000000000000000000  28.937500000000000000  78.994659419009252588
464.000000000000000000  29.000000000000000000  79.019692797640445292
465.000000000000000000  29.062500000000000000  79.044569717655193131
466.000000000000000000  29.125000000000000000  79.069291156919845776
467.000000000000000000  29.187500000000000000  79.093858087189090611
468.000000000000000000  29.250000000000000000  79.118271474144165722
469.000000000000000000  29.312500000000000000  79.142532277430760246
470.000000000000000000  29.375000000000000000  79.166641450696815241
471.000000000000000000  29.437500000000000000  79.190599941629955083
472.000000000000000000  29.500000000000000000  79.214408691994762535
473.000000000000000000  29.562500000000000000  79.238068637669798022
474.000000000000000000  29.625000000000000000  79.261580708684363117
475.000000000000000000  29.687500000000000000  79.284945829255079275
476.000000000000000000  29.750000000000000000  79.308164917822239204
477.000000000000000000  29.812500000000000000  79.331238887085845590
478.000000000000000000  29.875000000000000000  79.354168644041564562
479.000000000000000000  29.937500000000000000  79.376955090016309668
480.000000000000000000  30.000000000000000000  79.399599120703712174

Let's try plotting the results. We first need to import the required libraries and methods


In [3]:
%matplotlib inline
import numpy as np
from matplotlib import pyplot

Next, we create numpy arrays to store the (x,y) values


In [4]:
my_time = np.zeros(nsteps)
my_temp = np.zeros(nsteps)

We have to re write the loop to store the values in the arrays. Remember that numpy arrays start from 0.


In [5]:
T = T0
my_temp[0] = T0
for i in range(1,nsteps):
    T = T - r*(T-Ts)*dt
    my_time[i] = i*dt
    my_temp[i] = T

In [6]:
pyplot.plot(my_time, my_temp, color='#003366', ls='-', lw=3)
pyplot.xlabel('time')
pyplot.ylabel('temperature');


We could have saved effort by defining


In [7]:
my_time = np.linspace(0.,tmax,nsteps)

In [8]:
pyplot.plot(my_time, my_temp, color='#003366', ls='-', lw=3)
pyplot.xlabel('time')
pyplot.ylabel('temperature');


Alternatively, and in order to re use code in future problems, we could have created a function.


In [9]:
def euler(y, f, dx):
    """Computes y_new = y + f*dx
    
    Parameters
    ----------
    y  : float
        old value of y_n at x_n
    f  : float
        first derivative f(x,y) evaluated at (x_n,y_n)
    dx : float
        x step
    """
    
    return y + f*dx

In [10]:
T = T0
for i in range(1,nsteps):
    T = euler(T, -r*(T-Ts), dt)
    my_temp[i] = T

Actually, for this particularly simple case, calling a function may introduce unecessary overhead, but it is a an example that we will find useful for future applications. For a simple function like this we could have used a "lambda" function (more about lambda functions here).


In [11]:
euler = lambda y, f, dx: y + f*dx

Now, let's study the effects of different time steps on the convergence:


In [12]:
dt = 1.
#my_color = ['#003366','#663300','#660033','#330066']
my_color = ['red', 'green', 'blue', 'black']
for j in range(0,4):
    nsteps = int(tmax/dt)    #the arrays will have different size for different time steps
    my_time = np.linspace(dt,tmax,nsteps) 
    my_temp = np.zeros(nsteps)
    T = T0
    for i in range(1,nsteps):
        T = euler(T, -r*(T-Ts), dt)
        my_temp[i] = T
        
    pyplot.plot(my_time, my_temp, color=my_color[j], ls='-', lw=3)
    dt = dt/2.

pyplot.xlabel('time');
pyplot.ylabel('temperature');
pyplot.xlim(8,10)
pyplot.ylim(48,58);


Challenge 1.1

To properly study convergence, one possibility it so look at the result at a given time, for different time steps. Modify the previous program to print the temperature at $t=10$ as a function of $\Delta t$.


In [ ]: