Module 12: Differential Equations

1203 : Constructing Numerical Solutions

O B J E C T I V E

Explicit solutions to many differential equations are not possible.In this project we will learn how construct approximate solutions using special plots, animations, numeric tables, and then improve on the idea using the Runge-Kutta method.

S E T U P

In this project we will use the following command packages. Type and execute this line before begining the project below. If you re-enter the worksheet for this project, be sure to re-execute this statement before jumping to any point in the worksheet.

> restart; with(DEtools): with(plots):

Warning, the name changecoords has been redefined

___________________________________________________________________________________

A. The Idea

___________________________________________________________________________________

We will begin with a simple idea of using a starting point and slope to find other points and develop it further. This method leads to Eulers method of approximating the solution of differential equations.

THE BASIC CONCEPT

The basic idea comes from elementary algebra. Given a point (x1,y1) and a slope m, we can find another point on the same line. Solving m = Dy/Dx for Dy = mDx, we get the second point (x1 + Dx, y1 + m*Dx). In other words, x2 = x1 + Dx, and y2 = y1 + m*Dx. Here is an example.

> x1 := 5; y1:= 2; m := 5/3; delta := 2;
x2 := x1 + delta; y2 := y1 + m*delta;

> display( plot( [[x1,y1], [x2,y1]], x = 0..10, y=0..10, color = gold ),
plot( [[x2,y1], [x2, y2]], x = 0..10,color = green),
plot( [[x1,y1], [x2,y2]] , x = 0..10, thickness= 2, color = red),
plottools[disk]( [x1,y1], .3, color=blue),
plottools[disk]( [x2,y2], .3, color=yellow), scaling = constrained);

The blue dot is (x1,y1) and the yellow dot is (x2,y2). The gold horizontal segment represents Dx, the change in x, and the green vertical segment represents Dy, and the change in y.

We can continue this process and generate a line.

> delta := .5;

> plot({seq( [ [x1 + (k-1)*delta, y1+(k-1)*delta*m], [x1 + k*delta,y1 + k*delta*m]],
k = 1..12) }, x = -1..15, y = -2..15);

Each segment of different color is another extension using this method.

ANIMATED SOLUTIONS

We can animate this process. Here is a custom plotting command.

> y_construct := proc( df, a, b, num_iter, delta)
local n, a1, a2, x1,x2, y1, y2, across, up, diag, p, i;
n := num_iter * 2; a1 := a-1; a2 := round(a + n*delta + 1);
x2 := a; y2 := b; p[0] := pointplot( [x2,y2]):
for i from 1 to n by 2 do
x1 := x2; x2 := evalf(x1 + delta);
y1 := y2; y2 := evalf(y1 + df(x1,y1)*delta);
across := plot( [[x1,y1], [x2,y1]], x = a1..a2, color = gold ):
up := plot([[x2,y1], [x2,y2]], x = a1..a2, color = green):
diag:= plot([[x1,y1], [x2,y2]], x=a1..a2, thickness= 2, color = red):
p[i]:= plots[display]({p[i-1],across,up});
p[i+1] := plots[display]( {p[i] , diag }); od:
plots[display](seq( p[i], i = 0..n), insequence=true ); end:

We can use this function by defining the slope now as a function of x and y.

> m := (x,y) -> 4/3 ; y_construct( m, 5, 2, 6 , 0.4) ;

As we can see, the result is a line when the slope is a constant. The parameters are the slope (the derivative function), x1, y1, the number of iterations of the process, and the change in x.

The situation is far more interesting when the slope is not a constant, but a function that depends on x alone or x and y. In such cases, where we know the slope m = df(x,y), but not the function it comes from y = f(x), we are plotting a solution to a differential equation. Note the we need a starting point also.

> m := (x,y) -> .5*x^2 - x; y_construct( m, 2, 3, 6 , 0.4) ;

> m := (x,y) -> sin(x*y); y_construct( m, 2, 4, 10 , 0.5) ;

___________________________________________________________________________________

B. Accuracy of Solutions

___________________________________________________________________________________

The method above will yield graphs which are approximations to solutions. However, these solutions are only truly accurate at the initial point. The further away from the initial point,the less accurate the approximation.

> slope_plot := proc( df, a, b, delta)
local x1,x2,y1,y2,i,p1,p2,p3,n,a1,a2;
x2 := a; y2 := b; n := 12;
a1 := x2-1;
a2 := round(x2 + n*delta + 1);
for i from 1 to n do
x1 := x2;
x2 := evalf(x1 + delta);
y1 := y2;
y2 := evalf(y1 + df(x1,y1)*delta);
p1[i] := plot( [[x1,y1], [x2,y1]],x = a1..a2, color = blue):
p2[i] := plot( [[x2,y1], [x2,y2]],x = a1..a2, color = green):
p3[i] := plot( [[x1,y1], [x2,y2]],x = a1..a2,thickness= 2, color = gold):
od:
plots[display](seq( p1[i], i = 1..n),
seq( p2[i], i = 1..n),
seq( p3[i], i = 1..n) );
end:

Lets put this command to work.

> m := x -> x^2 - 2*x;
F := x -> (1/3)*x^3 - x^2 + 8/3;

> slope_plot ( m ,1, 2 , 0.25);

If we plot a solution along with the true solution, we can see what is happening.

The red graph is the true solution, and the gold graph is the approximate solution for delta = .25. Note that the two get further an further apart for x values to the left of 1.

> m := x -> x^2 - 2*x;
F := x -> (1/3)*x^3 - x^2 + 8/3;

> display( plot( F(x), x = 0..4, y = 0..10,
color = red), slope_plot ( m ,1,
2 , 0.25), plots[pointplot]
([1,2], color = blue) );

Here is another example.

> F := x -> .2*x^4 - 2*x^3 + x + 7; a := 2;

> display( slope_plot ( D(F) ,a, F(a) , 0.7 ), plot( F(x), x = -2..10, thickness = 2, color = red));

___________________________________________________________________________________

C. Numerical Solutions

___________________________________________________________________________________

Using these same ideas, we can construct numerical solution to differential equations.

COMPARING TRUE SOLUTION TO APPROXIMATION

> m := (x,y) -> x^2 + x + 1;

> F := x-> (1/3)*x^3 +(1/2)*x^2 + x - 85 ;

> a := 3; b := 2; delta := .5;

> y_approx := proc( df, a, b, delta, n)
local i,x,y1,y2; x := a; y2 := b;
for i from 1 to n do y1 := y2;
y2 := y2 + df(x,y1)*delta; x := a + i*delta; od;
evalf(y2); end:

Here is the numerical approximation to the solution for initial point (6,11) with delta = .5 and a list showing what happens as we take additional steps repeating the process.

The columns are Step, x, approximation to y, true y, error between true and approximate values.

> y_approx( m, 6, 11 , .5, 1);

> array( [ seq( [ k, a + delta*k ,
y_approx( m,a,b,delta,k), F(a + delta*k),
abs( y_approx( m, a,b,delta,k) - F(a +
delta*k))], k = 0..10)]);

APPROXIMATING A SOLUTION WHEN NONE EXISTS

We can use the same technique when we don't the true solution. In fact, even if there is no way to find to the true the solution, we can still approximate it.

The columns are Step, x, and the approximation to y.

> m := (x,y) -> sin(x)/x;

> a := Pi/2; b := 1; delta := .2;

> array( [ seq( [ k, a + delta*k , y_approx( m, a,b,delta,k) ], k = 0..10)]);

>