Author: Frank Plate Subj: solving differential equations Date: 1 April 1992 [Note: received by mail on a floppy. Thanx, Frank! -jkh-] This program solves every common differential equation and system (every potential and degree also with variable coefficients) on repetitive use of the Runge-Kutta-algorithm. * Attention : ÉÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ» º - avoid discontinuity on the sampled x-range ! º º - periodic functions require stepsize << one period! º ÈÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍͼ * It generates the following variables : - 'degree' - 'functions' ( list of a deq-system ) - 'initvalues' - 'xmax' 'xmin' - 'h' ( stepsize ) * y0 in this program equals y, y1 equals y', y2 equals y'' ect. * The progam is controlled by menu and will take the value of the last calculation when invalid or no input-data is entered at any step. * the result can be found in äDAT : x y y' y'' y''' ... 1. values [[ ] 2. values [ ] .... ... [ ]] * As result the second column of äDAT ( = y(x) ) is drawn. y', y'' etc. can be plotted by changing YCOL. ÉÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ» º Example #1 : º ÈÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍͼ ------------------------------------------------------------------------- deq: y'' = -y' - 1.25 * y ³ y (3/2*ã)=1 ³ y'(3/2*ã)=0 solution: y = e^(0.25*(3*ã-2*x)) * ( 0.5*cos(x) - sin(x) ) ------------------------------------------------------------------------- - press "DEQ" - enter 2 when asked for degree - enter highest derivative : :ëX(Y1)=Y2: '-Y1-1.25*Y0' in 'functions' stored : { :ëX(Y0): 'Y1' :ëX(Y1)=Y2: '-Y1-1.25*Y0' } - enter initial values : : x0 : '3/2*ã' :Y0(x0): 1 :Y1(x0): 0 stored in 'initvalues' : [ 4.71238898038 1 0 ] - enter x-range : :xmin: 0 :xmax: '2*ã' - enter stepsize : 'ã/4' ÉÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ» º Example #2 : º ÈÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍͼ -------------------------------------------------------------------------- deq: y''= ( 2 / y ) * (y')ý - ( 2 / x ) * y' ³ y (1) = 1 ³ y'(1) = 0.5 solution : y = 2 * x / ( 1 + x ) -------------------------------------------------------------------------- - press "DEQ" - degree : 2 - enter the highest derivative : :ëX(Y1)=Y2: '2/Y0*Y1^2-2/X*Y1' stored in 'functions' : { :ëX(Y0): 'Y1' :ëX(Y1)=Y2: '2/Y0*Y1^2-2/X*Y1' } - enter initial values : : x0 : 1 :Y0(x0): 1 :Y1(x0): 0.5 stored in 'initvalues' : [ 1 1 0.5 ] - enter x-range : :xmin: 0 :xmax: 5 - enter stepsize : 1.1 see what's happening when you change the x-range : - press "DEQ" - press "ENTER" three times ( 'degree' 'functions' 'initvalues' remain unchanged ) - enter x-range : :xmin: -2 :xmax: 3 - press "ENTER" ( leave 'h' unchanged ) Calculated values are quite accurate untill the program passes the discontinuity x = -1. Values beyond x = -1 are calculated and stored in äDAT but are not correct anymore. If you now change stepsize 'h' to 0.5 you'll get in stack: 1: program killed: "infinite result" and stuff in other levels. This error is caused by y'' when y-values for x = 0 were calculated. But values already stored in äDAT are correct. ÉÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ» º Example #3 : º ÈÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍͼ ------------------------------------------------------------------------- deq-system : u''(t) = -v'(t) ³ u(0) = 0 ; u'(0) = 1 v''(t) = 4 * u'(t) ³ v(0) = 0 ; v'(0) = 0 solution : u(t) = 0.5 * sin(2*t) v(t) = 1 - cos(2*t) ------------------------------------------------------------------------- substitution for the program : Replace all functions except the derivatives of highest degree. t -> X u (t) -> Y0 u'(t) -> Y1 v (t) -> Y2 v'(t) -> Y3 with these new variables you get : original deq-system ³ new deq-system for calculator ÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÅÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ ët(u) = u' ³ ëX (Y0) = Y1 ët(u') = u'' = - v' ³ ëX (Y1) = - Y3 ët(v) = v' ³ ëX (Y2) = Y3 ët(v') = v'' = 4 * u' ³ ëX (Y3) = 4 * Y1 - press "DEQ" - degree = 4 - enter the deq-system : :ëX(Y3)=Y4: '4*Y1' :ëX(Y2) : 'Y3' ( use PgDn ) :ëX(Y1) : '-Y3' :ëX(Y0) : 'Y1' - enter initial values : : x0 : 0 :Y0(x0): 0 ( use PgDn ) :Y1(x0): 1 :Y2(x0): 0 :Y3(x0): 0 - enter x-range : :xmin: 0 :xmax: 'ã' - enter stepsize : 'ã/10' u(t) is plotted at the end of the program and v(t) can be found in column 4 of äDAT. ********************************************************************** As I've got only some of the GOODIES disks I don't know if someone has already written a similar program. [Note: See also DIFFEQ on Goodies Disk #1 for another programmer's approach to differential equations. -jkh-] Although I didn't take much time to speed up the program I hope it can be useful to someone. Don't forget this is only a numerical guess process and therefore only useful to control the manually calculated result ( how to find out discontinuities or periodic times before getting to know y(x) ??). But nonsymbolic results can be calculated with any run-of-the-mill calculator as well. Frank Plate Im Hachegrund 3a 2808 Syke Germany - W