From: Pat Thoyts Date: Tue, 5 Aug 2014 14:03:19 +0000 (+0100) Subject: C# port of the BASIC second-order non-linear equation solver. X-Git-Url: http://privyetmir.co.uk/gitweb.cgi?a=commitdiff_plain;h=df2a8d500089e77891ff76d3e4f8e5fc4fae5f22;p=MervsSolver C# port of the BASIC second-order non-linear equation solver. Validated by running the original code on the VICE emulator. A unit test compares the generated data against the emulator generated values. Signed-off-by: Pat Thoyts --- diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..50b22b1 --- /dev/null +++ b/.gitignore @@ -0,0 +1,9 @@ +/MervsSolver/bin +/MervsSolver/obj +/SolverTest/bin +/SolverTest/obj +/SolverTest/test-results/ +*.userprefs +*.prg +*~ + diff --git a/MervsSolver.sln b/MervsSolver.sln new file mode 100644 index 0000000..acf8d77 --- /dev/null +++ b/MervsSolver.sln @@ -0,0 +1,26 @@ + +Microsoft Visual Studio Solution File, Format Version 11.00 +# Visual Studio 2010 +Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "MervsSolver", "MervsSolver\MervsSolver.csproj", "{BFBAA38C-7F6F-44BB-A4DA-C943818FDD3A}" +EndProject +Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "SolverTest", "SolverTest\SolverTest.csproj", "{459AC64B-AA17-492D-80E4-A4A91AF14296}" +EndProject +Global + GlobalSection(SolutionConfigurationPlatforms) = preSolution + Debug|x86 = Debug|x86 + Release|x86 = Release|x86 + EndGlobalSection + GlobalSection(ProjectConfigurationPlatforms) = postSolution + {459AC64B-AA17-492D-80E4-A4A91AF14296}.Debug|x86.ActiveCfg = Debug|Any CPU + {459AC64B-AA17-492D-80E4-A4A91AF14296}.Debug|x86.Build.0 = Debug|Any CPU + {459AC64B-AA17-492D-80E4-A4A91AF14296}.Release|x86.ActiveCfg = Release|Any CPU + {459AC64B-AA17-492D-80E4-A4A91AF14296}.Release|x86.Build.0 = Release|Any CPU + {BFBAA38C-7F6F-44BB-A4DA-C943818FDD3A}.Debug|x86.ActiveCfg = Debug|x86 + {BFBAA38C-7F6F-44BB-A4DA-C943818FDD3A}.Debug|x86.Build.0 = Debug|x86 + {BFBAA38C-7F6F-44BB-A4DA-C943818FDD3A}.Release|x86.ActiveCfg = Release|x86 + {BFBAA38C-7F6F-44BB-A4DA-C943818FDD3A}.Release|x86.Build.0 = Release|x86 + EndGlobalSection + GlobalSection(MonoDevelopProperties) = preSolution + StartupItem = MervsSolver\MervsSolver.csproj + EndGlobalSection +EndGlobal diff --git a/MervsSolver/Coord.cs b/MervsSolver/Coord.cs new file mode 100644 index 0000000..8c20672 --- /dev/null +++ b/MervsSolver/Coord.cs @@ -0,0 +1,15 @@ +using System; + +namespace MervsSolver +{ + public class Coord + { + public double X { get; set;} + public double Y { get; set; } + public Coord(double x, double y) + { + X = x; + Y = y; + } + } +} diff --git a/MervsSolver/MervsSolver.csproj b/MervsSolver/MervsSolver.csproj new file mode 100644 index 0000000..0ff94ec --- /dev/null +++ b/MervsSolver/MervsSolver.csproj @@ -0,0 +1,44 @@ + + + + Debug + x86 + 10.0.0 + 2.0 + {BFBAA38C-7F6F-44BB-A4DA-C943818FDD3A} + Exe + MervsSolver + MervsSolver + + + true + full + false + bin\Debug + DEBUG; + prompt + 4 + true + x86 + .1 .1 .1 .1 + + + full + true + bin\Release + prompt + 4 + true + x86 + + + + + + + + + + + + \ No newline at end of file diff --git a/MervsSolver/Program.cs b/MervsSolver/Program.cs new file mode 100644 index 0000000..09156eb --- /dev/null +++ b/MervsSolver/Program.cs @@ -0,0 +1,49 @@ +using System; +using System.Collections.Generic; +using System.Diagnostics; +using System.IO; + +namespace MervsSolver +{ + class MainClass + { + public static void Main(string[] args) + { + double x0, v0, mu, steps; + if (Double.TryParse(args[0], out x0) + && Double.TryParse(args[1], out v0) + && Double.TryParse(args[2], out mu) + && Double.TryParse(args[3], out steps)) + { + Solver solver = new Solver(x0, v0, mu, steps); + Process gnuplot = Gnuplot(mu); + //Console.WriteLine("{0} results", solver.Result.Count); + int i = 0; + foreach (Coord coord in solver.Result) + { + gnuplot.StandardInput.WriteLine("{0:0.000}\t{1:0.000}", coord.X, coord.Y); + ++i; + } + gnuplot.StandardInput.Close(); + gnuplot.WaitForExit(); + gnuplot.Close(); + } + else + { + Console.Error.WriteLine("usage: {0} x0 v0 mu stepsize", "solver"); + } + } + public static Process Gnuplot(double mu) + { + string arg = String.Format("-persist -e \"plot \\\"-\\\" using 1:2 with lines title \\\"mu {0:0.000}\\\"\"", mu); + ProcessStartInfo si = new ProcessStartInfo(); + si.FileName = "gnuplot"; + si.Arguments = arg; + si.CreateNoWindow = true; + si.UseShellExecute = false; + si.RedirectStandardInput = true; + Process process = Process.Start(si); + return process; + } + } +} diff --git a/MervsSolver/Properties/AssemblyInfo.cs b/MervsSolver/Properties/AssemblyInfo.cs new file mode 100644 index 0000000..a7d591e --- /dev/null +++ b/MervsSolver/Properties/AssemblyInfo.cs @@ -0,0 +1,22 @@ +using System.Reflection; +using System.Runtime.CompilerServices; + +// Information about this assembly is defined by the following attributes. +// Change them to the values specific to your project. +[assembly: AssemblyTitle("MervsSolver")] +[assembly: AssemblyDescription("")] +[assembly: AssemblyConfiguration("")] +[assembly: AssemblyCompany("")] +[assembly: AssemblyProduct("")] +[assembly: AssemblyCopyright("pat")] +[assembly: AssemblyTrademark("")] +[assembly: AssemblyCulture("")] +// The assembly version has the format "{Major}.{Minor}.{Build}.{Revision}". +// The form "{Major}.{Minor}.*" will automatically update the build and revision, +// and "{Major}.{Minor}.{Build}.*" will update just the revision. +[assembly: AssemblyVersion("1.0.*")] +// The following attributes are used to specify the signing key for the assembly, +// if desired. See the Mono documentation for more information about signing. +//[assembly: AssemblyDelaySign(false)] +//[assembly: AssemblyKeyFile("")] + diff --git a/MervsSolver/Solver.cs b/MervsSolver/Solver.cs new file mode 100644 index 0000000..9f5db28 --- /dev/null +++ b/MervsSolver/Solver.cs @@ -0,0 +1,83 @@ +using System; +using System.Collections.Generic; + +namespace MervsSolver +{ + public class Solver + { + List mResult = new List(); + public List Result { get { return mResult; } } + + public Solver(double initialDisplacement, double initialVelocity, double continuousParameter, double stepSize) + { + double x0 = initialDisplacement; + double v0 = initialVelocity; + double mu = continuousParameter; + double d = stepSize; + + double[] A = new double[300]; + double[] E = new double[300]; + double b = 0, c = 0, z = 0, s = 0, h = 0, x1 = 0, v1 = 0; + int i = 0; + + if (v0 < 0) + goto Line320; + if (v0 > 0) + goto Line300; + if (x0 > 0) + goto Line320; + Line300: + z = d; + goto Line330; + Line320: + z = -d; + Line330: + x1 = x0 + z; + b = equation(v0, z, mu, x1); + c = z * x1; + s = (b * b) - 4 * c; + if (s < 0) + goto Line540; + double k = z / Math.Abs(z); + v1 = (-b + k * Math.Sqrt((b * b) - 4 * c)) / 2; + Line430: + i = i + 1; + A[i] = x1; + E[i] = v1; + h = h + 1; + if (h == 50) + goto Line590; + Line480: + v0 = v1; + x0 = x1; + if (v0 == 0) + goto Line520; + goto Line330; + Line520: + if (i == 0 || i > 500) + return; + if (E[i - 1] < 0) + goto Line300; + goto Line320; + Line540: + z = z / 10; + if (Math.Abs(z) < 1.0e-3) + goto Line570; + goto Line330; + Line570: + v1 = 0; + goto Line430; + Line590: + for (i = 1; i != h; ++i) + mResult.Add(new Coord(0 + (A[i] * -75), (E[i] * -75))); + h = 0; i = 0; + goto Line480; + } + + static double equation(double v0, double z, double mu, double x1) + { + double b = -(v0 + z * (mu * (1 - (x1*x1)))); + return b; + } + } +} diff --git a/README.txt b/README.txt new file mode 100644 index 0000000..691b84f --- /dev/null +++ b/README.txt @@ -0,0 +1,19 @@ +This is a non-linear 2nd order equation solver for + x' - mu * (1-x**2) * x' + x = 0 + +solver.bas is the original BASIC 3.5 program + +solver2.bas is modified to print the results to a printer which under +the emulator yields a text file with the result as columns of numbers. + +To run the program, convert the basic program to a .prg file: + petcat -w3 -o solver2.prg solver2.bas + +Then launch the emulator: + xplus4 -device4 1 -prtxtdev1 output.txt solver2.prg + +This command attaches the emulator printer to the local text file +"output.txt" which we can then use for plotting with gnuplot: + plot "output.txt" 1:2 with lines + + diff --git a/SolverTest/SolverTest.csproj b/SolverTest/SolverTest.csproj new file mode 100644 index 0000000..dddebaa --- /dev/null +++ b/SolverTest/SolverTest.csproj @@ -0,0 +1,52 @@ + + + + Debug + AnyCPU + 10.0.0 + 2.0 + {459AC64B-AA17-492D-80E4-A4A91AF14296} + Library + SolverTest + SolverTest + + + true + full + false + bin\Debug + DEBUG; + prompt + 4 + false + + + full + true + bin\Release + prompt + 4 + false + + + + + False + + + + + + {BFBAA38C-7F6F-44BB-A4DA-C943818FDD3A} + MervsSolver + + + + + PreserveNewest + + + + + + \ No newline at end of file diff --git a/SolverTest/Test.cs b/SolverTest/Test.cs new file mode 100644 index 0000000..7517495 --- /dev/null +++ b/SolverTest/Test.cs @@ -0,0 +1,32 @@ +using NUnit.Framework; +using System; +using System.Collections.Generic; +using System.Linq; +using System.IO; +using MervsSolver; + +namespace SolverTest +{ + [TestFixture()] + public class Test + { + [Test()] + public void solver_regression_test() + { + Solver solver = new Solver(0.1, 0.1, 0.1, 0.1); + string[] lines = File.ReadAllLines("test.data"); + Assert.That(solver.Result.Count == lines.Length, + "generated less results than expected ({0} should be {1})", solver.Result.Count, lines.Length); + for (int n = 0; n < lines.Length; ++n) + { + string[] bits = lines[n].Split(new char[0], StringSplitOptions.RemoveEmptyEntries); + double x, y; + Assert.IsTrue(double.TryParse(bits[0], out x)); + Assert.IsTrue(double.TryParse(bits[1], out y)); + Assert.That(Math.Abs(x - solver.Result[n].X) < 0.001); + Assert.That(Math.Abs(y - solver.Result[n].Y) < 0.001); + } + } + } +} + diff --git a/SolverTest/test.data b/SolverTest/test.data new file mode 100644 index 0000000..364537f --- /dev/null +++ b/SolverTest/test.data @@ -0,0 +1,343 @@ +-8.250 -6.643 +-9.000 -5.486 +-9.750 -3.425 +-9.825 -3.202 +-9.900 -2.958 +-9.975 -2.687 +-10.050 -2.378 +-10.125 -2.006 +-10.200 -1.506 +-10.275 0.000 +-2.775 4.952 +-2.025 5.313 +-1.275 5.560 +-0.525 5.704 +0.225 5.749 +0.975 5.696 +1.725 5.537 +2.475 5.259 +3.225 4.834 +3.975 4.198 +4.725 3.147 +4.800 3.036 +4.875 2.918 +4.950 2.793 +5.025 2.658 +5.100 2.514 +5.175 2.356 +5.250 2.184 +5.325 1.990 +5.400 1.769 +5.475 1.503 +5.550 1.148 +5.625 0.000 +4.875 -1.950 +4.125 -3.042 +3.375 -3.785 +2.625 -4.316 +1.875 -4.691 +1.125 -4.937 +0.375 -5.067 +-0.375 -5.087 +-1.125 -4.993 +-1.875 -4.773 +-2.625 -4.401 +-3.375 -3.812 +-4.125 -2.770 +-4.200 -2.659 +-4.275 -2.540 +-4.350 -2.412 +-4.500 -2.122 +-4.575 -1.954 +-4.650 -1.764 +-4.725 -1.541 +-4.800 -1.264 +-4.875 -0.831 +-4.950 0.000 +-4.200 1.813 +-3.450 2.809 +-2.700 3.468 +-1.950 3.916 +-1.200 4.205 +-0.450 4.357 +0.300 4.381 +1.050 4.272 +1.800 4.010 +2.550 3.546 +3.300 2.706 +3.375 2.616 +3.450 2.521 +3.525 2.419 +3.600 2.310 +3.675 2.192 +3.750 2.063 +3.825 1.921 +3.900 1.763 +3.975 1.582 +4.050 1.367 +4.125 1.091 +4.200 0.000 +3.450 -1.646 +2.700 -2.524 +1.950 -3.074 +1.200 -3.413 +0.450 -3.582 +-0.300 -3.595 +-1.050 -3.441 +-1.800 -3.077 +-2.550 -2.332 +-2.625 -2.252 +-2.700 -2.166 +-2.775 -2.073 +-2.850 -1.972 +-2.925 -1.862 +-3.000 -1.740 +-3.075 -1.603 +-3.150 -1.448 +-3.225 -1.264 +-3.300 -1.031 +-3.450 0.000 +-2.700 1.461 +-1.950 2.201 +-1.200 2.619 +-0.450 2.814 +0.300 2.809 +1.050 2.579 +1.800 1.967 +1.875 1.901 +1.950 1.828 +2.025 1.749 +2.100 1.662 +2.175 1.565 +2.250 1.457 +2.325 1.333 +2.400 1.189 +2.475 1.014 +2.550 0.774 +2.625 0.000 +1.875 -1.224 +1.125 -1.774 +0.375 -1.991 +-0.375 -1.919 +-1.125 -1.385 +-1.200 -1.324 +-1.275 -1.256 +-1.350 -1.177 +-1.425 -1.086 +-1.500 -0.979 +-1.575 -0.847 +-1.650 -0.670 +-1.725 0.000 +-0.975 0.893 +-0.225 1.119 +-0.150 1.137 +-0.075 1.149 +0.000 1.157 +0.075 1.159 +0.150 1.157 +0.225 1.150 +0.300 1.137 +0.375 1.120 +0.450 1.097 +0.525 1.067 +0.600 1.031 +0.675 0.987 +0.750 0.935 +0.825 0.871 +0.900 0.793 +1.050 0.564 +1.125 0.000 +0.375 -0.569 +0.300 -0.613 +0.225 -0.647 +0.150 -0.671 +0.075 -0.687 +0.000 -0.694 +-0.075 -0.694 +-0.150 -0.685 +-0.225 -0.667 +-0.300 -0.639 +-0.375 -0.600 +-0.450 -0.546 +-0.525 -0.469 +-0.600 -0.347 +-0.675 0.000 +-0.600 0.216 +-0.525 0.339 +-0.450 0.426 +-0.375 0.491 +-0.300 0.540 +-0.225 0.577 +-0.150 0.603 +-0.075 0.620 +0.000 0.627 +0.075 0.626 +0.150 0.615 +0.225 0.594 +0.300 0.561 +0.375 0.514 +0.450 0.446 +0.525 0.336 +0.600 0.000 +0.525 -0.202 +0.450 -0.316 +0.375 -0.395 +0.300 -0.452 +0.225 -0.494 +0.150 -0.523 +0.075 -0.541 +0.000 -0.548 +-0.075 -0.546 +-0.150 -0.532 +-0.225 -0.506 +-0.300 -0.465 +-0.375 -0.403 +-0.450 -0.297 +-0.525 0.000 +-0.375 0.291 +-0.300 0.361 +-0.225 0.410 +-0.150 0.443 +-0.075 0.462 +0.000 0.470 +0.075 0.465 +0.150 0.448 +0.225 0.415 +0.300 0.359 +0.375 0.258 +0.450 0.000 +0.375 -0.171 +0.300 -0.264 +0.225 -0.324 +0.150 -0.362 +0.075 -0.384 +0.000 -0.392 +-0.075 -0.385 +-0.150 -0.361 +-0.225 -0.315 +-0.300 -0.221 +-0.375 0.000 +-0.300 0.154 +-0.225 0.234 +-0.150 0.281 +-0.075 0.307 +0.000 0.314 +0.075 0.303 +0.150 0.269 +0.225 0.186 +0.300 0.000 +0.225 -0.134 +0.150 -0.198 +0.075 -0.230 +0.000 -0.237 +-0.075 -0.219 +-0.150 -0.154 +-0.225 0.000 +-0.150 0.110 +-0.075 0.154 +0.000 0.161 +0.075 0.123 +0.150 0.000 +0.075 -0.079 +0.000 -0.086 +-0.075 0.000 +0.000 0.007 +0.075 0.000 +-0.075 0.000 +0.000 0.007 +0.075 0.000 +0.000 -0.008 +-0.075 0.000 +0.000 0.007 +0.075 0.000 +0.000 -0.008 +-0.075 0.000 +0.000 0.007 +0.075 0.000 +0.000 -0.008 +-0.075 0.000 +0.000 0.007 +0.075 0.000 +0.000 -0.008 +-0.075 0.000 +0.000 0.007 +0.075 0.000 +0.000 -0.008 +-0.075 0.000 +0.000 0.007 +0.075 0.000 +0.000 -0.008 +-0.075 0.000 +0.000 0.007 +0.075 0.000 +0.000 -0.008 +-0.075 0.000 +0.000 0.007 +0.075 0.000 +0.000 -0.008 +-0.075 0.000 +0.000 0.007 +0.075 0.000 +0.000 -0.008 +-0.075 0.000 +0.000 0.007 +0.075 0.000 +0.000 -0.008 +-0.075 0.000 +0.000 0.007 +0.075 0.000 +0.000 -0.008 +-0.075 0.000 +0.000 0.007 +0.075 0.000 +0.000 -0.008 +-0.075 0.000 +0.075 0.000 +0.150 0.000 +0.225 0.000 +0.300 0.000 +0.375 0.000 +0.450 0.000 +0.525 0.000 +0.600 0.000 +0.675 0.000 +0.750 0.000 +0.825 0.000 +0.900 0.000 +0.975 0.000 +1.050 0.000 +1.125 0.000 +1.200 0.000 +1.275 0.000 +1.350 0.000 +1.425 0.000 +1.500 0.000 +1.575 0.000 +1.650 0.000 +1.725 0.000 +1.800 0.000 +1.875 0.000 +1.950 0.000 +2.025 0.000 +2.100 0.000 +2.175 0.000 +2.250 0.000 +2.325 0.000 +2.400 0.000 +2.475 0.000 +2.550 0.000 +2.625 0.000 +2.700 0.000 +2.775 0.000 +2.850 0.000 +2.925 0.000 +3.000 0.000 +3.075 0.000 +3.150 0.000 +3.225 0.000 +3.300 0.000 +3.375 0.000 +3.450 0.000 +3.525 0.000 +3.600 0.000 +3.675 0.000 diff --git a/solver2.bas b/solver2.bas new file mode 100644 index 0000000..a890382 --- /dev/null +++ b/solver2.bas @@ -0,0 +1,76 @@ +10 rem y.h.ku non-linear equation solver +20 rem enter equation at line 800 +30 rem color 1,1:color 0,2,4:color 4,1:scnclr +40 rem open3,6:open5,3,3:open2,6,2:open1,6,1 +45 rem print#1,"h" +46 open 4,4,7 : rem open printer lowercase mode +50 print"***************************************" +60 print"* *" +70 print"*non-linear 2nd. order equation solver*" +80 print"* *" +90 print"* using y.h.ku's acc. phase plane *" +100 print"* *" +110 print"* m.k.hobden,plus-4 version may 89 *" +120 print"* *" +130 print"***************************************" +135 print +140 print" enter d.e at line 800." +150 print +160 print" enter initial velocity and displacement" +170 print" at the program prompts" +180 print +190 input" initial disp.(x0)";x0 +200 print +210 input" initial vel. (v0)";v0 +220 print +230 input" cont. param. (mu)";mu +240 print +245 input" step size (<1)";d +250 dim a(300),e(300) +260 print:print" calculating......" +270 if v0<0 then 320 +280 if v0>0 then 300 +290 if x0>0 then 320 +300 z=d +310 goto 330 +320 z=-d +330 x1=x0+z +340 gosub 800 +350 c=z*x1 +360 s=(b*b)-4*c +370 if s<0 then 540 : rem if elliptic +380 k=z/abs(z) +390 v1=(-b+k*sqr((b*b)-4*c))/2 : rem (b+sqrt(b**2-4ac))/2a +400 if g%=1 then 430 +410 if int(h/50)<>h/50 then 430 +430 i=i+1 +440 a(i)=x1 +450 e(i)=v1 +460 h=h+1 +470 if h=50 then 590 +480 v0=v1 +490 x0=x1 +500 if v0=0 then 520 +510 goto 330 +520 if e(i-1)<0 then 300 +530 goto 320 +540 z=z/10 +550 if abs(z)<1.0e-3 then 570 +560 goto 330 +570 v1=0 +580 goto 430 +590 if j>0 then 660 +600 rem scnclr:print:print +610 print#4,q$;" for mu =";mu +660 rem if g%=1 then 720 +720 if g%=0 then gosub 1000 +730 for i=1 to h +740 print#4,0+(a(i)*-75),(e(i)*-75) +750 next i +760 h=0:i=0:j=j+1:goto 480 +800 b=-(v0+z*(mu*(1-(x1*x1)))) +810 q$="x'-mu(1-(x*x))x'+x=0" +820 return +1000 rem print#4,"r",0+(a(1)*-75),(e(1)*-75) +1100 g%=1:return +1200 end