'********************************************************************
'* Differential equations of order N *
'* by Runge-Kutta method of order 4 *
'* ---------------------------------------------------------------- *
'* Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc *
'* DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991" *
'* [BIBLI 03]. *
'* *
'* Basic version by J-P Moreau, Paris *
'* (www.jpmoreau.fr) *
'* ---------------------------------------------------------------- *
'* SAMPLE RUN: *
'* *
'* Example: integrate y"=(2y'y'+yy)/y from x=4 to x=6 *
'* with initial conditions: y(4)=2 and y'(4)=-2tan(1) *
'* *
'* Exact solution is: y = 2cos(1)/cos(x-5) *
'* *
'* DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER N *
'* of type y(n) = f(y(n-1),...,y',y,x) *
'* *
'* order of equation: ? 2 *
'* begin value x : ? 4 *
'* end value x : ? 6 *
'* y0 value at x0 : ? 2 *
'* y1 value at x0 : ? -3.114815 *
'* number of points : ? 11 *
'* finesse : ? 20 *
'* *
'* X Y *
'* --------------------------- *
'* 4.000000 2.000000 *
'* 4.200000 1.551018 *
'* 4.400000 1.309291 *
'* 4.600000 1.173217 *
'* 4.800000 1.102583 *
'* 5.000000 1.080605 *
'* 5.200000 1.102583 *
'* 5.400000 1.173217 *
'* 5.600000 1.309291 *
'* 5.800000 1.551018 *
'* 6.000000 2.000000 *
'* *
'********************************************************************
defint i-n
defdbl a-h,o-z
'ifi,i,ndata,iordre: INTEGER
dim yi(10), t(50)
dim ta(10),tb(10),tc(10),td(10),y(10),z(10)
Cls
print
print " DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER N"
print " of type y(n) = f(y(n-1),...,y',y,x)"
print
print " order of equation: "; : input iordre
print
print " begin value x : "; : input xi
print " end value x : "; : input xf
for i=0 to iordre-1
print " y";i;" value at x0 : "; : input yi(i)
next i
print " number of points : "; : input m
print " finesse : "; : input ifi
'call subroutine equadiffn
gosub 2000
END
'Example: y"=(2y'y'+yy)/y
1000 'FUNCTION fp(x,y())
if abs(y(0))<1e-12 then y(0)=1e-12
fp=(2*y(1)*y(1)+y(0)*y(0))/y(0)
return
'***************************************************************************
'* SOLVING DIFFERENTIAL EQUATIONS WITH 1 VARIABLE OF ORDER N *
'* of type y(n) = f(y(n-1),y(n-2),...,y',y,x) *
'* ----------------------------------------------------------------------- *
'* INPUTS: *
'* fp Equation to integrate *
'* xi, xf Begin, end values of variable x *
'* Yi Begin values at xi (of f and derivatives) *
'* m Number of points to calculate *
'* n Order of differential equation *
'* fi finesse (number of intermediate points) *
'* *
'* OUTPUTS: *
'* t real vector storing m results for function y *
'* ----------------------------------------------------------------------- *
'* EXAMPLE: y" = (2 y'y' + yy) / y with y(4) = 2, y'(4) = -2*tan(1) *
'* Exact solution: y = (2 cos(1))/cos(x-5) *
'***************************************************************************
2000 'Subroutine Equadiffn
'h,x,a,b,c,d : double
'ta,tb,tc,td,y,z : Table;
'i,j,k,ni,n1,n2 : integer
n=iordre
if ifi<1 then return
h = (xf - xi) / ifi / (m-1)
n1=n-1 : n2=n-2
t(1)=Yi(0)
for k=0 to n1
y(k)=Yi(k) : z(k)=Yi(k)
next k
for i=1 to m
ni=(i-1)*ifi-1
for j=1 to ifi
x=xi+h*(ni+j)
for k=0 to n1
y(k)=z(k)
next k
gosub 1000
a=h*fp
for k=0 to n2
ta(k)=h*y(k+1) : y(k)=z(k)+ta(k)/2#
next k
y(n1)=z(n1)+a/2#
x=x+h/2
gosub 1000
b=h*fp
for k=0 to n2
tb(k)=h*y(k+1) : y(k)=z(k)+tb(k)/2#
next k
y(n1)=z(n1)+b/2#
gosub 1000
c=h*fp
for k=0 to n2
tc(k)=h*y(k+1) : y(k)=z(k)+tc(k)
next k
y(n1)=z(n1)+c
x=x+h/2
gosub 1000
d=h*fp
for k=0 to n2
z(k)=z(k)+(ta(k)+2#*tb(k)+2#*tc(k)+h*y(k+1))/6#
next k
z(n1)=z(n1)+(a+b+b+c+c+d)/6#
next j
t(i+1)=z(0)
next i
'call subroutine Affiche
gosub 3000
Return
3000 'Subroutine Affiche
h=(xf-xi)/(m-1)
x=xi-h
Cls
print
print " X Y "
print "----------------------"
for i=1 to m
x=x+h
print using " ##.###### ##.######"; x; t(i)
next i
return
'End of file teqdifn.bas
Crear cuenta
