summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/examples/slope.asy
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/utils/asymptote/examples/slope.asy')
-rw-r--r--Build/source/utils/asymptote/examples/slope.asy83
1 files changed, 83 insertions, 0 deletions
diff --git a/Build/source/utils/asymptote/examples/slope.asy b/Build/source/utils/asymptote/examples/slope.asy
new file mode 100644
index 00000000000..832f8c13fc7
--- /dev/null
+++ b/Build/source/utils/asymptote/examples/slope.asy
@@ -0,0 +1,83 @@
+import ode;
+import graph;
+import math;
+size(200,200,IgnoreAspect);
+
+real f(real t, real y) {return cos(y);}
+//real f(real t, real y) {return 1/(1+y);}
+typedef real function(real,real);
+
+real a=0;
+real b=1;
+real y0=0;
+
+real L[]={1,2};
+
+int M=L.length; // Number of modes.
+
+//real Y0[]=array(M,y0);
+real Y0[]=new real[] {-1,2};
+
+real[] F(real t, real[] y) {
+ return sequence(new real(int m) {return f(t,y[M-m-1]);},M);
+ // return new real[] {exp((L[1]-1)*t)*y[1],
+ // -exp(-(L[1]-1)*t)*y[0]};
+ // return new real[]{-y[0]^2};
+}
+
+real[] G(real t, real[] y) {
+ return F(t,y)-sequence(new real(int m) {return L[m]*y[m];},M);
+}
+
+real lambda=sqrt(0.5);
+real[] tau,error,error2;
+int n=25;
+
+real order=3;
+
+for(int i=0; i < n-1; ++i) {
+ real dt=(b-a)*lambda^(n-i);
+ Solution S=integrate(Y0,L,F,a,b,dt,dynamic=false,0.0002,0.0004,ERK3BS,verbose=false);
+ real maxnorm=0;
+
+ Solution E=integrate(Y0,G,a,b,1e-2*dt,dynamic=false,0.0002,0.0004,RK5);
+ real[] exact=E.y[E.y.length-1];
+
+ // real[] exact=new real[] {exp(-b)*sin(b),exp(-L[1]*b)*cos(b)};
+ for(int m=0; m < M; ++m)
+ maxnorm=max(maxnorm,abs(S.y[S.y.length-1][m]-exact[m]));
+ if(maxnorm != 0) {
+ tau.push(dt);
+ // error.push(dt^-(order+1)*maxnorm);
+ error.push(maxnorm);
+ }
+}
+
+/*
+ for(int i=0; i < n-1; ++i) {
+ real dt=(b-a)*lambda^(n-i);
+ real maxnorm=0;
+ for(int m=0; m < M; ++m) {
+ solution S=integrate(Y0[m],L[m],f,a,b,dt,dynamic=false,0.000,1000,RK4_375,verbose=false);
+ maxnorm=max(maxnorm,abs(S.y[S.y.length-1]-exact[m]));
+ }
+ error2.push(dt^-order*maxnorm);
+ }
+*/
+
+//scale(Log,Log);
+scale(Log,Linear);
+
+//draw(graph(tau,error),marker(scale(0.8mm)*unitcircle,red));
+//draw(graph(tau,error2),marker(scale(0.8mm)*unitcircle,blue));
+
+int[] index=sequence(error.length-1);
+real[] slope=log(error[index+1]/error[index])/log(tau[index+1]/tau[index]);
+real[] t=sqrt(tau[index]*tau[index+1]);
+//write(t,slope);
+draw(graph(t,slope),red);
+
+
+
+xaxis("$\tau$",BottomTop,LeftTicks);
+yaxis("$e/\tau^"+string(order)+"$",LeftRight,RightTicks);