An Introduction to MathCAD You can always find a solution !!

roots When things go wrong x 5 79 x 4 2180 x 3 133764 x 2 d d x f x 5 x 4 316 x 3 6540 x 2 16039296 3098448 56 34 78 12 9 “Result is too large to Display” MathCAD #6  dpl 2001 2

Can we fix it ?

Yes we can !!

1.0 x 5 79 x 4 2180 x 3 133764 x 2 d d x f x 5.0 x 4 316 x 3 6540 x 2 3098448 16039296 roots 56.

34.

78.

12.

9.

turnings 27.033163438232028292

10.481803779012762137

31.678249774297679276

69.036717442947111153

Force solver to work with real numbers MathCAD #6  dpl 2001 3

What if its still broken ?

x 5 79 x 4 2180 x 3 133764 x 2 d d x f x 5 x 4 316 x 3 6540 x 2 16039296 3098448 roots 56 34 78 12 9 turnings 27.033163438446261373

10.481803778652327208

31.678249774081398831

69.036717443017189749

Risky !! – OK in this case MathCAD #6  dpl 2001 4

Still can’t solve it ?

 Use polyroots to find roots  Use coeffs keyword on solver to get coefficients of x cfsx 3098448 267528 6540 316 5 polyroots cfsx )  27.033

10.482

31.678

69.037

MathCAD #6  dpl 2001 5

Still can’t solve it ?

 Graph shows maxima & minima  => there must be solutions  Try guessing & using roots(f1(x),x) MathCAD #6  dpl 2001 6

Solving ODEs numerically  Produce numeric solution to system of ODEs.

 Must have initial conditions  Manipulate equations  Use one of several different solvers  Produces matrix of solutions MathCAD #6  dpl 2001 7

First order linear ODE #1  Radioactive decay, Newton’s law of cooling etc d d t A  A is amount of material temperature difference etc  k is rate constant MathCAD #6  dpl 2001 8

First order linear ODE #2  Define initial conditions as a vector  1 st order so only 1 element in vector  Can’t use units in ODE solver  Call vector ‘ic’  Element 0 = A at t=0 ic 0 100 MathCAD #6  dpl 2001 9

First order linear ODE #3  Now define ODE & manipulate for mathCAD d d t A k 0.1

) MathCAD #6  dpl 2001 10

First order linear ODE #4  Now define range of solution t0 t1 N 0 5 1000 Start time Finish Time Number of Points  And solve using rkfixed Soln (  t1  N  D )  Creates matrix ‘Soln’ containing solution MathCAD #6  dpl 2001 11

First order linear ODE #5 The Solution Matrix Soln  6 7 8 9 10 11 12 13 14 15 0 1 2 3 4 5 0 0 1 100 99.501

0.01 99.005

0.015 98.511

0.02

98.02

0.025 97.531

0.03 97.045

0.035 96.561

0.04 96.079

0.045

95.6

0.05 95.123

0.055 94.649

0.06 94.176

0.065 93.707

0.07 93.239

0.075 92.774

0.08 92.312

Column 0 holds t values Column 1 holds A values 1 row of matrix per timestep (0..N) MathCAD #6  dpl 2001 12

i First order linear ODE #6  Plot the results  Use M <> to extract columns  Use subscripting to extract rows 100 Soln 1 i 50 0 0 2 Soln 0 i 4 MathCAD #6  dpl 2001 6 13

Second Order ODE  Same steps as for first order  Slightly less obvious manipulation  Replace ODE by system of 1 st orders  Can use symbolic solver to formulate equations MathCAD #6  dpl 2001 14

Second Order ODE Example #1  Damped SHM – LCR circuit – Damped Pendulum L  d d t i q C 0 Substitute for i & divide through by L d d t q i Gives: d 2 d t 2 q R L  d d t q 1 LC  q 0 MathCAD #6  dpl 2001 15

Second Order ODE Example #2 d 2 d t 2 q R  L d d t q 1 LC  q 0 rewrite equation (by hand) substituting d 2 q q2 d t 2 d d t q q q0 q1 Gives manipulated equation: q2 R  q1 L 1  q0 0 MathCAD #6  dpl 2001 16

Second Order ODE Example #3 q2 R  q1 L 1  q0 0 q2  Use solver to solve for q2 R  q1 L 1  q0 ( R q1  C ( ) q0 )  Now ready to create D(t,q) function MathCAD #6  dpl 2001 17

Second Order ODE Example #4  Equations to create D(t,q) d d t q q1 d 2 q q2 (  C d t 2 (  2 row vector to hold D(t,q) ) q0 )  Change subscripts for suffixes ) q 1 R q 1  C ( ) q 0 d d t q d 2 d t 2 q MathCAD #6  dpl 2001 18

Second Order ODE Example #5  Specify initial conditions  For second order need 2 elements in ic vector ic 1 0 Charge (q) at t=0 Current (dq/dt) at t=0  Also need to specify constants in D(t,q) equations  Watch for lack of units !!

L 10 2 C 10 3 R 1 MathCAD #6  dpl 2001 19

Second Order ODE Example #6  Define times and number of points as before  Call rkfixed to solve t1 0.1

t0 0 N 1000 Soln  t1  N  D )  Matrix Soln filled with solution points MathCAD #6  dpl 2001 20

Second Order ODE Example #7 Soln  11 12 13 14 15 0 1 2 3 4 5 6 7 8 9 10 • The solution matrix 0 0 1 1 1 0.998

0.996

0.992

0.988

0.982

0.976

0.969

0.961

0.952

0.942 -102.078

0.932 -110.386

0.92 -118.501

0.908 -126.417

0.895 -134.129

0.881 -141.631

2 0 -9.949

-19.788

-29.51

-39.106

-48.568

-57.887

-67.055

-76.066

-84.912

-93.585

Column 0 holds t values Column 1 holds q values Column 2 holds dq/dt values MathCAD #6  dpl 2001 21

Second Order ODE Example #8  Graphing the solution Time Charge Soln   Soln   1 0.5

Charge 0 0.5

1 0 0.05

Time MathCAD #6  dpl 2001 0.1

22

Surely nothing can go wrong ?

 Solution relies on numeric integration which divides timestep up into smaller chunks for integration  If system is changing much faster than timestep, solution will fail  Clues to look for: – “found a number >10^307” – Singularity at t finish MathCAD #6  dpl 2001 23

Can we fix it ?

 Check time constants  Try more point in solution (N)  Try smaller interval  Use a different solver – ‘Stiff’ systems give problems – Special solvers for stiff systems – See Quicksheet & help system MathCAD #6  dpl 2001 24