sys := [(diff(c(t), t))*(diff(a(t), t))*(diff(b(t), t))+(diff(c(t), t))*t^2, (diff(c(t), t))*(diff(a(t), t))*t+(diff(c(t), t))*t*(diff(b(t), t)), (diff(c(t), t))*(diff(a(t), t))*(diff(b(t), t))+(diff(c(t), t))*(diff(a(t), t))*t+(diff(c(t), t))*t^2]
DEplot(sys, [a(t), b(t), c(t)], t = 0 .. 2, a = -15 .. 15, b = -15 .. 15, c = -15 .. 15, color = magnitude, title = `Stable Limit Cycles`, arrows = curve, dirfield = 800, axes = none);
odeplot(sys, [t, a(t), b(t), c(t)], -4 .. 4, color = orange);

mysolve.mw

Download mysolve.mw

with(RegularChains);
with(ConstructibleSetTools);
source1 := PolynomialRing([a, b, c, d]);
target1 := PolynomialRing([e1, e2, e3, e4]);
source1list := [eq2a, eq3a, eq4a, 0];
target1list := [eq2b, eq3b, eq4b, eq5b];
cs := PolynomialMapPreimage(target1list, source1list, source1, target1)

return
cs := constructible_set

then i type
Elements(constructible_set)

it return 

Elements(constructible_set)

 since source has less number of variables, in order to use this function

i add a dummy variable d to this and use 0 for the fourth equation

dsol := dsolve([eq2=diff(s(t),t), eq3=diff(s(t),t), eq4=diff(s(t),t)], [a(t),b(t),c(t)]);
eval(subs(t(t)=t,subs(_C1(t) = _C1, subs(_C2(t) = _C2, subs(_C3(t) = _C3, sol1)))));

after subs t(t) = t , it still  have t(t)

hello all

i have two equations

(1-pi*x/2)*(f/2/E)*yc-1/3=int(t^4/sqrt((t^2-(f/2/E)*yc)^2+ya^2),t=0..1);

x^2*yc*(1-pi*x/2)-1/3=int(t^4/sqrt((t^2-x^2*yc)^2+ya^2),t=0..1);

i want to plot ya and yc vs. x that x change beetwin -3..3 in some steps like 100 steps. and another members like E can be any number.

i think maby i shoud use loop! but i dont know how!

plz help me

Hi,

I want to plot two density functions of norm(1,1) and norm(4,1) in one figure.

But it appears the function DensityPlot can only plot one at a time.

The downloaded worksheet below displays 3 points on the unit sphere which define a solid angle with a triangular face. The sides of the solid angle's are red arcs on the surface of the sphere and red radii which outline the planar sides within the sphere.

Three questions:

1. Is there a way to make the surfaces of the solid triangle more apparent by filling them with color?

2. Is there a way to calculate the area of the face on the surface of the sphere?

3. Is there a way to calculate the volume of the solid triangle?

Download Mechanics;_irregular_solid_angle.mw

I am trying to solve system of partial differential equation, but i got some trouble. This is system of heat equations

Here are the equations:

restart;with(PDEtools):
U1 := diff_table(u1(x,t)):
pde[1] := (1/kappa1)*U1[t]=U1[x,x];

bc[1]:=u1(-L, t)=V;

U2 := diff_table(u2(x,t)):
pde[2] := (1/kappa2)*U2[t]=U2[x,x];

bc[2]:=u2(a, t)=0;

sys:=[pde[1],bc[1],pde[2],bc[2], u1(x=0,t)=u2(x=0,t), 

D[2](u1)(x,0)=D[2](u2)(x,0)

]

pdsolve(sys);

a and L is arbitrary constanta

anyone can help? The last boundary conditions I mean du1/dx=du2/dx at x=0

How can i solve this problem? 

So given a list [ [ [1,2], [2,1] ], [ [3,4], [4,3] ], [ [5] ] ],

I would like to produce a list 

[ [1,2,3,4,5], [1,2,4,3,5], [2,1,3,4,5], [2,1,4,3,5] ]

so it makes a combination of sublist of every lists in the given list

is there a way to do it ?

Thanks,

Hi all

Assume that we have a cost functional, namely J, which we want to be minimized but we have to set of constraint CC1 and CC2.

what should we do it?

The code which I have written is attached(Minimize J under CC1[i]=0 and CC2[i]=0)

 any suggestion or guide is praiseworthy.

Thanks in advance

ExNew.mws

Mahmood   Dadkhah

Ph.D Candidate

Applied Mathematics Department

Dear all,


I'm having trouble defining a Maple procedure containing for loops. I have one outer for loop which I terminate at the end and two separate for loops of which the first one get's terminated before the second one starts. Nethertheless Maple complains, that my outer loop doesn't get terminated and points to to a point before the second inner for loop starts. I based this procedure on a procedure I defined in another worksheet and the latter one is working just fine, although the structure ist the same in both procedures. As far as I can see, the structure is

for n from 0 to T do
for i from 1 to J do
....
end do;
for i from J-2 to 1 by -1 do
end do;

end do;

I have several if statements in the procedure, too, but as far as I can see, everything is terminated fine. Perhaps there is something else wrong with my procedure and I'm just getting a strange error? I would really appreciateit, if somebody could take a short look at the problem. Since the procedure takes up several pages I'm uploading my Maple work sheet. The relevant procedure is the one defined just after the red comment "#  actual code" which is about 60% through the work sheet.

Thank you all.

BD_collapse_2.mw.

Dear Collagues

