code10.cpp

ordinary differential equations

		{
			str.Format("%lf",pt[i].z); 
			table.SetItemText(i,3,str);
		}
	}
}

void CCode10::OnButton()
{
	if (fStatus==1)
	{
		for (int i=1;i<=nInputItems;i++)
			input[i].ed.GetWindowText(input[i].item);
		switch(fMenu)
		{
			case 1:
				ODE1Taylor(); fStatus=2; break;
			case 2:
				ODE1RK2(); fStatus=2; break;
			case 3:
				ODE1RK4(); fStatus=2; break;
			case 4:
				ODE1AB(); fStatus=2; break;
			case 5:
				ODE1System(); fStatus=2; break;
			case 6:				
				ODE2toODE1System(); fStatus=2; break;
			case 7:				
				ODE2FD1(); fStatus=2; break;
			case 8:				
				ODE2FD2(); fStatus=2; break;
		}
	}
	if (fStatus=2)
	{
		ShowTable();			
		InvalidateRect(curve.rc);
	}
}

void CCode10::ODE1System()
{
	
	int i,psi[6];
	double h,psv[6],tmp,max;
	double k1,k2,k3,k4;
	double K1,K2,K3,K4;
	pt[0].x=atof(input[3].item); 
	pt[0].y=atof(input[4].item); 
	pt[0].z=atof(input[5].item);
	h=atof(input[6].item); 
	tmp=atof(input[7].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; psi[3]=25;
	for (i=0;i<=m;i++)
	{
		psv[1]=pt[i].x;
		psv[2]=pt[i].y;
		psv[3]=pt[i].z;
		k1=h*parse(input[1].item,3,psv,psi); 
		K1=h*parse(input[2].item,3,psv,psi);

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

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

		psv[1]=pt[i].x+h;
		psv[2]=pt[i].y+k3;
		psv[3]=pt[i].z+K3;
		k4=h*parse(input[1].item,3,psv,psi); 
		K4=h*parse(input[2].item,3,psv,psi);

		if (i<m)
		{
			pt[i+1].y=pt[i].y+(k1+2*k2+2*k3+k4)/6;
			pt[i+1].z=pt[i].z+(K1+2*K2+2*K3+K4)/6;
			pt[i+1].x=pt[i].x+h;
		}
	}
}

void CCode10::ODE1Taylor()
{
	double psv[6],tmp,max;
	int psi[6];
	int i;
	double h,u,v,w,z;

	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;
	psi[3]=20; psi[4]=21; //psi[5]=22;
	for (i=0;i<=m;i++)
	{
		psv[1]=pt[i].x; 
		psv[2]=pt[i].y; 
		u=parse(input[1].item,2,psv,psi); 
		psv[3]=u; 
		v=parse(input[6].item,3,psv,psi); 
		psv[4]=v; 
		w=parse(input[7].item,4,psv,psi);
		if (i<m)
		{
			pt[i+1].y=pt[i].y+h*u+pow(h,2)/2*v+pow(h,3)/6*w;
			pt[i+1].x=pt[i].x+h;
		}
	}
}

void CCode10::ODE1AB()
{
	
	int i,psi[6];
	double h,f0,f1,f2,f3,fp,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++)
	{
		if (i<3)
		{
			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);

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

				psv[1]=pt[i-2].x;
				psv[2]=pt[i-2].y;
				f1=parse(input[1].item,2,psv,psi);

				psv[1]=pt[i-1].x;
				psv[2]=pt[i-1].y;
				f2=parse(input[1].item,2,psv,psi);

				psv[1]=pt[i].x; 
				psv[2]=pt[i].y; 
				f3=parse(input[1].item,2,psv,psi);
				pt[i+1].p=pt[i].y+h/24*(-9*f0+37*f1-59*f2+55*f3);

				psv[1]=pt[i+1].x;
				psv[2]=pt[i+1].p;
				fp=parse(input[1].item,2,psv,psi);
				pt[i+1].y=pt[i].y+h/24*(f1-5*f2+19*f3+9*fp);
			}
		}
	}
}

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

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

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

		psv[1]=pt[i].x+h; 
		psv[2]=pt[i].y+k3;
		psv[3]=pt[i].z+K3;
		k4=h*(pt[i].z+K3); K4=h*parse(input[1].item,3,psv,psi);

		if (i<m)
		{
			pt[i+1].y=pt[i].y+(k1+2*k2+2*k3+k4)/6;
			pt[i+1].z=pt[i].z+(K1+2*K2+2*K3+K4)/6;
			pt[i+1].x=pt[i].x+h;
		}
	}
}

