Di erential equations using Scilab

Predator-Prey by Euler for different step sizes, Black = 0.25, Red = 0.1 Predator Prey As step size increases, the solution diverges more from the act...

0 downloads 73 Views 201KB Size
Differential equations using Scilab Kannan M. Moudgalya IIT Bombay [email protected]

Scilab Workshop 13 September 2010 Kannan Moudgalya

Differential equations using Scilab

1/31

Outline I

I

Weight reduction ODE model - analytical solution Numerical integration I I

I

Functions in Scilab Euler’s method

Predator-prey system I I I

Modelling Euler method - user created integrator Backward difference method - built-in function

Kannan Moudgalya

Differential equations using Scilab

2/31

Weight Reduction Model I I I I

Weight of person = x kg Tries to reduce weight Weight loss per month = 10% of weight Starting weight = 100 kg dx dt

= −0.1x

Initial conditions: x = 100

at t=0

Determine x(t) as a function of t. Kannan Moudgalya

Differential equations using Scilab

3/31

Analytical Solution of Simple Model I Recall the model: dx

= −0.1x dt x(t = 0) = 100 dx Cross multiplying, = −0.1dt x Integrating both sides from 0 to t, Z

dx

Z

= −0.1 dt x C + ln x(t) = −0.1t Kannan Moudgalya

Differential equations using Scilab

4/31

Analytical Solution of Simple Model II Using initial conditions, C = − ln 100 Thus, the final solution is, ln

x(t) 100

= −0.1t

or

Kannan Moudgalya

x(t) = 100e−0.1t

Differential equations using Scilab

5/31

Solution, Continued I I I I

Weight of person = x kg Tries to reduce weight Weight loss per month = 10% of weight Starting weight = 100 kg x(t) = 100e−0.1t

Compute, plot for 2 years, i.e. for 24 months: 1 2 3 4 5

T=0:0.1:24; p l o t 2 d (T, 1 0 0 ∗ e x p ( −0.1∗T ) ) ; x t i t l e ( ’ Weight v s . month ’ , . . . ’ Time i n months ’ , . . ’ Weight ( kg ) ’ ) Kannan Moudgalya

Differential equations using Scilab

6/31

Weight vs. month Weight (kg) 100 90 80 70 60 50 40 30 20 10

Time, in months

0 0

4

8

Kannan Moudgalya

12

16

20

Differential equations using Scilab

24

7/31

Need for Numerical Solution I I I

Exact solution is ok for simple models What if the model is complicated? Consider integrating the more difficult problem: dx dt

I

I

=

2 + 18t + 68t2 + 180t3 + 250t4 + 250t5 x2

with initial condition, x(t = 0) = 1 Analytical (i.e. exact) solution difficult to find Let us look for a numerical solution Kannan Moudgalya

Differential equations using Scilab

8/31

Simple numerical solution: Explicit Euler I

I

I

Suppose that we want to integrate dx = g(x, t) dt with initial condition: x(t = 0) = x0 Approximate numerical method - divide time into equal intervals: t0 , t1 , t2 , etc. xn − xn−1 = g(xn−1 , tn−1 ) ∆t Simplifying, xn − xn−1 = ∆t g(xn−1 , tn−1 ) xn = xn−1 + ∆t g(xn−1 , tn−1 ) Kannan Moudgalya

Differential equations using Scilab

9/31

Example revisited

Recall the problem statement for numerical solution: dx 2 + 18t + 68t2 + 180t3 + 250t4 + 250t5 = dt x2 with initial condition, x(t = 0) = 1 Recall the Euler method: dx = g(x, t) dt Solution for initial condition, x(t = 0) = x0 is, Kannan Moudgalya

Differential equations using Scilab

10/31

Scilab Code I

1 2 3 4 5 6 7 8

exec ( ” d i f f 1 . s c i ” ) ; exec ( ” Euler . s c i ” ) ; x0 =1; t 0 =0; T = 0 : 1 : 2 4 ; s o l = E u l e r ( x0 , t0 , T , d i f f 1 ) ; //

s o l

=

o d e ( x0 , t0 , T , d i f f 1 ) ;

