Preben Alsholm

13728 Reputation

22 Badges

20 years, 250 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@STHence If X(0) is not known, but is just

X0:= <seq(x0[i],i=1..7)>;

then you can still compute E.X0, but then you are faced with having to solve Y = C.X w.r.t. X0 in some "best" way.

@STHence If X(0) is not known, but is just

X0:= <seq(x0[i],i=1..7)>;

then you can still compute E.X0, but then you are faced with having to solve Y = C.X w.r.t. X0 in some "best" way.

@STHence Somehow angular brackets often disappear when you paste them into the editor.

X was supposed to be
X:=<seq(x[i](t),i=1..7)>;

that is the same as

X:=Vector([seq(x[i](t),i=1..7)]);

which is safer with this editor, but not as convenient. I have corrected it above. I hope it stays OK.

@STHence Somehow angular brackets often disappear when you paste them into the editor.

X was supposed to be
X:=<seq(x[i](t),i=1..7)>;

that is the same as

X:=Vector([seq(x[i](t),i=1..7)]);

which is safer with this editor, but not as convenient. I have corrected it above. I hope it stays OK.

@STHence The LeastSquare solution minimizes the 2-norm of the residual A.X-Y. So in what sense do you expect a better solution than that?

X:=LinearAlgebra:-LeastSquares(A,Y);
A.X-Y:
LinearAlgebra:-Norm(%,2);
#Returns 0.1478777205
#LeastSquare does the same as the following (not codewise):
use LinearAlgebra in
   LinearSolve(Transpose(A).A,Transpose(A).Y)
end use;

@STHence The LeastSquare solution minimizes the 2-norm of the residual A.X-Y. So in what sense do you expect a better solution than that?

X:=LinearAlgebra:-LeastSquares(A,Y);
A.X-Y:
LinearAlgebra:-Norm(%,2);
#Returns 0.1478777205
#LeastSquare does the same as the following (not codewise):
use LinearAlgebra in
   LinearSolve(Transpose(A).A,Transpose(A).Y)
end use;

@acer You are quite right. The following splits a plot structure containing a single curve, so is less general than the evalindets idea in that respect. Also individual curve options (like thickness and color are gone):

restart;
u:=sqrt(x^2-5)+3;
plot(u,view=0..15,thickness=6); p:=%:
M:=plottools:-getdata(p)[-1];
S:=ListTools:-Split(hastype,convert(M,listlist),undefined):
plots:-display(plot~([S]),p);

You don't need the p in the last line except for keeping the view option.

#Here is an updated version that seems to work well:

R:=proc(p) local f;
   f:=proc(curve)
      evalindets(curve,Or(listlist,Matrix),m->ListTools:-Split(hastype,convert(m,listlist),undefined))
   end proc;
   evalindets(p,specfunc(anything,CURVES),f);
end proc;

R(p); #Returns the plot hoped for

@acer You are quite right. The following splits a plot structure containing a single curve, so is less general than the evalindets idea in that respect. Also individual curve options (like thickness and color are gone):

restart;
u:=sqrt(x^2-5)+3;
plot(u,view=0..15,thickness=6); p:=%:
M:=plottools:-getdata(p)[-1];
S:=ListTools:-Split(hastype,convert(M,listlist),undefined):
plots:-display(plot~([S]),p);

You don't need the p in the last line except for keeping the view option.

#Here is an updated version that seems to work well:

R:=proc(p) local f;
   f:=proc(curve)
      evalindets(curve,Or(listlist,Matrix),m->ListTools:-Split(hastype,convert(m,listlist),undefined))
   end proc;
   evalindets(p,specfunc(anything,CURVES),f);
end proc;

R(p); #Returns the plot hoped for

Here is another way of removing HFloat(undefined) (as unfortunately seems necessary in Maple 16.02):

restart;
u:=sqrt(75-3*x^2);
plot(u); p:=%:
evalindets(p,Matrix,m->remove(hastype,convert(m,listlist),undefined));

#Turning the data matrix into a listlist makes it easier to remove undefined elements.
#To catch cases where the data (or some of it) in fact is already a listlist as in
http://www.mapleprimes.com/questions/142416-Problems-With-DEplot
#we could do
evalindets(p,Or(Matrix,listlist),m->remove(hastype,convert(m,listlist),undefined));

Here is another way of removing HFloat(undefined) (as unfortunately seems necessary in Maple 16.02):

restart;
u:=sqrt(75-3*x^2);
plot(u); p:=%:
evalindets(p,Matrix,m->remove(hastype,convert(m,listlist),undefined));

#Turning the data matrix into a listlist makes it easier to remove undefined elements.
#To catch cases where the data (or some of it) in fact is already a listlist as in
http://www.mapleprimes.com/questions/142416-Problems-With-DEplot
#we could do
evalindets(p,Or(Matrix,listlist),m->remove(hastype,convert(m,listlist),undefined));

@Axel Vogt Amazingly simple! And thank you!

@Axel Vogt Amazingly simple! And thank you!

@sanaz Why not settle for a numerical solution?

de := diff(x(t), t, t)+30*x(t)+.1*x(t)^3 = 29*cos(t);
#Series solution to order 15 (highest term order 14):
resser14:=dsolve({de, x(0) = 0, D(x)(0) = 0},x(t),series,order=15);
pol14:=convert(rhs(resser14),polynom);
resnum:=dsolve({de, x(0) = 0, (D(x))(0) = 0},numeric);
#The numerical solution seems to me so much more interesting than a series solution:
plots:-odeplot(resnum,[t,x(t)],0..10);
#odeplot allows you to show both together, if you so desire:
plots:-odeplot(resnum,[[t,x(t)],[t,pol14]],0..1,view=0..5);

#You may also take a look at
resser16:=dsolve({de, x(0) = 0, D(x)(0) = 0},x(t),series,order=17);
pol16:=convert(rhs(resser16),polynom);
plots:-odeplot(resnum,[[t,x(t)],[t,pol14],[t,pol16]],0..1,view=-5..5,legend=["Numerical","pol14","pol16"]);


@sanaz Why not settle for a numerical solution?

de := diff(x(t), t, t)+30*x(t)+.1*x(t)^3 = 29*cos(t);
#Series solution to order 15 (highest term order 14):
resser14:=dsolve({de, x(0) = 0, D(x)(0) = 0},x(t),series,order=15);
pol14:=convert(rhs(resser14),polynom);
resnum:=dsolve({de, x(0) = 0, (D(x))(0) = 0},numeric);
#The numerical solution seems to me so much more interesting than a series solution:
plots:-odeplot(resnum,[t,x(t)],0..10);
#odeplot allows you to show both together, if you so desire:
plots:-odeplot(resnum,[[t,x(t)],[t,pol14]],0..1,view=0..5);

#You may also take a look at
resser16:=dsolve({de, x(0) = 0, D(x)(0) = 0},x(t),series,order=17);
pol16:=convert(rhs(resser16),polynom);
plots:-odeplot(resnum,[[t,x(t)],[t,pol14],[t,pol16]],0..1,view=-5..5,legend=["Numerical","pol14","pol16"]);


It would probably help if you gave us the differential equation from which you obtain f(x).

As you give us f(x) above you end with  + ...  indicating that there are more terms (infinitely many, I suspect).

First 183 184 185 186 187 188 189 Last Page 185 of 230