pdeeul0.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:");i:=i+1;asc(unknown[i-1])!=0},{});
if(i==1,{abort("You must enter at least one unknown function.");},{});
unknown[i-1]:=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):={"+
GetUserInput("Enter the right hand side of the equation "+unknown[i]+"'(t,x)=...")+"};"+chr(13)+chr(10)+
unknown[i]+"_(x):={"+GetUserInput("Enter the initial function for the unknown function "+unknown[i]+" :")+"};"+chr(13)+chr(10)+
unknown[i]+"_L(t):={"+GetUserInput("Enter the left hand boundary condition for the function "+unknown[i]+" :")+"};"+chr(13)+chr(10)+
unknown[i]+"_R(t):={"+GetUserInput("Enter the right hand 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 interval :")+";"+chr(13)+chr(10)+
"x1:="+GetUserInput("Enter the right end of the interval :")+";"+chr(13)+chr(10)+
"n:="+GetUserInput("Enter the number of intervals the (x0,x1) interval will be cut :")+";"+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 interval.");},{});
if((valid(x1)!=-1)&&(valid(x1)!=-2),
{autoconfig("You must define the variable x1 as the right end of the 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;
t:=t0;i:=0;dx:=(x1-x0)/n;dx2:=dx*2.0;dx4:=dx*dx;
while({valid(unknown[i])==-4},
{
	execute(unknown[i]+ "_C[n+2]:=("+unknown[i]+ "_D[n+2]:=0);" );
	if((valid(execute(unknown[i]+"_(x0)"))==-1)||(valid(execute(unknown[i]+"_(x0)"))==-2),
	{
		j:=0;
		while({j<=n},{execute(unknown[i]+ "_C[j+1]:="+unknown[i]+ "_(x0+dx*j)*1.0;");j:=j+1;});
		execute(unknown[i]+"_C[0]:="+unknown[i]+"_C[1];"+unknown[i]+"_C[n+2]:="+unknown[i]+"_C[n+1];");
		execute(unknown[i]+"(t,x):={xx:=(x-x0)/dx+1;xy:=TRUNC(xx);xz:=xx-xy;" +unknown[i]+ "_C[xy]*(1.0-xz)+"+unknown[i]+"_C[xy+1]*xz};"+unknown[i]+ "'x(t,x):={("+unknown[i]+"_C[j+1]-"+unknown[i]+"_C[j-1])/dx2};"+unknown[i]+ "'xx(t,x):={("+unknown[i]+"_C[j+1]+"+unknown[i]+"_C[j-1]-2.0*"+unknown[i]+"_C[j])/dx4};");
	},
	{
autoconfig("The initial function for "+unknown[i]+" is not defined or gives a non-numeric value.");
	});
	if((valid(execute(unknown[i]+"_L(t0)"))!=-1)&&(valid(execute(unknown[i]+"_L(t0)"))!=-2),
	{
autoconfig("The left hand boundary function for "+unknown[i]+" is not defined or gives a non-numeric value.");
	},{});
	if((valid(execute(unknown[i]+"_R(t0)"))!=-1)&&(valid(execute(unknown[i]+"_R(t0)"))!=-2),
	{
autoconfig("The right hand 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)"))==-1)||(valid(execute(unknown[i]+"'(t0,x0)"))==-2),
	{
		computeall1:=computeall1+unknown[i]+"_C[0]:=cc*("+unknown[i]+"_C[1]:="+unknown[i]+"_L(t))+dd*"+unknown[i]+"_C[2]+ee*"+unknown[i]+"_C[3];"+unknown[i]+"_C[n+2]:=cc*("+unknown[i]+"_C[n+1]:="+unknown[i]+"_R(t))+dd*"+unknown[i]+"_C[n]+ee*"+unknown[i]+"_C[n-1];";
		computeall2:=computeall2+unknown[i]+"_D[0]:=cc*("+unknown[i]+"_D[1]:="+unknown[i]+"_L(t))+dd*"+unknown[i]+"_D[2]+ee*"+unknown[i]+"_D[3];"+unknown[i]+"_D[n+2]:=cc*("+unknown[i]+"_D[n+1]:="+unknown[i]+"_R(t))+dd*"+unknown[i]+"_D[n]+ee*"+unknown[i]+"_D[n-1];";
		computeall3:=computeall3+unknown[i]+"_D[j]:=(y_seg"+numtos(i)+"[j]:=dt*"+unknown[i]+"'(t,x))+"+unknown[i]+"_C[j];";
		computeall4:=computeall4+unknown[i]+"_C[j]:=(y_seg"+numtos(i)+"[j]+dt*"+unknown[i]+"'(t,x))/2.0+"+unknown[i]+"_C[j];";
		computeall5:=computeall5+unknown[i]+"(t,x):={xx:=(x-x0)/dx+1;xy:=TRUNC(xx);xz:=xx-xy;" +unknown[i]+ "_C[xy]*(1.0-xz)+"+unknown[i]+"_C[xy+1]*xz};"+unknown[i]+ "'x(t,x):={("+unknown[i]+"_C[j+1]-"+unknown[i]+"_C[j-1])/dx2};"+unknown[i]+ "'xx(t,x):={("+unknown[i]+"_C[j+1]+"+unknown[i]+"_C[j-1]-2.0*"+unknown[i]+"_C[j])/dx4};";
		computeall6:=computeall6+unknown[i]+"(t,x):={xx:=(x-x0)/dx+1;xy:=TRUNC(xx);xz:=xx-xy;" +unknown[i]+ "_D[xy]*(1.0-xz)+"+unknown[i]+"_D[xy+1]*xz};"+unknown[i]+ "'x(t,x):={("+unknown[i]+"_D[j+1]-"+unknown[i]+"_D[j-1])/dx2};"+unknown[i]+ "'xx(t,x):={("+unknown[i]+"_D[j+1]+"+unknown[i]+"_D[j-1]-2.0*"+unknown[i]+"_D[j])/dx4};";
		computeall7:=computeall7+"x_seg:="+unknown[i]+"_C[j];"+unknown[i]+"_C[j]:=-(y_seg"+numtos(i)+"+"+unknown[i]+"_C[j+2])/6.0+((y_seg"+numtos(i)+":=x_seg"+numtos(i)+")+"+unknown[i]+"_C[j+1])/1.5;x_seg"+numtos(i)+":=x_seg;";
		computeall8:=computeall8+"y_seg"+numtos(i)+":="+unknown[i]+"_C[0];x_seg"+numtos(i)+":="+unknown[i]+"_C[1];";
	},
	{
autoconfig("The equation for "+unknown[i]+" is not defined or gives a non-numeric value.");
	});
	i:=i+1;
});
j:=1;cc:=3;dd:=-3;ee:=1;
computeall1:="{"+computeall5+computeall1+
"j:=0;while({j<=n},{x:=x0+j*dx;j:=j+1;"+computeall3+"});t:=t+dt;"+
computeall6+computeall2+
"j:=0;while({j<=n},{x:=x0+j*dx;j:=j+1;"+computeall4+"});"+
computeall1+"updateviews(0);}";
computeall:=execute(computeall1);
execviewtemplates(0);
rr(t):={r:=0;o:=0;while({r<6.283},{o:=max(o,abs(u(t,r)-sin(r)));r:=r+dx;});o};
cont(mt):={maxt:=mt;while({t<=maxt},computeall);};
cont(maxt);