code10.cpp

ordinary differential equations

#include "Code10.h"

CMyWinApp MyApplication;

BOOL CMyWinApp::InitInstance()
{
    CCode10* pFrame = new CCode10;
    m_pMainWnd = pFrame;
    pFrame->ShowWindow(SW_SHOW);
    pFrame->UpdateWindow();
    return TRUE;
}

BEGIN_MESSAGE_MAP(CCode10,CFrameWnd)
	ON_WM_PAINT()
	ON_WM_LBUTTONDOWN()
    ON_BN_CLICKED(IDC_BUTTON,OnButton)
END_MESSAGE_MAP()

CCode10::CCode10()
{
	Create(NULL,"Code10: Ordinary Differential Equations",
		WS_OVERLAPPEDWINDOW,CRect(0,0,800,635),NULL);
	Arial80.CreatePointFont(80,"Arial");
	pt=new PT [M+1];
	fMenu=0; fStatus=0; idc=301;
	curve.hm=CPoint(330,320); curve.end=CPoint(750,570);
	curve.rc=CRect(curve.hm.x-10,curve.hm.y-10,
		curve.end.x+30,curve.end.y+10);
	menu[1].item="ODE 1: Taylor (Order 3)"; 
	menu[2].item="ODE 1: Runge-Kutta 2";
	menu[3].item="ODE 1: Runge-Kutta 4"; 
	menu[4].item="ODE 1: Adams-Bashforth";
	menu[5].item="ODE 1: System"; 
	menu[6].item="ODE 2: Initial Value Problem";
	menu[7].item="ODE 2: Finite-Difference 1"; 
	menu[8].item="ODE 2: Finite-Difference 2";
	for (int i=1;i<=nMenuItems;i++)
	{
		menu[i].hm=CPoint(20,30+(i-1)*30);
		menu[i].rc=CRect(menu[i].hm.x,menu[i].hm.y,
			menu[i].hm.x+150,menu[i].hm.y+20);
	}
	input[0].hm=CPoint(200,10);
	input[0].rc=CRect(input[0].hm.x,input[0].hm.y,
		input[0].hm.x+560,input[0].hm.y+280);
	for (i=1;i<=maxInput;i++)
		input[i].hm=CPoint(input[0].hm.x+10,input[0].hm.y+30+(i-1)*25);
}

CCode10::~CCode10()
{	
	delete pt;
}

void CCode10::OnPaint()
{
	CPaintDC dc(this);
	CString str;
	dc.SelectObject(Arial80);
	dc.SetBkColor(RGB(150,150,150));
	dc.SetTextColor(RGB(255,255,255));
	for (int i=1;i<=nMenuItems;i++)
	{
		dc.FillSolidRect(&menu[i].rc,RGB(150,150,150));
		dc.TextOut(menu[i].hm.x+5,menu[i].hm.y+5,menu[i].item);
	}
	dc.SetBkColor(RGB(255,255,255));
	dc.SetTextColor(RGB(100,100,100));
	dc.Rectangle(input[0].rc);
	dc.TextOut(input[0].hm.x+150,input[0].hm.y+5,menu[fMenu].item);
	if (fMenu>0)
	{
		for (i=1;i<=nInputItems;i++)
			dc.TextOut(input[i].hm.x+10,input[i].hm.y,input[i].label);
		if (fStatus==2)
		{
			DrawCurve();
			fStatus=1;
		}
	}
	if (fMenu==2)
		dc.TextOut(input[6].hm.x+10,input[6].hm.y,input[6].label);
	if (fMenu==7)
		dc.TextOut(input[0].hm.x+300,input[0].hm.y+5,
		"p(x) y\"+q(x) y'+r(x) y = w(x),  x[0]<=x<=x[m]");
	if (fMenu==8)
		dc.TextOut(input[0].hm.x+300,input[0].hm.y+5,
		"p(x) y\"+q(x) y'+r(x) y = w(x),  x[0]<=x<=x[m]");
}

