Initial Value Problems

1. Initial Value Problems Yi Lin Winter, 2007 Initial Value Problems

2. Ordinary Differential Equations In many areas of application we can measure how things change in some process From these measurements we would like to find the function that describes the process Examples: • Change in concentrations during chemical reaction • Heating or cooling of objects • Current flow in electrical circuits • Population dynamics Initial Value Problems

3. ODE’s • An ordinary differential equation takes the form G(t, x(t), x’(t), x”(t), ...) = 0 for all t, • Example • Here are some examples of ordinary differential equations: • x’(t) - 1 = 0 for all t (first-order) • x’(t) + x(t) - 1 = 0 for all t (first-order) • x”(t) - 1 = 0 for all t (second-order) • x”(t) - 2 t x(t) = 0 for all t (second-order) Initial Value Problems

4. Solution to ODE • A solution of an ordinary differential equation is a functionx that satisfies the equation for all values of t. • Example: • A solution of the equation x’(t) - 1 = 0, for example, is the function x defined by x(t) = t for all t. • There are many solutions. In general: x(t) = t + c (c is a constant) Initial Value Problems

5. Initial value problem • If x’(t) - 1 = 0 is given an initial value x(0)=4, then the only solution is x(t) = t + 4 • The additional conditions on the values of the function are referred to as initial conditions (though these conditions may specify the value of x or the value of its derivatives at any value of t, not necessarily the "first" value). • A differential equation together with a set of initial conditions is called an initial value problem Initial Value Problems

6. Analytic Solutions Some ODE’s have analytic solutions dy/dx = x + y -1 Has the solution y(x) = e^x – x Others have no analytic solution. For example: dy/dx = x^2 + y^2 Initial Value Problems

7. Initial Value Problems We use numerical methods to approximate the solution The methods we consider are called initial value problems since we must know the value of the function, y0, at some initial point x0 Initial Value Problems

8. The Euler Method • We “grow” the function from the starting value one step at a time • Think of the dy/dx in terms of discrete steps, delta_y and delta_x and multiply by delta_x. • When we step the independent variable by one stepsize, delta_x, we get the change in delta_y and the new value of the function Initial Value Problems

9. Euler Method We want to find an approximate solution to: dy/dx = f(x,y) y(x0) = y0 Now f(x0,y0) is the slope of the function at (x0,y0) Approximate the function value at x0+h by y0 + h*f(x0,y0) Repeat this process so that x_(n+1) = x_n + h y_(n+1) = y_n + hf(x_n,y_n) Initial Value Problems

10. #include <stdio.h> #include <math.h> typedef double (*DfDD)(double, double); // This is dy/dx = f(x,y) double dydx(double x, double y){ return x + y -1; //- (y * y); } // This is the euler method for each step. void euler(double x, double *y, DfDD dydx, double h){ *y = *y + h * dydx(x, *y); } int main(){ DfDD derivative = dydx; double x=0, y_euler=1.0, y_analy; double h; int i=0, n=10; printf("Enter the step h:\n"); scanf("%lf", &step); printf("x \t\ty(Analy)\ty(Euler)\n"); for(i=0; i<n; i++){ x += step; y_analy = exp(x) - x; euler(x, &y_euler, derivative, h); printf("%lf\t%lf\t%lf\n", x, y_analy, y_euler); } } A Very Simple Example Initial Value Problems

11. #define • After the #include commands, we can also define names using #define • Syntax • No ; at end • No = between name and value • In the program, the name is replaced by the expression • This can cause problems if complex expressions are used Initial Value Problems

12. Writing to a File In this version we output the results to a file. As in Fortran, we must give the file a name and open it. We name it by declaring it to be of type FILE and use fopen to open it. FILE * myfile = fopen("test2.csv","w+"); Then we can write to the file using fprintf fprintf (myfile," %f, %f, %f\n", x, y,exp(x)-x); When we are finished we must close the file fclose(myfile); Initial Value Problems

