Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@samiyare Notice that I have added a faster version in my answer above.

You will have a better chance of getting a response if you provide the code. You could upload a worksheet.

@samiyare You could try minimizing the sum of squares of pC and pT using LSSolve from the Optimization package.

Input something like

Optimization:-LSSolve([pC,pT]);

It worked for me in Maple 12, which is what I have available this week.

@samiyare You could try minimizing the sum of squares of pC and pT using LSSolve from the Optimization package.

Input something like

Optimization:-LSSolve([pC,pT]);

It worked for me in Maple 12, which is what I have available this week.

@samiyare The algorithm in the link you give could perhaps be adjusted to the present situation like this:

solpar:=dsolve(subs(f(x)=F(x),{eq2=0,eq3=0, C(0)=1,D(C)(0)=c1, T(0)=1, D(T)(0) = t1}), parameters=[c1,t1],numeric,output=listprocedure,known=F);
Cn,Tn:=op(subs(solpar,[C(x),T(x)]));
#Making procedure pC and pT which compute C(b) and T(b), i.e. Cn(b) and Tn(b)
pC:=proc(c1,t1) if not type([c1,t1],list(numeric)) then return 'procname(_passed)' end if;
solpar(parameters=[c1,t1]);
Cn(b)
end proc;

pT:=proc(c1,t1) if not type([c1,t1],list(numeric)) then return 'procname(_passed)' end if;
solpar(parameters=[c1,t1]);
Tn(b)
end proc;

#Testing on previously found
sol(0);
pC(-2.11068134505884,-.922230659910390);
pT(-2.11068134505884,-.922230659910390);

Then you would have to use fsolve on the system {pC(c1,t1),pT(c1,t1)}.



@samiyare The algorithm in the link you give could perhaps be adjusted to the present situation like this:

solpar:=dsolve(subs(f(x)=F(x),{eq2=0,eq3=0, C(0)=1,D(C)(0)=c1, T(0)=1, D(T)(0) = t1}), parameters=[c1,t1],numeric,output=listprocedure,known=F);
Cn,Tn:=op(subs(solpar,[C(x),T(x)]));
#Making procedure pC and pT which compute C(b) and T(b), i.e. Cn(b) and Tn(b)
pC:=proc(c1,t1) if not type([c1,t1],list(numeric)) then return 'procname(_passed)' end if;
solpar(parameters=[c1,t1]);
Cn(b)
end proc;

pT:=proc(c1,t1) if not type([c1,t1],list(numeric)) then return 'procname(_passed)' end if;
solpar(parameters=[c1,t1]);
Tn(b)
end proc;

#Testing on previously found
sol(0);
pC(-2.11068134505884,-.922230659910390);
pT(-2.11068134505884,-.922230659910390);

Then you would have to use fsolve on the system {pC(c1,t1),pT(c1,t1)}.



@AliKhan In the new version you changed the definition of C3. Unfortunately with the consequence that the formal variable n is not seen after the arrow. Since the body of the procedure is not evaluated at the time of definition this means that something like C3(k) or C3(4) does not result in the expected. The solution is simple, either keep the original definition or use unapply:

C3 := unapply(C1+C2, n):

Btotal can be defined like this:

Btotal:=add(add(Bar[m,n],n=1..5),m=1..5):

@AliKhan In the new version you changed the definition of C3. Unfortunately with the consequence that the formal variable n is not seen after the arrow. Since the body of the procedure is not evaluated at the time of definition this means that something like C3(k) or C3(4) does not result in the expected. The solution is simple, either keep the original definition or use unapply:

C3 := unapply(C1+C2, n):

Btotal can be defined like this:

Btotal:=add(add(Bar[m,n],n=1..5),m=1..5):

@Markiyan Hirnyk

limit(J,a=R/sqrt(2),left) assuming R>0;

Answer: 5/3*sqrt(2)*Pi*R^3-4/3*R^3*Pi

(Now done as it happens in Maple 16).

@Markiyan Hirnyk

