### Computational Physics Differential Equations - Guy Tel-Zur

```Computational Physics
Differential Equations
Dr. Guy Tel-Zur
Autumn Colors, by Bobby Mikul, http://www.publicdomainpictures.net
Version 10-11-2010 18:30
Agenda
• MHJ Chapter 13 & Koonin Chapter 2
• How to solve ODE using Matlab
• Scilab
Topics
•
•
•
•
•
Defining the scope of the discussion
Simple methods
Multi-Step methods
Runge-Kutta
Case Studies - Pendulum
The scope of the discussion
For a higher order ODE  a set of coupled 1st order equations:
Simple methods
Euler method:
Integration using higher order accuracy:
Taylor series expansion:
Local error!
Better than Euler’s method but useful only when it is easy to
differentiate f(x,y)
An Example
Let’s solve:
Boundary
conditiion
FORTRAN code - I
Euler’s method
Demo: virtual box, folder: ~/fortran, program: chap2a.f
Compilation: [email protected]/* <![CDATA[ */!function(t,e,r,n,c,a,p){try{t=document.currentScript||function(){for(t=document.getElementsByTagName('script'),e=t.length;e--;)if(t[e].getAttribute('data-cfhash'))return t[e]}();if(t&&(c=t.previousSibling)){p=t.parentNode;if(a=c.getAttribute('data-cfemail')){for(e='',r='0x'+a.substr(0,2)|0,n=2;a.length-n;n+=2)e+='%'+('0'+('0x'+a.substr(n,2)^r).toString(16)).slice(-2);p.replaceChild(document.createTextNode(decodeURIComponent(e)),c)}p.removeChild(t)}}catch(u){}}()/* ]]> */:~/fortran\$ fort77 -o chap2a chap2a.f
FORTRAN code - II
Taylor’s series method
FUNC(X,Y)=-X*Y
C-------scientific computing course, lecture 05 - differential
equations
C
Guy Tel-Zur, April 2011
C
example chap2a from Koonin
C
compile under ubuntu using: fort77 -o chap2a_taylor
chap2a_taylor.f
20
PRINT *,'ENTER STEP SIZE'
IF (H.LE.0) STOP
NSTEP=3./h
Y=1.
DO 10 IX=0,NSTEP-1
X=IX*H
Y=Y+H*FUNC(X,Y)+0.5*H*H*(-Y-FUNC(X,Y)*X)
DIFF=EXP(-0.5*(X+H)**2)-Y
PRINT *,IX,X+H,Y,DIFF
10
CONTINUE
GOTO 20
END
The Results
Euler’s
IX
0
1
2
3
4
5
X
0.5
1.0
1.5
2.0
2.5
3.0
Y
DIFF
1.
-.117503099
.75
-.143469334
.375
-.0503475331
.09375 .0415852815
0.
.0439369343
0.
.0111089963
Taylor’s
Y
DIFF
.875
.00749690272
.57421875
.0323119089
.287109375
.0375430919
.116638184
.0186970998
.0437393188 .000197614776
.0177690983 -.00666010194
Multi-Step methods
2 steps:
4 steps:
(So far) Explicit methods
Future = Function(Present && Past)
Implicit methods
Future = Function(Future && Present && Past)
Let’s calculate dy/dx at a mid way
between lattice points:
Let’s replace:
Rearrange:
This is a recursion relation!
A simplification occurs if f(x,y)=y*g(x), then the recursion equation becomes:
An example, suppose g(x)=-x
yn+1=(1-xnh/2)/(1+xn+1h/2)yn
This can be easily
calculated, for example:
Calculate y(x=1) for h=0.5
X0=0, y(0)=1
X1=0.5, y(0.5)=?
x2=1.0, y(1.0)=?
The solution:
0.5
1 − 0.5 2
1
1 =
( )
0.5
2
1+1 2
1
1 − 0 ∙ 0.5/2

=
∙ 1 = 0.8889
2
1 + 1 ∙ 0.5/2
1 = 0.62222
1  =  −1/2 = 0.60653
Error=-0.01569
Predictor-Corrector method
Predictor-Corrector method
Runge-Kutta
Proceed to:
Physics examples: Ideal harmonic oscillator – section 13.6.1
Physics Project – The pendulum, 13.7
I use a modified the C++ code from:
http://www.fys.uio.no/compphys/cp/programs/FYS3150/chapter13/cpp/program2.cpp
fout.close  fout.close()
Demo on folder:
\Lectures\05\CPP
ODEs in Matlab
function dydt = odefun(t,y)
a=0.001;
b=1.0;
dydt =b*t*sin(t)+a*t*t;
Usage:
[t1, y1]=ode23(@odefun,[0 100],0);
[t2, y2]=ode45(@odefun,[0 100],0);
plot(t1,y1,’r’);
hold on
plot(t2,y2,’b’);
hold off
Demo folder:
C:\Users\telzur\Documents\Weizmann\ScientificComputing\SC2011B\Lectures\05\Matlab
Output
http://www.scilab.org
Parallel tools for Multi-Core and
Distributed Parallel Computing
In preparation
Parallel execution
A new function (parallel_run) allows parallel computations and leverages multicore
architectures and their capacities.
In future Scilab versions:
http://help.scilab.org/docs/5.3.1/en_US/parallel_run.html
Xcos demo
```