void CCode10::DrawCurve()
{
	CClientDC dc(this);
	int i;
	double m1,m2,c1,c2;
	CString str;
	CPoint px,px2;
	PT top,bottom,max,min,left,right;
	Clear(curve.rc);
	left.x=pt[0].x; right.x=pt[m].x;
	left.y=pt[0].y;
	max.y=pt[0].y; min.y=pt[0].y;
	
	if (fMenu==5)
	{
		max.z=pt[0].z; 
		min.z=pt[0].z;
	}
	for (i=0;i<=m;i++)
	{
		if (max.y<pt[i].y)
			max.y=pt[i].y;	
		if (min.y>pt[i].y)
			min.y=pt[i].y;
		if (fMenu==5)
		{
			if (max.z<pt[i].z)
				max.z=pt[i].z;
			if (min.z>pt[i].z)
				min.z=pt[i].z;
		}
	}
	top.y=max.y; bottom.y=min.y;
	if (fMenu==5)
	{
		if (max.y<max.z)
			top.y=max.z;
		if (min.y>min.z)
			bottom.y=min.z;
	}

	// Cartesian-Windows conversion coordinates
	m1=(double)(curve.end.x-curve.hm.x)/(right.x-left.x);
	c1=(double)curve.hm.x-left.x*m1;
	m2=(double)(curve.hm.y-curve.end.y)/(top.y-bottom.y);
	c2=(double)curve.end.y-bottom.y*m2;

	// Draw & label the x,y axes
	CPen pGray(PS_SOLID,1,RGB(100,100,100));
	dc.SelectObject(pGray);
	dc.SelectObject(Arial80);
	px=CPoint(m1*0+c1,m2*bottom.y+c2); dc.MoveTo(px);
	px=CPoint(m1*0+c1,m2*top.y+c2); dc.LineTo(px);
	
	px=CPoint(m1*left.x+c1,m2*0+c2); dc.MoveTo(px);
	str.Format("%.0lf",left.x); dc.TextOut(px.x,px.y,str);
	px=CPoint(m1*right.x+c1,m2*0+c2); dc.LineTo(px);
	str.Format("%.0lf",right.x); dc.TextOut(px.x-10,px.y,str);

	// draw the curve
	CPen pDark(PS_SOLID,2,RGB(50,50,50));
	dc.SelectObject(pDark);
	//str.Format("(%.1lf,%.1lf)",left.x,left.y);
	//dc.TextOut((int)(m1*left.x+c1),(int)(m2*left.y+c2),str);
	for (i=0;i<=m;i++)
	{
		px=CPoint((int)(m1*pt[i].x+c1),(int)(m2*pt[i].y+c2));
		if (i==0)
			dc.MoveTo(px);
		else
			dc.LineTo(px);
		if (pt[i].y==top.y)
		{
			str.Format("%.3lf",top.y); 
			dc.TextOut(px.x,px.y-10,str);
		}
		if (pt[i].y==bottom.y)
		{
			str.Format("%.3lf",bottom.y); 
			dc.TextOut(px.x,px.y,str);
		}
	}
	if (fMenu==5)
	{
		CPen pLight(PS_SOLID,2,RGB(150,150,150));
		dc.SelectObject(pLight);
		for (i=0;i<=m;i++)
		{
			px2=CPoint((int)(m1*pt[i].x+c1),(int)(m2*pt[i].z+c2));
			if (i==0)
			{
				dc.MoveTo(px2);
				dc.TextOut(px2.x+15,px2.y-10,"z=g(x)");
			}
			else
				dc.LineTo(px2);
		}
	}
}

void CCode10::OnLButtonDown(UINT nFlags,CPoint pt)
{
	int i,k;
	for (k=1;k<=nMenuItems;k++)
		if (menu[k].rc.PtInRect(pt))
		{
			fMenu=k; fStatus=1;
			switch(k)
			{
				case 1: 
					nInputItems=7; 
					input[1].label="y'=g(x,y)";
					input[2].label="x[0]"; 
					input[3].label="x[m]";
					input[4].label="h"; 
					input[5].label="y[0]"; 
					input[6].label="y'' (u=y')";
					input[7].label="y'''(v=y'')";
					break;
				case 2: 
					nInputItems=6; 
					input[1].label="y'=g(x,y)"; 
					input[2].label="x[0]"; 
					input[3].label="x[m]";
					input[4].label="h"; 
					input[5].label="y[0]"; 
					input[6].label="r (0<r<=1)";
					break;
				case 3: 
					nInputItems=5; 
					input[1].label="y'=g(x,y)"; 
					input[2].label="x[0]"; 
					input[3].label="x[m]";
					input[4].label="h"; 
					input[5].label="y[0]"; 
					break;
				case 4: 
					nInputItems=5; 
					input[1].label="y'=g(x,y)"; 
					input[2].label="x[0]"; 
					input[3].label="x[m]";
					input[4].label="h"; 
					input[5].label="y[0]"; 
					break;
				case 5: 
					nInputItems=7; 
					input[1].label="y'=f(x,y,z)"; 
					input[2].label="z'=g(x,y,z)"; 
					input[3].label="x[0]"; 
					input[4].label="y[0]";	
					input[5].label="z[0]";
					input[6].label="h"; 
					input[7].label="x[m]";
					break;
				case 6: 
					nInputItems=6; 
					input[1].label="y\"=f(x,y,z=y')"; 
					input[2].label="x[0]"; 
					input[3].label="y[0]"; 
					input[4].label="z[0]=y'[0]";
					input[5].label="h"; 
					input[6].label="x[m]";
					break;
				case 7: 
					nInputItems=9; 
					input[1].label="p(x)"; 
					input[2].label="q(x)"; 
					input[3].label="r(x)"; 
					input[4].label="w(x)";
					input[5].label="h"; 
					input[6].label="x[0]"; 
					input[7].label="x[m]";
					input[8].label="y[0]"; 
					input[9].label="y[m]";
					break;
				case 8: 
					nInputItems=9; 
					input[1].label="p(x)"; 
					input[2].label="q(x)"; 
					input[3].label="r(x)"; 
					input[4].label="w(x)";
					input[5].label="h"; 
					input[6].label="x[0]";
					input[7].label="x[m]"; 
					input[8].label="y'[0]";
					input[9].label="y'[m]";
					break;
			}
			for (i=1;i<=maxInput;i++)
				input[i].ed.DestroyWindow();
			for (i=1;i<=nInputItems;i++)
				input[i].ed.Create(WS_CHILD| WS_VISIBLE| WS_BORDER, 
					CRect(input[i].hm.x+100,input[i].hm.y,
					input[i].hm.x+520,input[i].hm.y+20),this,idc++);
			InvalidateRect(input[0].rc);
			Clear(curve.rc);
			table.DestroyWindow();
			btn.DestroyWindow();
			btn.Create("Compute",WS_CHILD| WS_VISIBLE| BS_DEFPUSHBUTTON,
				CRect(CPoint(input[0].hm.x+10,input[0].hm.y+5),
				CSize(100,20)),this,IDC_BUTTON);
		}
}