limit(J,a=R/sqrt(2),left) assuming R>0;

Answer: 5/3*sqrt(2)*Pi*R^3-4/3*R^3*Pi

(Now done as it happens in Maple 16).

@Christopher2222 Here it is:

restart;
L := .22:
lambda := 3.4:
rho := 2950:
cp := 1050:
a := lambda/(rho*cp):
h := 9:
T0 := 298:
sigma := 5.6697*10^(-8): epsilon:= .85:
PDE := diff(T(x, t), t) = lambda*(diff(T(x, t), x, x))/(rho*cp);
IBC := {-lambda*(D[1](T))(L, t) = h*(T(L, t)-T0)+sigma*epsilon*(T(L, t)^4-T0^4), T(0, t) = 1573, T(x, 0) = 298};
pde:=diff(u(x,t),t)=(1+u(x,t)^2)/(1+t^2);
ibc:={u(x,0)=0};
SOL:=pdsolve({PDE,pde}, IBC union ibc, numeric,timestep=100);

#SOL := pdsolve(PDE, IBC, numeric, timestep = 100);
SOL:-plot(T,t = 80*3600, thickness = 3);
SOL:-animate(T,t = 0 .. 80*3600, frames = 200, thickness = 3);
p1 := SOL:-value(T);
p1(.22, 80*3600);

# In the uploaded worksheet I have also given another method: Add the term T(x,t)^2*piecewise(x<0,1,0) to the right hand side of the PDE to make the equation nonlinear. But since x>0 in our problem this doesn't hurt.

pde_workaround.mw

@Christopher2222 Here it is:

restart;
L := .22:
lambda := 3.4:
rho := 2950:
cp := 1050:
a := lambda/(rho*cp):
h := 9:
T0 := 298:
sigma := 5.6697*10^(-8): epsilon:= .85:
PDE := diff(T(x, t), t) = lambda*(diff(T(x, t), x, x))/(rho*cp);
IBC := {-lambda*(D[1](T))(L, t) = h*(T(L, t)-T0)+sigma*epsilon*(T(L, t)^4-T0^4), T(0, t) = 1573, T(x, 0) = 298};
pde:=diff(u(x,t),t)=(1+u(x,t)^2)/(1+t^2);
ibc:={u(x,0)=0};
SOL:=pdsolve({PDE,pde}, IBC union ibc, numeric,timestep=100);

#SOL := pdsolve(PDE, IBC, numeric, timestep = 100);
SOL:-plot(T,t = 80*3600, thickness = 3);
SOL:-animate(T,t = 0 .. 80*3600, frames = 200, thickness = 3);
p1 := SOL:-value(T);
p1(.22, 80*3600);

# In the uploaded worksheet I have also given another method: Add the term T(x,t)^2*piecewise(x<0,1,0) to the right hand side of the PDE to make the equation nonlinear. But since x>0 in our problem this doesn't hurt.

pde_workaround.mw

In the procedure `pdsolve/default/construct_method` (lines  176 - 181 in Maple 12) there is a check on the value of a local table INFO produced by `pdsolve/numeric` and given as argument to `pdsolve/default/construct_method`.
If INFO["linear"] = true then .... else .... end if;
The problem seems to be that somehow INFO["linear"] has got the value 'true' (in Maple 12) at that point, whereas at the similar point in Maple 16 it has the value 'false'.

If you introduce a peaceful nonlinear pde quite independent of the other pde it appears that you are able to work your way around the problem

restart;
PDE := diff(T(x, t), t) = diff(T(x, t), x, x);
IBC_4:={D[1](T)(1, t) = -T(1,t)^4, T(0, t) = 1, T(x, 0) = 2};
pde:=diff(u(x,t),t)=(1+u(x,t)^2)/(1+t^2);
pdsolve(pde);
ibc:={u(x,0)=0};
SOL:=pdsolve({PDE,pde}, IBC_4 union ibc, numeric,spacestep=.005);
SOL:-plot(T,t=2);
SOL:-plot(u,t=2);
SOL:-value(T)(1,1);