p l o t 2 d (T , s o l ) , p a u s e p l o t 2 d (T,1+2∗T+5∗Tˆ 2 , 5 ) x t i t l e ( ’ x v s . t : P o l y n o m i a l Problem ’ , ’ t ’ , ’ x

Kannan Moudgalya

Differential equations using Scilab

11/31

Scilab Code II 1 2 3 4 5 6 7

1 2 3

f u n c t i o n x = E u l e r ( x0 , t0 , t , g ) n = l e n g t h ( t ) , x = x0 ; f o r j = 1 : n−1 x0 = x0 + ( t ( j +1)− t ( j ) ) ∗ g ( t ( j ) , x0 ) ; x = [ x x0 ] ; end ; endfunction

f u n c t i o n xdot = d i f f 1 ( t , x ) x d o t = (2+18∗ t +68∗ t ˆ2+180∗ t ˆ3+250∗ t ˆ4+250∗ t endfunction

Kannan Moudgalya

Differential equations using Scilab

12/31

Comparison of numerical, exact solution x vs. t: Polynomial Problem x 3200

2800

2400

2000

1600

1200

800

400 t 0 0

4

8

Kannan Moudgalya

12

16

20

Differential equations using Scilab

24

13/31

Predator-Prey Problem I I I

I

I

I

Population dynamics of predator-prey Prey can find food, but gets killed on meeting predator Examples: parasites and certain hosts; wolves and rabbits x1 (t) - number of prey; x2 (t) - number of predator at time t Prey, if left alone, grows at a rate proportional to x1

Kannan Moudgalya

Differential equations using Scilab

14/31

Predator-Prey Problem II I

I

I

I

Predator, on meeting prey, kills it ⇒ proportional to x1 x2 dx1 = 0.25x1 − 0.01x1 x2 dt Predator, if left alone, decrease by natural causes Predators increase their number on meeting prey dx2 = −x2 + 0.01x1 x2 dt Determine x1 (t), x2 (t) when x1 (0) = 80, x2 (0) = 30 Kannan Moudgalya

Differential equations using Scilab

15/31

Explicit Euler for a System of Equations I

dx1 dt dxN dt

= g1 (x1 , . . . , xn , t) .. . = gn (x1 , . . . , xn , t)

Kannan Moudgalya

Differential equations using Scilab

16/31

Explicit Euler for a System of Equations II     x1 g1 (x1 , . . . , xn , t − 1) d . ..  ..  =   . dt xn gn (x1 , . . . , xn , t − 1)

      x1 x1 g1 ((x1 , . . . , xn )|t−1 , t − 1) ..  ...  =  ...   + ∆t  . xn t xn t−1 gN ((x1 , . . . , xn )|t−1 , t − 1) Solution in vector form: xt = xt−1 + ∆tg(xt−1 ) Kannan Moudgalya

Differential equations using Scilab

17/31

Scilab Code for Predator-Prey Problem I

4

exec ( ” pred . s c i ” ) ; exec ( ” Euler . s c i ” ) ; x0 = [ 8 0 , 3 0 ] ’ ; t 0 =0; T = 0 : 0 . 1 : 2 0 ; T=T ’ ; s o l = E u l e r ( x0 , t0 , T , p r e d ) ;

5

//

1 2 3

6 7 8 9

s o l

=

o d e ( x0 , t0 , T , p r e d ) ;

clf (); p l o t 2 d (T , s o l ’ ) x s e t ( ’ window ’ , 1 ) plot2d ( sol (2 ,:) , sol (1 ,:))

Kannan Moudgalya

Differential equations using Scilab

18/31

Scilab Code for Predator-Prey Problem II 1 2 3 4 5 6 7

1 2 3 4

f u n c t i o n x = E u l e r ( x0 , t0 , t , g ) n = l e n g t h ( t ) , x = x0 ; f o r j = 1 : n−1 x0 = x0 + ( t ( j +1)− t ( j ) ) ∗ g ( t ( j ) , x0 ) ; x = [ x x0 ] ; end ; endfunction f u n c t i o n xdot = pred ( t , x ) xdot ( 1 ) = 0 . 2 5 ∗ x (1) −0.01∗ x ( 1 ) ∗ x ( 2 ) ; x d o t ( 2 ) = −x ( 2 ) + 0 . 0 1 ∗ x ( 1 ) ∗ x ( 2 ) ; endfunction

Kannan Moudgalya

Differential equations using Scilab

19/31

Predator-Prey Problem: Solution by Euler Predator−Prey by Euler for different step sizes, Black = 0.25, Red = 0.1 Prey 140

130

120

110

100

90

80

70 Predator 60 10

14

18

22

26

30

34

38

42

46

As step size increases, the solution diverges Kannan Moudgalya

Differential equations using Scilab

20/31

General method to handle stiff systems I

I

I

I

The predator-prey problem is an example of a stiff system Results because of sudden changes in the derivative Approximation of using previous time values does not work General approach to solve this problem: xn = xn−1 + ∆t g(xn , tn )

I

I

Requires solution by trial and error, as g(xn , tn ) is unknown Scilab has state of the art methods (ode) to solve such systems Kannan Moudgalya

Differential equations using Scilab

21/31

General method to handle stiff systems

I

Derived from ODEPACK I I I

FOSS In use for thirty years Bugs have been removed by millions of users

Kannan Moudgalya

Differential equations using Scilab

22/31

Predator-Prey Problem by Scilab Integrator I Execute the following code, after commenting out Euler and uncommenting ode: 1 2 3 4 5 6 7 8

exec ( ” pred . s c i ” ) ; exec ( ” Euler . s c i ” ) ; x0 = [ 8 0 , 3 0 ] ’ ; t 0 =0; T = 0 : 0 . 1 : 2 0 ; T=T ’ ; s o l = E u l e r ( x0 , t0 , T , p r e d ) ; / /

s o l

=

o d e ( x0

, t0

, T , p r e d

) ;

clf (); p l o t 2 d (T , s o l ’ ) x s e t ( ’ window ’ , 1 ) Kannan Moudgalya

Differential equations using Scilab

23/31

Predator-Prey Problem by Scilab Integrator II

9

1 2 3 4

plot2d ( sol (2 ,:) , sol (1 ,:)) f u n c t i o n xdot = pred ( t , x ) xdot ( 1 ) = 0 . 2 5 ∗ x (1) −0.01∗ x ( 1 ) ∗ x ( 2 ) ; x d o t ( 2 ) = −x ( 2 ) + 0 . 0 1 ∗ x ( 1 ) ∗ x ( 2 ) ; endfunction

Kannan Moudgalya

Differential equations using Scilab

24/31

Predator-Prey Problem: Scilab Integrator Predator−Prey by Scilab integrator Prey 130

120

110

100

90

80 Predator 70 15

19

23

27

31

35

39

Use the Scilab built-in integrator to get the correct solution Kannan Moudgalya

Differential equations using Scilab

25/31

Partial Differential Equations

Kannan Moudgalya

Differential equations using Scilab

26/31

Parabolic Differential Equations I

Heat conduction equation

I

Diffusion equation ∂u(t, x)

=c

∂t I

0≤x≤1

Boundary conditions: u(t, 0) = α,

I

∂x2

Initial condition: u(0, x) = g(x),

I

∂ 2 u(t, x)

u(t, 1) = β,

t≥0

Let um j be approximate solution at xj = j∆x, tm = m∆t um+1 − um j j

c

m m (um j−1 − 2uj + uj+1 ) 2 ∆tKannan Moudgalya (∆x)Differential equations using Scilab

=

27/31

Finite Difference Approach um+1 − um j j

=

c (∆x)2

m m (um j−1 − 2uj + uj+1 ),

∆t m m m um+1 = um j j + µ(uj−1 − 2uj + uj+1 )

µ=

c∆t (∆x)2

m m = µum j−1 + (1 − 2µ)uj + µuj+1

Write this equation at every spatial grid: m m = µum um+1 1 0 + (1 − 2µ)u1 + µu2 m m um+1 = µum 2 1 + (1 − 2µ)u2 + µu3 .. . m m um+1 = µum N N−1 + (1 − 2µ)uN + µuN+1 Kannan Moudgalya

Differential equations using Scilab

28/31

Finite Difference Approach - Continued m m um+1 = µum 1 0 + (1 − 2µ)u1 + µu2 m m um+1 = µum 2 1 + (1 − 2µ)u2 + µu3 .. . m m um+1 = µum N N−1 + (1 − 2µ)uN + µuN+1

In matrix form,

m+1    m  µ u1 1 − 2µ µ u1    u2   µ   1 − 2µ µ u  .    .2  +  = ..  ..      ..  . µu uN uN µ 1 − 2µ 

Kannan Moudgalya

Differential equations using Scilab

29/31

Conclusions

I

I I

I I

Scilab is ideal for educational institutions, including schools Built on a sound numerical platform It has good integrators for differential equations It is free Also suitable for industrial applications

Kannan Moudgalya

Differential equations using Scilab

30/31

Thank you

Kannan Moudgalya

Differential equations using Scilab

31/31