Ma 449: Numerical Applied Mathematics Model Solutions to Homework Assignment 11 Prof. Wickerhauser Exercise 4 of Section 9.1, p.463. (a) If y(t) = C exp(-2t) + t exp(-2t), then y'(t)= -2 C exp(-2t) + exp(-2t) -2t exp(-2t), while f(t,y)=exp(-2t)-2y= exp(-2t)-2C exp(-2t)-2t exp(-2t). These are evidently equal, so y'=f. (b) f(t,y)=exp(-2t)-2y is differentiable in y at all (t,y), and f_y(t,y)=-2, so by theorem 9.1, L=2 is the smallest Lipschitz constant for f on any rectangle R. Exercise 6 of Section 9.1, p.463. The slope field consists of vectors tangent to the circles x^2+y^2=r^2 for various radii r, and the solutions are quarter-circles of radius sqrt(C). This may be plotted in Octave with the following commands: [t y] = meshgrid(0:1:4); h = 0.2; % small h keeps the flow field arrows small dt = h*ones(size(t)); % horizontal displacement of an arrow dy = h*(-t./y); % vertical displacement of an arrow quiver(t,y,dt,dy) % plot the arrows tt=0:0.1:4; % for plotting the solution curves yy1= sqrt(1-tt.*tt); % C=1 solution yy2= sqrt(2-tt.*tt); % C=2 solution yy4= sqrt(4-tt.*tt); % C=4 solution yy9= sqrt(9-tt.*tt); % C=9 solution hold on % plot the solutions on the flow field graph plot(tt, [yy1' yy2' yy4' yy9']) % transpose to get columns yy* See the file ./9-1-6.pdf for the resulting graph. Exercise 7 of Section 9.2, p.473. The Euler's method relation y_{k+1} = y_k + h f(t_k) for this IVP may be reapplied to get y_{k+1} = y_{k-1} + h f(t_{k-1}) + h f(t_k) y_{k+1} = y_{k-2} + h f(t_{k-2}) + h f(t_{k-1}) + h f(t_k) and so on. After k reapplications, the f terms may be grouped to get y_{k+1} = y_0 + h Sum_{j=0 to k} f(t_j) Substituting n for k+1, n-1 for k, and 0 for y_0=y(a), gives the desired formula: y_n = h Sum_{j=0 to n-1} f(t_j) Exercise 6 of Section 9.3, p.480. The Heun's method relation y_{k+1}=y_k +(h/2)[f(t_k)+f(t_{k+1))] for this IVP may be reapplied to get y_{k+1} = y_{k-1} + (h/2)[f(t_{k-1})+f(t_{k))] + (h/2)[f(t_k)+f(t_{k+1))] y_{k+1} = y_{k-2} + (h/2)[f(t_{k-2})+f(t_{k-1))] + (h/2)[f(t_{k-1})+f(t_{k))] + (h/2)[f(t_k)+f(t_{k+1))] and so on. After k reapplications, the f terms may be grouped to get y_{k+1} = y_0 + (h/2) Sum_{j=0 to k} [f(t_j)+f(t_{j+1})] Substituting n for k+1, n-1 for k, and 0 for y_0=y(a), gives the desired formula: y_n = (h/2) Sum_{j=0 to n-1} [f(t_j)+f(t_{j+1})] Problem H11.1. Write a computer program to solve the following initial value problem on the interval [0,1]: y'(t) = 2*t*y(t), for 0 #include #include /* Define the function used in the right hand side of the first-order differential equation of the initial value problem: */ double f ( double t, /* Independent variable */ double y) /* Dependent variable */ { return (2.0 * t * y); /* f(t,y) = 2 t y */ } /* Compute one new value of the approximate solution by Euler's method: */ double euler ( /* Next value by first-order Euler method */ double t, /* Current value of independent variable */ double y, /* Current value of dependent variable */ double h) /* Step size */ { return ( y + h*f(t,y) ); /* Euler increment */ } int main() { double h; /* Step size, or grid spacing */ double t; /* Independent variable */ double y; /* Dependent variable */ int k; /* Step counter */ int N; /* Total number of steps */ /* Get N */ printf("N-point Euler's method for y'(t)=2ty(t). Please enter N: "); scanf("%d", &N); assert(N>0); h = 1/(double)N; /* Compute */ y = 1; /* Initial y = 1 */ for ( k=0; k #include #include /* Define the function used in the right hand side of the first-order differential equation of the initial value problem: */ double f ( double t, /* Independent variable */ double y) /* Dependent variable */ { return (2.0 * t * y); /* f(t,y) = 2 t y */ } /* Compute one new value of the approximate solution by Euler's method: */ double euler ( /* Next value by first-order Euler method */ double t, /* Current value of independent variable */ double y, /* Current value of dependent variable */ double h) /* Step size */ { return ( y + h*f(t,y) ); /* Euler increment */ } /* Compute one new value of the approximate solution by Heun's method: */ double heun ( /* Next value by second-order Heun method */ double t, /* Current value of independent variable */ double y, /* Current value of dependent variable */ double h) /* Step size */ { return ( y + 0.5*h*( f(t,y) + f(t+h, euler(t,y,h)) ) ); } int main() { double h; /* Step size, or grid spacing */ double t; /* Independent variable */ double y; /* Dependent variable */ int k; /* Step counter */ int N; /* Total number of steps */ /* Get N */ printf("N-point Heun's method for y'(t)=2ty(t). Please enter N: "); scanf("%d", &N); assert(N>0); h = 1/(double)N; /* Compute */ y = 1; /* Initial y = 1 */ for ( k=0; k