/*************************************************************** * View the solutions of a differential equations system having * * the form: y'=f(x,y,z), z'=g(x,y,z) * * * * Visual C++ in API style version by J-P Moreau * * ------------------------------------------------------------ * * REFERENCE: * * "Graphisme dans le plan et dans l'espace avec Turbo Pascal * * 4.0 by R. Dony - MASSON, Paris 1990, pages 254-260 ". * * ------------------------------------------------------------ * * TUTORIAL: * * It can be proved that any differential equation of order > 1 * * may be replaced by a 1st order differential equations system.* * For example, the 2nd order equation y" = f(x,y,y') is equi- * * valent to the system: * * y' = z * * z' = f(x,y,z) * * * * This program draws the solutions (integral curves) of a sys- * * tem with the form: * * y' = f(x,y,z) * * z' = g(x,y,z) * * * * The proposed example is: y' = x*x*y*y - 1 * * z' = x*x+y*y - 4 * * * * The steps to follow are: * * * * 1. choose a physical coordinates window x1,x2,y1,y2 * * 2. choose an increment step, dx * * 3. for each point of plane (x,y) inside previous window, * * calculate by Runge-Kutta method and draw the solution * * curve y=f(x) passing through this point. * * * * Obviously, these curves will eventually go out of screen. It * * can be managed by the clipping capability of module Gr2D. * * However to limit computation time, a curve is stopped when: * * * * - the following point will exit graph window, * * - the number of segments already drawn > 300, * * - the curve is approaching a critical point where * * - f(x,y) and g(x,y) ==> 0. * * * * Other possible examples: * * * * 1) y' = x + y z' = x*y * * 2) y' = y(y*y-1) z' = x(x*x-1) * * 3) y' = -x z' = x*x + y*y - 1 * * 4) y' = cos(x)-x*sin(y) z' = x*sin(y) + y*sin(x) * * * * ------------------------------------------------------------ * * NOTE: the variable switch (in Runge-Kutta procedure) is used * * to "switch" direction. For each start point in plane * * (x,y), the integral curve can be drawn in two opposite * * directions. When switch=1, the curve is drawn in one * * direction until a stop condition is met; at this moment* * switch is forced to -1 and the other part of the curve * * is drawn until a new stop condition is met. Then switch* * is again put to 1 and the process will continue with a * * new starting point. * ***************************************************************/ //Use with gr2d.cpp and rk.mak #include #include #include #include HDC hdc; RECT rect; HPEN hpen, hpenOld; //Functions used here of module Gr2D.cpp void Fenetre(double,double,double,double); void Cloture(int,int,int,int); void Bordure(); void Axes(); void Gradue(double,double); void Grille(double,double); void MoveXY(double,double); void LineXY(double,double); //"home made" graphic commands for hdc environment used by above functions void DrawPixel(int ix,int iy) { //sorry, no other available command found Rectangle(hdc,rect.left+ix,rect.top+iy, rect.left+ix+2,rect.top+iy+1); } void Swap(int *i1,int *i2) { int it; it=*i1; *i1=*i2; *i2=it; } void DrawLine(int ix1,int iy1,int ix2,int iy2) { int i,il,ix,iy; float dx,dy; if (ix2=f1 && y1<=f2 && z1>=f3 && z1<=f4 && seg1e-4); iswitch=-iswitch; } while (iswitch!=1); } //RungeKutta() void Initgraph() { extern MaxX,MaxY; HPEN redpen; redpen = CreatePen(PS_SOLID, 1, RGB(255, 0, 0)); Fenetre(f1,f2,f3,f4); Cloture(60,MaxX-10,40,MaxY-10); Gradue((f2-f1)/10,(f4-f3)/10); SelectObject(hdc,redpen); Grille((f2-f1)/10,(f4-f3)/10); SelectObject(hdc,hpen); Axes(); Bordure(); } //visualize Runge-Kutta solutions void Draw_RungeKutta() { extern MaxX, MaxY; //client drawing zone in pixels MaxX=(rect.right-rect.left); MaxY=(rect.bottom-rect.top); Data(); Initgraph(); i=0.0; while (i < 1) { i+=h; xpos=f1+i*(f2-f1); j=0; while (j < 1) { j+=h; ypos=f3+j*(f4-f3); RungeKutta(xpos,ypos); } } } //Draw_RungeKutta() //Handle window Paint and Destroy messages long FAR PASCAL /*_export*/ WndProc(HWND hwnd, UINT message, UINT wParam, LONG lParam) { PAINTSTRUCT ps; switch (message) { case WM_PAINT: hdc = BeginPaint(hwnd, &ps); GetClientRect(hwnd, &rect); hpen = CreatePen(PS_SOLID, 1, RGB(0, 0, 255)); hpenOld = SelectObject(hdc,hpen); Draw_RungeKutta(); SelectObject(hdc,hpenOld); DeleteObject(hpen); EndPaint(hwnd, &ps); return 0; case WM_DESTROY: PostQuitMessage(0); return 0; } return DefWindowProc(hwnd, message, wParam, lParam); } // main program int PASCAL WinMain(HANDLE hInstance, HANDLE hPrevInstance, LPSTR lpszCmdParam, int nCmdShow) { static char szAppName[] = "RK"; HWND hwnd; MSG msg; WNDCLASS wndclass; if (!hPrevInstance) { wndclass.style = CS_HREDRAW | CS_VREDRAW; wndclass.lpfnWndProc = WndProc; wndclass.cbClsExtra = 0; wndclass.cbWndExtra = 0; wndclass.hInstance = hInstance; wndclass.hIcon = LoadIcon(NULL, IDI_APPLICATION); wndclass.hCursor = LoadCursor(NULL, IDC_ARROW); wndclass.hbrBackground = GetStockObject(WHITE_BRUSH); wndclass.lpszMenuName = NULL; wndclass.lpszClassName = szAppName; RegisterClass(&wndclass); } hwnd = CreateWindow(szAppName, // window class name " RUNGE KUTTA", // window caption WS_OVERLAPPEDWINDOW, // window style CW_USEDEFAULT, // initial x position CW_USEDEFAULT, // initial y position CW_USEDEFAULT, // initial x size CW_USEDEFAULT, // initial y size NULL, // parent window handle NULL, // window menu handle hInstance, // program instance handle NULL); // creation parameters ShowWindow(hwnd, nCmdShow); UpdateWindow(hwnd); while (GetMessage(&msg, NULL, 0, 0)) { TranslateMessage(&msg); DispatchMessage(&msg); } return msg.wParam; } //end of file rk.cpp