void CCode10::ODE2FD1()
{
	
	int i,j,psi[6];
	double psv[6],h,tmp;
	double **a,*b;
	double *p,*q,*r,*w;
	b=new double [M+1];
	p=new double [M+1];
	q=new double [M+1];
	r=new double [M+1];
	w=new double [M+1];
	a=new double *[M+1];
	h=atof(input[5].item);
	m=(int)(atof(input[7].item)-atof(input[6].item))/h; m=((m<M)?m:M);
	pt[0].x=atof(input[6].item); 
	pt[m].x=pt[0].x+(double)m*h;
	pt[0].y=atof(input[8].item); 
	pt[m].y=atof(input[9].item);
	psi[1]=23;
	for (i=1;i<=m-1;i++)
	{
		pt[i].x=pt[i-1].x+h;
		psv[1]=pt[i].x; 
		p[i]=parse(input[1].item,1,psv,psi);
		q[i]=parse(input[2].item,1,psv,psi);
		r[i]=parse(input[3].item,1,psv,psi);
		w[i]=parse(input[4].item,1,psv,psi);
	}
	for (i=0;i<=M;i++)
		a[i]=new double [M+1];
	for (i=1;i<=m-1;i++)
	{
		for (j=1;j<=m-1;j++)
			a[i][j]=0;
		b[i]=0;
	}
	b[1]=w[1]-(p[1]/(h*h)-q[1]/(2*h))*pt[0].y;
	b[m-1]=w[m-1]-(p[m-1]/(h*h)+q[m-1]/(2*h))*pt[m].y;
	for (i=1;i<=m-1;i++)
	{
		a[i][i]=-2*p[i]/(h*h)+r[i];
		if (i<m-1)
			a[i+1][i]=p[i+1]/(h*h)-q[i+1]/(2*h);
		if (i>1)
			a[i-1][i]=p[i-1]/(h*h)+q[i-1]/(2*h);
		if (i>1 && i<m-1)
			b[i]=w[i];
	}
	SolveSLE(a,b);
	for (i=0;i<=M;i++)
		delete a[i];
	delete p,q,r,w,a,b;
}

void CCode10::SolveSLE(double **a,double *b)
{
	int i,j,k,lo,hi;
	double m1,Sum;
	if (fMenu==7)
	{
		lo=1; hi=m-1;
	}
	if (fMenu==8)
	{
		lo=0; hi=m;
	}
	for (k=lo;k<=hi-1;k++)
		for (i=k+1;i<=hi;i++)
		{
			m1=a[i][k]/a[k][k];
			for (j=lo;j<=hi;j++)
				a[i][j] -= m1*a[k][j];
			b[i] -= m1*b[k];
		}
	for (i=hi;i>=lo;i--)
	{
    	Sum=0;
		pt[i].y=0;
		for (j=i;j<=hi;j++)
       		Sum += a[i][j]*pt[j].y;
		pt[i].y=(b[i]-Sum)/a[i][i];
	}
}

void CCode10::ODE2FD2()
{
	
	int i,j,psi[6];
	double psv[6],tmp;
	double h,alpha,beta;
	double **a,*b;
	double *p,*q,*r,*w;
	b=new double [M+1];
	p=new double [M+1];
	q=new double [M+1];
	r=new double [M+1];
	w=new double [M+1];
	a=new double *[M+1];
	for (i=0;i<=M;i++)
		a[i]=new double [M+1];
	h=atof(input[5].item);
	m=(int)(atof(input[7].item)-atof(input[6].item))/h; m=((m<M)?m:M);
	pt[0].x=atof(input[6].item); 
	pt[m].x=pt[0].x+(double)m*h;
	alpha=atof(input[8].item); beta=atof(input[9].item); 
	for (i=0;i<=m;i++)
	{
		if (i<m)
			pt[i+1].x=pt[i].x+h;
		psv[1]=pt[i].x; psi[1]=23;
		p[i]=parse(input[1].item,1,psv,psi);
		q[i]=parse(input[2].item,1,psv,psi);
		r[i]=parse(input[3].item,1,psv,psi);
		w[i]=parse(input[4].item,1,psv,psi);
	}
	for (i=0;i<=m;i++)
	{
		for (j=0;j<=m;j++)
			a[i][j]=0;
		b[i]=0;
	}
	for (i=0;i<=m;i++)
	{
		a[i][i]=-2*p[i]/(h*h)+r[i];
		if (i>0 && i<m)
			a[i][i+1]=p[i]/(h*h)+q[i]/(2*h);
		if (i<m-1)
			a[i+1][i]=p[i+1]/(h*h)-q[i+1]/(2*h);
		if (i>0 && i<m)
			b[i]=w[i];
	}
	a[0][1]=2*p[0]/m;
	a[m][m-1]=2*p[m]/(h*h);
	b[0]=w[0]+(p[0]/(h*h)-q[0]/(2*h))*2*h*alpha;
	b[m]=w[m]-(p[m]/(h*h)+q[m]/(2*h))*2*h*beta;
	SolveSLE(a,b);
	for (i=0;i<=M;i++)
		delete a[i];
	delete p,q,r,w,a,b;
}

void CCode10::Clear(CRect rc)
{
	CClientDC dc(this);
	CBrush whiteBrush(RGB(255,255,255));
	dc.FillRect(&rc,&whiteBrush);
}