1

 Solving Differential Equations using MATHCAD

2

 Previously worked through different numerical solutions by hand
 MathCAD has builtin solution methods
 RungeKutta 4^{th} Order
 Other solvers available

3

 Scale equations, parameters & initial conditions to remove units
 Manipulate equations to give vector of derivatives
 Give vector of initial conditions
 Call solver
 Plot results

4

 Solvers produce matrix of numbers
 1 Row per timestep
 1 Column per derivative +
1 Column for time

5

 Radioactive decay, Newton’s law of cooling etc
 A is amount of material, temperature difference, etc
 k is rate constant

6

 “The derivative of A with respect to t is –k times A”

7

 First order ODE
 1 member of D(t,A) vector
 1 member of ic vector

8

 TStart and TFinish times
 Number of points, N

9

 Use rkfixed() function
 Return matrix with 1 column per DE + 1 for independent variable
 Use column operator to strip off columns

10


11

 Found number greater than 10^307
 Use more points
 Reduce TFinish
 Can’t have anything with dimensions here
 Strip units from system before solution

12

 Double decay
 A => B => C
 K1 K2

13


14

 System of 2 DEs
 2 members in D(t,F) vector
 2 members in ic vector
 F_{0} refers to A => ic_{0}
 F_{1} refers to B => ic_{1}

15

 Solution parameters as before
 Call rkfixed() as before to create matrix

16


17

 Pendulum with drag
 Mass on spring
 Molecules on bonds
 LCR circuit
 Rewrite as system of first order equations
 Solve as before

18


19


20


21


22

 Have expressions for q1 & q2
 Now fill in D(t,q) vector
 Replace with vector subscripts
 q0 ® q_{0}
,q1 ® q_{1 },q2
® q_{2}

23

 Start with 1V across capacitor & no current flowing

24

 Examine over 0.1s
 Use 1000 points

25


26

 So far all systems have been ‘Relaxation to steady state’
 Can also model systems driven by ‘Forcing Function’

27

 Enables us to see resonance
 Put voltage source in loop

28

 Use symbolic solver to solve for q2 as before
 Compare with undriven case…

29

 Drive with sinusoidal waveform
 Define f(t) before D(t,q)

30

 No initial charge or current

31


32


33

 Plot Current (q1) vs Charge (q0)

34

 Timestep = (TFinishTStart)/N
 Small Timestep è
 More Accurate Solution
 Longer Computation Time
 Larger Solution Matrix
 How to choose the right timestep?

35

 How to choose the right timestep?
 Depends on time constants of system
 Experiment with changing N

36

 Rkadapt() function
 Adjusts timestep to get acceptable error
 Good for systems with wildly differing timescales

37


38

 Programming construct to simplify nested if construct

39


40


41

 Builtin powerful ODE solver
 Need to form derivative vector
 Experiment to find optimum timestep
 Use adaptive timestep with Rkadapt() function