I wrote a code to solve a system of ode numerically. When the equations have been solved, I want to find the value of MTE which is in integral form. I check and see that the value of MTE is calculated wrong. I dont know why?

I calculated the values of MTE twice (another variable is masst below res1). Two different answers!!

I would be most grateful if you could help me.

Thank you

Here is my code

restart:
EPSILONE:=1:
#Digits:=10:

a[mu1]:=39.11:
b[mu1]:=533.9:
a[k1]:=7.47:
b[k1]:=0:

rhop:=5180:
rhobf:=998.2:
#gama1:=0.2;
mu1[bf]:=9.93/10000:
k1[bf]:=0.597:

rhost(eta):=1-phi(eta)+phi(eta)*rhop/rhobf;
#mu:=unapply(mu1[bf]/(1-phi(eta))^2.5,eta);
mu:=unapply(mu1[bf]*(1+a[mu1]*phi(eta)+b[mu1]*phi(eta)^2),eta);
#mu1[bf]:=9.93/10000:
k:=unapply(k1[bf]*(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2),eta);
mu(eta)/mu1[bf];


#eq1:=diff((mu(eta)/mu1[bf])*diff(u(eta),eta),eta)+rhost(eta)-Ha^2*u(eta);
eq1:=(a[mu1]*(diff(phi(eta), eta))+2*b[mu1]*phi(eta)*(diff(phi(eta), eta)))*(diff(u(eta), eta))+(1+a[mu1]*phi(eta)+b[mu1]*phi(eta)^2)*(diff(u(eta), eta, eta))+rhost(eta)-Ha^2*u(eta);
eq2:=diff((k(eta)/k1[bf])*diff(T(eta),eta),eta);
eq3:=diff(phi(eta),eta)+phi(eta)/(N_bt*(1+gama1*T(eta))^2)*diff(T(eta),eta):
eq2:=subs(phi(0)=phi0,eq2);
eq3:=subs(phi(0)=phi0,eq3);



eval(dsolve({eq3,phi(1)=phiv},phi(eta)),(T)(1)=1);
Phi:=normal(combine(%));
Teq:=isolate(eval(eq2,Phi),diff(T(eta),eta));
ueq1:=eval(eq1,Phi)=0:
ueq2:=subs(Teq,ueq1):


Agama1:=<<0.1>,<0.2>,<0.3>>:
ANBT:=<<0.2>,<1>,<10>>:
Aphiv:=<<0.02>,<0.06>,<0.1>>:
lambda:=0;
Ha:=0;
#NBT:=0.5:
N_bt:=cc*NBT+(1-cc)*10;
#phiv:=0.02:
 
kratio:=k1[p]/k1[bf];




eq1;
eq2;
eq3;

res1 := dsolve(subs(NBT=1,gama1=0.3,phiv=0.06,{eq1,eq2,eq3,u(0)=lambda*D(u)(0),D(u)(1)=0,T(0)=0,phi(1)=phiv,T(1)=1}), numeric,method=bvp[midrich],output=listprocedure,continuation=cc):
G0,G1,G2:=op(subs(subs(res1),[phi(eta),u(eta),diff(T(eta),eta)])):

masst:=evalf(int((1-G0(eta)+G0(eta)*rhop/rhobf)*G1(eta), eta = 0..1));
heatt:=(1+a[k1]*G0(0)+b[k1]*G0(0)^2)*G2(0);



for iNBT from 1 to 3 do

for iphi from 1 to 3 do

for igamma from 1 to 3 do

res := dsolve(subs(NBT=ANBT(iNBT),gama1=Agama1(igamma),phiv=Aphiv(iphi),{eq1,eq2,eq3,u(0)=lambda*D(u)(0),D(u)(1)=0,T(0)=0,phi(1)=phiv,T(1)=1}), numeric,method=bvp[midrich],output=listprocedure,continuation=cc):
F0,F1,F2:=op(subs(subs(res),[phi(eta),u(eta),diff(T(eta),eta)])):
MTE[iNBT,iphi,igamma]:=evalf(Int((1-F0(eta)+F0(eta)*rhop/rhobf)*F1(eta), eta = 0..1)):
HTC[iNBT,iphi,igamma]:=(1+a[k1]*F0(0)+b[k1]*F0(0)^2)*F2(0):

end do;

end do;

end do;

for iNBT from 1 to 3 do

for iphi from 1 to 3 do

for igamma from 1 to 3 do



print (convert(HTC[iNBT,iphi,igamma],float,9));

end do;

end do;

end do;
for iNBT from 1 to 3 do

for iphi from 1 to 3 do

for igamma from 1 to 3 do


print (convert(MTE[iNBT,iphi,igamma],float,9));

end do;

end do;

end do;

Amir

I want to add text to the left of the vertical axis so I need more background space such that the full text fits to the left

As it is now it will appear truncated. I do not want to change the graph itself but I want to increase the background area around the axes so I can enter a text block using the drawing tools. 

_C1, _C2, _C3 are constant, how to set them constant, to make diff(_C2) = 0 etc

eval(simplify(subs(a=_C1,subs(b=1/(diff(c(t), t)),subs(c=_C2+_C3*exp(-t),eq2)))));

(diff(_C2(t), t)+(diff(_C3(t), t))*(exp(-t))(t)+_C3(t)*(diff((exp(-t))(t), t)))*(-(diff(_C1(t), t))*(diff(diff((c(t))(t), t), t))/(diff((c(t))(t), t))^2+(diff(s(t), t))/(diff((c(t))(t), t)))
a