#This workaround also works on the original problem from 2010:

http://www.mapleprimes.com/questions/97391-Whats-Wrong-With-Maple-Solution

PS. Instead you could add the term T(x,t)^2*piecewise(x<0,1,0) to the right hand side of the PDE to make the equation nonlinear. Since x>0 in our problem this doesn't hurt.

@Alejandro Jakubi Here is one:

restart;
PDE := diff(T(x, t), t) = diff(T(x, t), x, x);
IBC_4:={D[1](T)(1, t) = -T(1,t)^4, T(0, t) = 1, T(x, 0) = 2};
IBC_0:={D[1](T)(1, t) = 0, T(0, t) = 1, T(x, 0) = 2};
IBC_1:={D[1](T)(1, t) = -T(1,t), T(0, t) = 1, T(x, 0) = 2};
SOL_4 := pdsolve(PDE, IBC_4, numeric,spacestep=.005);
SOL_0 := pdsolve(PDE, IBC_0, numeric,spacestep=.005);
SOL_1 := pdsolve(PDE, IBC_1, numeric,spacestep=.005);
p4:=SOL_4:-plot(t = 1, color = red):
p0:=SOL_0:-plot(t = 1, color = blue):
p1:=SOL_1:-plot(t = 1, color = black):
plots:-display(p0,p4,p1);
SOL_4:-value()(1,1);
SOL_0:-value()(1,1);
SOL_1:-value()(1,1);

#There is no difference between SOL_4 and SOL_0 in Maple 12. However, SOL_1 is clearly different from the other two.   In Maple 16 all 3 are different and SOL_0 and SOL_1 are as in Maple 12.
This issue came up in 2010:

http://www.mapleprimes.com/questions/97391-Whats-Wrong-With-Maple-Solution

What I said was that there was no convergence in the original problem after reformulation neither in Maple 12 nor in Maple 16.

The original code stripped to essentials:

restart;
L := .22:
lambda := 3.4:
rho := 2950:
cp := 1050:
a := lambda/(rho*cp):
h := 9:
T0 := 298:
sigma := 5.6697*10^(-8): epsilon:= .85:
PDE := diff(T(x, t), t) = lambda*(diff(T(x, t), x, x))/(rho*cp);
IBC := {-lambda*(D[1](T))(L, t) = h*(T(L, t)-T0)+sigma*epsilon*(T(L, t)^4-T0^4), T(0, t) = 1573, T(x, 0) = 298};
SOL := pdsolve(PDE, IBC, numeric, timestep = 100);
SOL:-plot(t = 80*3600, thickness = 3);
SOL:-animate(t = 0 .. 80*3600, frames = 200, thickness = 3);
p1 := SOL:-value();
p1(.22, 80*3600);
##That "worked", but the result is wrong in Maple 12, but fine in Maple 16.
####################################################
##And the reformulated version:
restart;
L := .22:
lambda := 3.4:
rho := 2950:
cp := 1050:
a := lambda/(rho*cp):
h := 9:
T0 := 298:
sigma := 5.6697*10^(-8): epsilon:= .85:
sys:= {diff(T(x, t), t) = lambda*(diff(T(x, t), x, x))/(rho*cp),T4(x,t)=T(x,t)^4};
IBC := {-lambda*(D[1](T))(L, t) = h*(T(L, t)-T0)+sigma*epsilon*(T4(L, t)-T0^4), T(0, t) = 1573, T(x, 0) = 298};
SOL := pdsolve(sys, IBC, numeric, timestep = 100);
SOL:-plot(t = 80*3600, thickness = 3);
##Here there is an error message in Maple 12 and 16:

Error, (in pdsolve/numeric/plot) unable to compute solution for t>HFloat(0.0):
Newton iteration is not converging
#The only difference between Maple  12 and 16 is that Maple 12 says 't>0.' instead of 't>HFloat(0.0)'.



First 198 199 200 201 202 203 204 Last Page 200 of 231