void CCode10::ODE1RK2()
{
	int i,psi[6];
	double h,r,psv[6],tmp,max;
	double k1,k2;
	pt[0].x=atof(input[2].item); 
	pt[0].y=atof(input[5].item);
	h=atof(input[4].item); 
	tmp=atof(input[3].item); 
	m=(tmp-pt[0].x)/h; m=((m<M)?m:M);
	max=pt[0].x+(double)m*h;
	tmp=(tmp<max)?tmp:max;
	pt[m].x=tmp;
	r=atof(input[6].item);
	psi[1]=23; psi[2]=24;
	for (i=0;i<=m;i++)
	{
		psv[1]=pt[i].x;
		psv[2]=pt[i].y;
		k1=h*parse(input[1].item,2,psv,psi);
		psv[1]=pt[i].x+r*h;
		psv[2]=pt[i].y+r*k1;
		k2=h*parse(input[1].item,2,psv,psi);
		if (i<m)
		{
			pt[i+1].y=pt[i].y+(1-1/(2*r))*k1+1/(2*r)*k2;
			pt[i+1].x=pt[i].x+h;
		}
	}
}

void CCode10::ODE1RK4()
{
	
	int i,psi[6];
	double h,psv[6],tmp,max;
	double k1,k2,k3,k4;
	pt[0].x=atof(input[2].item); 
	pt[0].y=atof(input[5].item);
	h=atof(input[4].item); 
	tmp=atof(input[3].item); 
	m=(tmp-pt[0].x)/h; m=((m<M)?m:M);
	max=pt[0].x+(double)m*h;
	tmp=(tmp<max)?tmp:max;
	pt[m].x=tmp;
	psi[1]=23; psi[2]=24;
	for (i=0;i<=m;i++)
	{
		psv[1]=pt[i].x;
		psv[2]=pt[i].y;
		k1=h*parse(input[1].item,2,psv,psi);

		psv[1]=pt[i].x+h/2;
		psv[2]=pt[i].y+k1/2;
		k2=h*parse(input[1].item,2,psv,psi);

		psv[1]=pt[i].x+h/2;
		psv[2]=pt[i].y+k2/2;
		k3=h*parse(input[1].item,2,psv,psi);

		psv[1]=pt[i].x+h;
		psv[2]=pt[i].y+k3;
		k4=h*parse(input[1].item,2,psv,psi);
		if (i<m)
		{
			pt[i+1].y=pt[i].y+(k1+2*k2+2*k3+k4)/6;
			pt[i+1].x=pt[i].x+h;
		}
	}
}

void CCode10::ShowTable()
{
	CString str;
	CPoint hTable=CPoint(20,310);
	CRect rcTable=CRect(hTable.x,hTable.y,hTable.x+280,hTable.y+290);
	table.DestroyWindow();
    table.Create(WS_VISIBLE | WS_CHILD | WS_DLGFRAME | LVS_REPORT
        | LVS_NOSORTHEADER,rcTable,this,idc++);
	table.InsertColumn(0,"i",LVCFMT_CENTER,25);
	table.InsertColumn(1,"x",LVCFMT_CENTER,70);
	table.InsertColumn(2,"y=f(x)",LVCFMT_CENTER,70);
	if (fMenu==4 || fMenu==5)
		table.InsertColumn(3,((fMenu==4)?"p":"z"),LVCFMT_CENTER,70);
	for (int i=0;i<=m;i++)
	{
		str.Format("%d",i); table.InsertItem(i,str,0);
		str.Format("%lf",pt[i].x); table.SetItemText(i,1,str);
		str.Format("%lf",pt[i].y); table.SetItemText(i,2,str);
		if (fMenu==4 && i>3)
		{
			str.Format("%lf",pt[i].p); 
			table.SetItemText(i,3,str);
		}
		if (fMenu==5)