13. Writing to a File When we open a file, the filename is assigned a nonzero value if the operation is successful and a value of zero if it fails. if (myfile) { . . . } else printf("Could not open the file"); Initial Value Problems

14. Computing the Next Value We can pass a function as an argument to a general function that computes the next iteration from the previous one. First give a name (DfDD) to the type of functions we want to use. typedef double (*DfDD) (double,double); Initial Value Problems

15. Computing the Next Value Then define a function that uses this argument to compute the next value. We want the values of y to change, so we pass pointer to it. void eurler(double x, double *y, DfDD f, double h) { *y = *y + h * f(x,*y); return; } Initial Value Problems

16. Computing the Next Value We can define a function that we want to pass to the Euler function. double func(double x, double y) { return x+y-1; } Functions are stored in the memory of the machine at a specific address. We can pass a pointer to the function. DfDD f = func; …… euler(x,&y, f, h); Or simply: euler(x, &y, func, h); Initial Value Problems

17. Part of the Output File Initial Value Problems

18. Excel Generated Graph Initial Value Problems

19. Runge-Kutta Method The Euler method is not very accurate since the error tends to keep growing. In the (fourth-order) Runge_Kutta method the derivative is evalutated four times • At the initial point • Twice at a trial midpoint • At a trial endpoint Initial Value Problems

20. Runge-Kutta Formula Use the following to compute the next step k_1 = h*f(x_n, y_n) k_2 = h*f(x_n+h/2, y_n+k_1/2) K_3 = h*f(x_n+h/2, y_n+k_2/2) k_4 = h*f(x_n+h,y_n+k_3) y_(n+1) = y_n + (k_1+2*k_2+2*k_3+k_4)/6 Initial Value Problems

21. Runge-Kutta program // This is the Rounge Kunta method for each step void rk4(double x, double *y, DfDD dydx, double h){ double k1, k2, k3, k4; double h_half = h/2.0; k1 = h * dydx(x, *y); k2 = h * dydx(x + h/2, *y + k1/2 ); k3 = h * dydx(x + h/2, *y + k2/2 ); k4 = h * dydx(x + h, *y + k3); *y = *y + (k1+2.0*k2+2.0*k3+k4) /6.0; } Initial Value Problems

22. Extend Euler to 2 variables • Problems: with more than 1 variable dy/dx = 2x dz/dx = -x3 + 3y + z • To solve it analytically (difficult): • y = f(x) = x2 • z = f(x) = x3 + ex • To solve it numerically? Initial Value Problems

23. ODE for 2 variables typedef double (*DfDDD)(double, double, double); double fy(double x, double y, double z) { return 2*x; } double fz(double x, double y, double z) { return –x*x*x+3*y+z; } Initial Value Problems

24. Euler for 2 variables void eurler(double x, double *y, double *z, DfDDD fy, DfDDD fz, double h) { double y_cpy = *y; *y = *y + h * fy(x,*y,*z); *z = *z + h * fz(x,y_cpy,*z); return; } Initial Value Problems

25. Runge-Kutta for 2 variables void rk4_2_var(double x, double *y, double *z, DfDDD dydx, DfDDD dzdx, double h){ double k1, k2, k3, k4; double j1, j2, j3, j4; k1 = h * dydx(x, *y, *z); j1 = h * dzdx(x, *y, *z); k2 = h * dydx(x + h/2, *y + k1/2, *z+j1/2 ); j2 = h * dzdx(x + h/2, *y + k1/2, *z+j1/2 ); k3 = h * dydx(x + h/2, *y + k2/2, *z+j2/2 ); j3 = h * dzdx(x + h/2, *y + k2/2, *z+j2/2 ); k4 = h * dydx(x + h, *y + k3, *z+j3); j4 = h * dzdx(x + h, *y + k3, *z+j3); *y = *y + (k1+2.0*k2+2.0*k3+k4) /6.0; *z = *z + (j1+2.0*j2+2.0*j3+j4) /6.0; } Initial Value Problems