16 lines
438 B
Text
16 lines
438 B
Text
rk4(f,dx,x,y)={
|
|
my(k1=dx*f(x,y), k2=dx*f(x+dx/2,y+k1/2), k3=dx*f(x+dx/2,y+k2/2), k4=dx*f(x+dx,y+k3));
|
|
y + (k1 + 2*k2 + 2*k3 + k4) / 6
|
|
};
|
|
rate(x,y)=x*sqrt(y);
|
|
go()={
|
|
my(x0=0,x1=10,dx=.1,n=1+(x1-x0)\dx,y=vector(n));
|
|
y[1]=1;
|
|
for(i=2,n,y[i]=rk4(rate, dx, x0 + dx * (i - 1), y[i-1]));
|
|
print("x\ty\trel. err.\n------------");
|
|
forstep(i=1,n,10,
|
|
my(x=x0+dx*i,y2=(x^2/4+1)^2);
|
|
print(x "\t" y[i] "\t" y[i]/y2 - 1)
|
|
)
|
|
};
|
|
go()
|