npdeeu0.prg

defequ

Execute(GetConfig(1));
autoconfig(x__):={if(AskUserYesNo("There is an error in your main input configuration:"+chr(13)+chr(10)+x__+chr(13)+chr(10)+"Do you want to configure the main input automatically?")==0,{abort(x__);},
{
i:=0;
while({unknown[i]:=GetUserInput("Enter the "+numtos(i+1)+"th unknown function:");asc(unknown[i])!=0},{i:=i+1;});
if(i==0,{abort("You must enter at least one unknown function.");},{});
unknown[i]:=0;
i:=0;computeall:="";
while({valid(unknown[i])==-4},{
computeall:=computeall+"unknown["+numtos(i)+"]:="+chr(34)+unknown[i]+chr(34)+";"+chr(13)+chr(10)+unknown[i]+"'(t,x,y):={"+
GetUserInput("Enter the right hand side of the equation "+unknown[i]+"'(t,x,y)=...")+"};"+chr(13)+chr(10)+
unknown[i]+"_(x,y):={"+GetUserInput("Enter the initial function for the unknown function "+unknown[i]+" :")+"};"+chr(13)+chr(10)+
unknown[i]+"_x0(t,y):={"+GetUserInput("Enter the x=x0 boundary condition for the function "+unknown[i]+" :")+"};"+chr(13)+chr(10)+
unknown[i]+"_x1(t,y):={"+GetUserInput("Enter the x=x1 boundary condition for the function "+unknown[i]+" :")+"};"+chr(13)+chr(10)+
unknown[i]+"_y0(t,x):={"+GetUserInput("Enter the y=y0 boundary condition for the function "+unknown[i]+" :")+"};"+chr(13)+chr(10)+
unknown[i]+"_y1(t,x):={"+GetUserInput("Enter the y=y1 boundary condition for the function "+unknown[i]+" :")+"};"+chr(13)+chr(10);
i:=i+1;
});
computeall:=computeall+"t0:="+GetUserInput("Enter the initial time:")+";"+
chr(13)+chr(10)+"maxt:="+GetUserInput("Enter the time extent:")+";"+chr(13)+chr(10)+
"dt:="+GetUserInput("Enter the time step:")+";"+chr(13)+chr(10)+
"x0:="+GetUserInput("Enter the left end of the x interval :")+";"+chr(13)+chr(10)+
"x1:="+GetUserInput("Enter the right end of the x interval :")+";"+chr(13)+chr(10)+
"y0:="+GetUserInput("Enter the left end of the y interval :")+";"+chr(13)+chr(10)+
"y1:="+GetUserInput("Enter the right end of the y interval :")+";"+chr(13)+chr(10)+
"n:="+GetUserInput("Enter the number of grids :")+";"+chr(13)+chr(10);
SetConfig(computeall);
Abort("The setup is complete."+chr(13)+chr(10)+"Use File/Execute to start again.");
});
};
if((valid(t0)!=-1)&&(valid(t0)!=-2),
{autoconfig("You must define the variable t0 as the initial time.");},{});
if((valid(dt)!=-1)&&(valid(dt)!=-2),
{autoconfig("You must define the variable dt as the step size.");},{});
if((valid(maxt)!=-1)&&(valid(maxt)!=-2),
{autoconfig("You must define the variable maxt as the maximum time.");},{});
if((valid(x0)!=-1)&&(valid(x0)!=-2),
{autoconfig("You must define the variable x0 as the left end of the x interval.");},{});
if((valid(x1)!=-1)&&(valid(x1)!=-2),
{autoconfig("You must define the variable x1 as the right end of the x interval.");},{});
if((valid(y0)!=-1)&&(valid(y0)!=-2),
{autoconfig("You must define the variable y0 as the left end of the y interval.");},{});
if((valid(y1)!=-1)&&(valid(y1)!=-2),
{autoconfig("You must define the variable y1 as the right end of the y interval.");},{});
if(valid(n)!=-1,
{autoconfig("You must define the variable n as the number of segments the interval will be divided.");},{});
if((valid(unknown[0])!=-4),
{autoconfig("You must define the vector unknown[ ] as the names of the unknown functions.");},{});
t0:=t0*1.0;dt:=dt*1.0;maxt:=maxt*1.0;
x1:=x1*1.0;x0:=x0*1.0;y1:=y1*1.0;y0:=y0*1.0;
t:=t0;i:=0;
n2:=n+3;n3:=n2*n2;
dx:=(x1-x0)/n;dy:=(y1-y0)/n;
dx2:=dx*2.0;dy2:=dy*2.0;
dx4:=dx*dx;dy4:=dy*dy;
dx8:=4*dx4;dy8:=4*dy4;
while({valid(unknown[i])==-4},
{
	execute(unknown[i]+ "_C[n3]:=("+unknown[i]+ "_D[n3]:=0);" );
	if((valid(execute(unknown[i]+"_(x0,y0)"))==-1)||(valid(execute(unknown[i]+"_(x0,y0)"))==-2),
	{
		j:=0;
		while({j<=n},{k:=0;while({k<=n},{
		execute(unknown[i]+ "_C[(j+1)*n2+k+1]:="+unknown[i]+ "_(x0+dx*j,y0+k*dy)*1.0;");
		k:=k+1;});
		execute(unknown[i]+"_C[(j+1)*n2]:="+unknown[i]+"_C[(j+1)*n2+1];");
		execute(unknown[i]+"_C[(j+2)*n2-1]:="+unknown[i]+"_C[(j+2)*n2-2];");
		j:=j+1;});
		k:=0;while({k<=n},{
		execute(unknown[i]+"_C[k+1]:="+unknown[i]+"_C[n2+k+1];");
		execute(unknown[i]+"_C[(n2-1)*n2+k+1]:="+unknown[i]+"_C[(n2-2)*n2+k+1];");
		k:=k+1;});
		execute(unknown[i]+"_C[0]:="+unknown[i]+"_C[1];"+unknown[i]+"_C[n+2]:="+unknown[i]+"_C[n+1];");
		execute(unknown[i]+"_C[(n2-1)*n2]:="+unknown[i]+"_C[(n2-1)*n2+1];"+unknown[i]+"_C[n2*n2-1]:="+unknown[i]+"_C[n2*n2-2];");
		execute(unknown[i]+"(t,x,y):={xx:=(x-x0)/dx+1;xy:=TRUNC(xx);xz:=xx-xy;yx:=(y-y0)/dy+1;yy:=TRUNC(yx);yz:=yx-yy;"+
		unknown[i]+ "_C[xy*n2+yy]*(1.0-xz)*(1.0-yz)+"+unknown[i]+"_C[xy*n2+yy+1]*(1.0-xz)*yz+"+
		unknown[i]+ "_C[(xy+1)*n2+yy]*xz*(1.0-yz)+"+unknown[i]+"_C[(xy+1)*n2+yy+1]*xz*yz};"+
		unknown[i]+ "'x(t,x,y):={("+unknown[i]+"_C[(j+1)*n2+k]-"+unknown[i]+"_C[(j-1)*n2+k])/dx2};"+
		unknown[i]+ "'y(t,x,y):={("+unknown[i]+"_C[j*n2+k+1]-"+unknown[i]+"_C[j*n2+k-1])/dx2};"+
		unknown[i]+ "'xx(t,x,y):={("+unknown[i]+"_C[(j+1)*n2+k]+"+unknown[i]+"_C[(j-1)*n2+k]-2.0*"+unknown[i]+"_C[j*n2+k])/dx4};"+
		unknown[i]+ "'yy(t,x,y):={("+unknown[i]+"_C[j*n2+k+1]+"+unknown[i]+"_C[j*n2+k-1]-2.0*"+unknown[i]+"_C[j*n2+k])/dx4};"+
		unknown[i]+ "'xy(t,x,y):={("+unknown[i]+"_C[(j+1)*n2+k+1]-"+unknown[i]+"_C[(j+1)*n2+k-1]-"+unknown[i]+"_C[(j-1)*n2+k+1]+"+unknown[i]+"_C[(j-1)*n2+k-1])/dx8};");
	},
	{
autoconfig("The initial function for "+unknown[i]+" is not defined or gives a non-numeric value.");
	});
	if((valid(execute(unknown[i]+"_x0(t0,y0)"))!=-1)&&(valid(execute(unknown[i]+"_x0(t0,y0)"))!=-2),
	{
autoconfig("The x=x0 boundary function for "+unknown[i]+" is not defined or gives a non-numeric value.");
	},{});
	if((valid(execute(unknown[i]+"_x1(t0,y0)"))!=-1)&&(valid(execute(unknown[i]+"_x1(t0,y0)"))!=-2),
	{
autoconfig("The x=x1 boundary function for "+unknown[i]+" is not defined or gives a non-numeric value.");
	},{});
	if((valid(execute(unknown[i]+"_y0(t0,x0)"))!=-1)&&(valid(execute(unknown[i]+"_y0(t0,x0)"))!=-2),
	{
autoconfig("The y=y0 boundary function for "+unknown[i]+" is not defined or gives a non-numeric value.");
	},{});
	if((valid(execute(unknown[i]+"_y1(t0,x0)"))!=-1)&&(valid(execute(unknown[i]+"_y1(t0,x0)"))!=-2),
	{
autoconfig("The y=y1 boundary function for "+unknown[i]+" is not defined or gives a non-numeric value.");
	},{});
	i:=i+1;
});
computeall1:=(computeall2:=(computeall3:=(computeall4:=(computeall5:=(computeall6:=(computeall7:=(computeall8:="")))))));
i:=0;j:=1;
while({valid(unknown[i])==-4},
{
	if((valid(execute(unknown[i]+"'(t0,x0,y0)"))==-1)||(valid(execute(unknown[i]+"'(t0,x0,y0)"))==-2),
	{
		computeall1:=computeall1+"j:=0;while({j<=n},{"+
		unknown[i]+"_C[(j+1)*n2]:=("+unknown[i]+"_C[(j+1)*n2+1]:="+unknown[i]+"_y0(t,x));"+
		unknown[i]+"_C[(j+2)*n2-1]:=("+unknown[i]+"_C[(j+2)*n2-2]:="+unknown[i]+"_y1(t,x));"+
		"j:=j+1;});k:=0;while({k<=n},{"+
		unknown[i]+"_C[k+1]:=("+unknown[i]+"_C[n2+k+1]:="+unknown[i]+"_x0(t,y));"+
		unknown[i]+"_C[(n2-1)*n2+k+1]:=("+unknown[i]+"_C[(n2-2)*n2+k+1]:="+unknown[i]+"_x1(t,y));"+
		"k:=k+1;});"+
		unknown[i]+"_C[0]:="+unknown[i]+"_C[1];"+unknown[i]+"_C[n+2]:="+unknown[i]+"_C[n+1];"+
		unknown[i]+"_C[(n2-1)*n2]:="+unknown[i]+"_C[(n2-1)*n2+1];"+unknown[i]+"_C[n2*n2-1]:="+unknown[i]+"_C[n2*n2-2];";
		computeall2:=computeall2+"j:=0;while({j<=n},{"+
		unknown[i]+"_D[(j+1)*n2]:=("+unknown[i]+"_D[(j+1)*n2+1]:="+unknown[i]+"_y0(t,x));"+
		unknown[i]+"_D[(j+2)*n2-1]:=("+unknown[i]+"_D[(j+2)*n2-2]:="+unknown[i]+"_y1(t,x));"+
		"j:=j+1;});k:=0;while({k<=n},{"+
		unknown[i]+"_D[k+1]:=("+unknown[i]+"_D[n2+k+1]:="+unknown[i]+"_x0(t,y));"+
		unknown[i]+"_D[(n2-1)*n2+k+1]:=("+unknown[i]+"_D[(n2-2)*n2+k+1]:="+unknown[i]+"_x1(t,y));"+
		"k:=k+1;});"+
		unknown[i]+"_D[0]:="+unknown[i]+"_D[1];"+unknown[i]+"_D[n+2]:="+unknown[i]+"_D[n+1];"+
		unknown[i]+"_D[(n2-1)*n2]:="+unknown[i]+"_D[(n2-1)*n2+1];"+unknown[i]+"_D[n2*n2-1]:="+unknown[i]+"_D[n2*n2-2];";
		computeall3:=computeall3+unknown[i]+"_D[j*n2+k]:=(y_seg"+numtos(i)+"[j*n2+k]:=dt*"+unknown[i]+"'(t,x,y))+"+unknown[i]+"_C[j*n2+k];";
		computeall4:=computeall4+unknown[i]+"_C[j*n2+k]:=(y_seg"+numtos(i)+"[j*n2+k]+dt*"+unknown[i]+"'(t,x,y))/2.0+"+unknown[i]+"_C[j*n2+k];";
		computeall5:=computeall5+unknown[i]+"(t,x,y):={xx:=(x-x0)/dx+1;xy:=TRUNC(xx);xz:=xx-xy;yx:=(y-y0)/dy+1;yy:=TRUNC(yx);yz:=yx-yy;"+
		unknown[i]+ "_C[xy*n2+yy]*(1.0-xz)*(1.0-yz)+"+unknown[i]+"_C[xy*n2+yy+1]*(1.0-xz)*yz+"+
		unknown[i]+ "_C[(xy+1)*n2+yy]*xz*(1.0-yz)+"+unknown[i]+"_C[(xy+1)*n2+yy+1]*xz*yz};"+
		unknown[i]+ "'x(t,x,y):={("+unknown[i]+"_C[(j+1)*n2+k]-"+unknown[i]+"_C[(j-1)*n2+k])/dx2};"+
		unknown[i]+ "'y(t,x,y):={("+unknown[i]+"_C[j*n2+k+1]-"+unknown[i]+"_C[j*n2+k-1])/dx2};"+
		unknown[i]+ "'xx(t,x,y):={("+unknown[i]+"_C[(j+1)*n2+k]+"+unknown[i]+"_C[(j-1)*n2+k]-2.0*"+unknown[i]+"_C[j*n2+k])/dx4};"+
		unknown[i]+ "'yy(t,x,y):={("+unknown[i]+"_C[j*n2+k+1]+"+unknown[i]+"_C[j*n2+k-1]-2.0*"+unknown[i]+"_C[j*n2+k])/dx4};"+
		unknown[i]+ "'xy(t,x,y):={("+unknown[i]+"_C[(j+1)*n2+k+1]-"+unknown[i]+"_C[(j+1)*n2+k-1]-"+unknown[i]+"_C[(j-1)*n2+k+1]+"+unknown[i]+"_C[(j-1)*n2+k-1])/dx8};";
		computeall6:=computeall6+unknown[i]+"(t,x,y):={xx:=(x-x0)/dx+1;xy:=TRUNC(xx);xz:=xx-xy;yx:=(y-y0)/dy+1;yy:=TRUNC(yx);yz:=yx-yy;"+
		unknown[i]+ "_D[xy*n2+yy]*(1.0-xz)*(1.0-yz)+"+unknown[i]+"_D[xy*n2+yy+1]*(1.0-xz)*yz+"+
		unknown[i]+ "_D[(xy+1)*n2+yy]*xz*(1.0-yz)+"+unknown[i]+"_D[(xy+1)*n2+yy+1]*xz*yz};"+
		unknown[i]+ "'x(t,x,y):={("+unknown[i]+"_D[(j+1)*n2+k]-"+unknown[i]+"_D[(j-1)*n2+k])/dx2};"+
		unknown[i]+ "'y(t,x,y):={("+unknown[i]+"_D[j*n2+k+1]-"+unknown[i]+"_D[j*n2+k-1])/dx2};"+
		unknown[i]+ "'xx(t,x,y):={("+unknown[i]+"_D[(j+1)*n2+k]+"+unknown[i]+"_D[(j-1)*n2+k]-2.0*"+unknown[i]+"_D[j*n2+k])/dx4};"+
		unknown[i]+ "'yy(t,x,y):={("+unknown[i]+"_D[j*n2+k+1]+"+unknown[i]+"_D[j*n2+k-1]-2.0*"+unknown[i]+"_D[j*n2+k])/dx4};"+
		unknown[i]+ "'xy(t,x,y):={("+unknown[i]+"_D[(j+1)*n2+k+1]-"+unknown[i]+"_D[(j+1)*n2+k-1]-"+unknown[i]+"_D[(j-1)*n2+k+1]+"+unknown[i]+"_D[(j-1)*n2+k-1])/dx8};";
	},
	{
autoconfig("The equation for "+unknown[i]+" is not defined or gives a non-numeric value.");
	});
	i:=i+1;
});
j:=1;
computeall1:="{"+computeall5+computeall1+
"j:=0;while({j<=n},{x:=x0+j*dx;j:=j+1;k:=0;while({k<=n},{y:=y0+k*dy;k:=k+1;"+
computeall3+"});});t:=t+dt;"+
computeall6+computeall2+
"j:=0;while({j<=n},{x:=x0+j*dx;j:=j+1;k:=0;while({k<=n},{y:=y0+k*dy;k:=k+1;"+
computeall4+"});});"+
computeall1+"updateviews(0);}";
computeall:=execute(computeall1);
execviewtemplates(0);
cont(mt):={maxt:=mt;while({t<=maxt},computeall);};
cont(maxt);