Preben Alsholm

13471 Reputation

22 Badges

20 years, 214 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

Starting from the odes and leaving out the data corresponding to initial conditions (they won't help)  we can do:
 

restart;
ode1 := diff(SS(t), t) = -k*SS(t)-SS(t)/B;
ode2 := diff(PP(t), t) = k*SS(t)-PP(t)/T1p;

init:=PP(0)=0,SS(0)=100000;
dsolve({ode1,ode2,init});
sol:=combine(expand(%));
PS:=subs(sol,[PP(t),SS(t)]);
T:=<2,4,6,8>:
S:=<86089.76983, 74114.4849, 63804.98946, 54929.56828>:
P:=<8932.499583, 15965.05975, 21410.35044, 25533.98716>:
###
Pfu,Sfu:=op(map(unapply,PS,t));
RP:=convert(P-Pfu~(T),list);
RS:=convert(S-Sfu~(T),list);
RPS:=[op(RP),op(RS)];
###
res:=Optimization:-LSSolve(RPS,k=0..1,B=0..1000,T1p=0..10);
SOL:=eval(sol,res[2]);
Psol,Ssol:=op(subs(SOL,[PP(t),SS(t)]));
plot(T,P,style=point,symbolsize=20); p1:=%:
plot(T,S,style=point,symbolsize=20); p2:=%:
plot(Psol,t=0..8,color=blue); p1s:=%:
plot(Ssol,t=0..8,color=blue); p2s:=%:
plots:-display(p1,p1s);
plots:-display(p2,p2s);

You have some syntax problems. I have corrected them below. Besides, X has 6 elements while Y has 5. I removed 0 from X.
By the assignment model:= 1.*10^5*exp(-t*k-*t/B)    I have assumed that you meant  model:= 1.*10^5*exp(-t*k-t/B).

Then I observe that exp(-t*k-t/B) has really only one parameter, namely  a=k+1/B.
Thus you can only expect unique results for the quantity k+1/B not for k and B individually.

## Maybe you wanted

model:= 1.*10^5*exp(-t*k-B/t);

which actually gives a good fit. ?

restart;
model:= 1.*10^5*exp(-t*k-t/B); #corrected
X:=<2,4,6,8,10>; #corrected
Y := <100000, 86089.76983, 74114.4849, 63804.98946, 54929.56828>;
res:=Statistics:-Fit(model,X,Y,t);
Statistics:-Fit(model,X,Y,t,output=[parametervalues]);
eval(model,op(%));

There is agreement between the last output and the first output from Fit.
### Is your mdel actually what I corrected it to be or what is it?
Plotting:
 

plot(X,Y,style=point, symbolsize=20); p1:=%:
plot(res,t=2..10,color=blue); p2:=%:
plots:-display(p1,p2);

 

I made a procedure way back when:
 

`convert/sn`:=proc(f) local SN;
SN:=t->if ilog10(t)=0 then t else 10^(-ilog10(t))*t*10^``(ilog10(t)) end if;
subsindets(f,float,SN)
end proc:
convert(evalf(1000*Pi),sn);
expand(%);

 

Express the right hand side of the last line with the help of piecewise:
p:=piecewise(m>=0, cos(m*theta), -sin(m*theta));

There are several calls to `evalf/AiryAi`, e.g. in line 19 of AiryAi:
evalf(('AiryAi')(n,z));
The unevaluation quotes are actually preventing recursion. They are also used in procedures for returning unevaluated. You also find those in AiryAi (e.g. in line 35).

To see the numerical handling look at:
showstat(`evalf/AiryAi`);

I'm just doing deq3.
 

restart;
deq3 := m*diff(v3(t), t) = -m*g-k*v3(t)*abs(v3(t));
sol3:=dsolve({deq3, v3(0) = v__0}, implicit) assuming v__0>0;
sol3p:=op(solve(sol3,{v3(t)})) assuming v3(t)>0;
sol3n:=op(solve(sol3,{v3(t)})) assuming v3(t)<0;
T:=solve(rhs(sol3p)=0,t);
solve(rhs(sol3n)=0,t); # The same
sol3A:=piecewise(t<T,rhs(sol3p),rhs(sol3n));
m := 0.258e-2; k := 3.1*10^(-4); g := 9.82; v__0 := 15;
sol3A;
plot(sol3A,t=0..10);

There are a couple of mismatches of parentheses in the definitions of the P's.
You need unevaluation quotes when plotting since your procedure assumes that t::numeric.
 

P[0] := (x, y) -> ((1/3)*y, (1/3)*x);

P[1] := (x, y) -> (-(1/3)*x+1/3, (1/3)*y+1/3); 
P[2] := (x, y) -> ((1/3)*x, (1/3)*y+2/3); 
P[3] := (x, y) -> ((1/3)*x+1/3, -(1/3)*y+1);
P[4] := (x, y) -> (2/3-(1/3)*y, 2/3-(1/3)*x); 
P[5] := (x, y) -> ( (1/3)*x+1/3, 1/3-(1/3)*y);
P[6] := (x, y) -> ((1/3)*x+2/3, (1/3)*y);
P[7] := (x, y) -> (-(1/3)*x+1, (1/3)*y+1/3);
P[8] := (x, y) -> ((1/3)*x+2/3, (1/3)*y+2/3);
###
peano := proc (t::numeric, depth::integer)
  local q, r; global P;
  if depth = 0 then return 0, 0 end if;
  q := floor(9*t); 
  r := 9*t-q;
  return P[q](peano(r, depth-1))
end proc;
## Test:
peano(0.45,5);
##
f:=(x,y)->x+y;
plot('f(peano(t,5))',t=0..1); #Notice unevaluation quotes

Actually, the following version works too, but I wouldn't recommend it:
(I have simplified the original version slightly):

peano := proc (t, depth::integer) 
local q, r,res; global P;
if depth = 0 then return 0, 0 end if;
q := floor(9*t); 
r := 9*t-q;
res:=peano(r, depth-1);
P[q](res[1],res[2])
end proc;
##Now try it with your f :
f:=(x,y)->x+y:
peano(0.45,5); #Test
peano(t,5); #Test
## Notice again a split to get two arguments: 
plot(f(peano(t,5)[1],peano(t,5)[2]),t=0..1);


################# A numeric version that can be used with fsolve:
 

peano := proc (t::numeric, depth::integer) option remember;
  local q, r; global P;
  if depth = 0 then return 0, 0 end if;
  q := floor(9*t); 
  r := 9*t-q;
  P[q](peano(r, depth-1))
end proc;
peano1:=proc(t,d) peano(_passed)[1] end proc;
peano2:=proc(t,d) peano(_passed)[2] end proc;

rcurry(f,x); #just to show the meaning of rcurry
rcurry(peano1,5); # and with peano1
fsolve(rcurry(peano1,5),0.1..1);
fsolve(rcurry(peano2,5),0.1..1);

 

Before resorting to other means by following the links Markiyan gave you, you should look into the possibility of combining the original system Dsys with the new ode to form one single system. That is simplest and best, when possible.
Thus if u(x) is the dependent variable in Dsys (or one of the dependent variables) and
ode:=diff(y(x),x)=y(x)-u(x);
then just form one system consisting of Dsys and ode and include all initial conditions in dsolve.

## Note: If you follow the links given by Markiyan you will see that the above approach won't work for the situation considered there. It would correspond to you having this new ode:
ode2:=diff(y(x),x)=y(x)-u(x/2);
differing from the first version in having u(x/2) instead of u(x).

Since it is winter, you might try freezing the vectors:

sys:={u+v = <-2,4> , u-v = <-3,6>};
sysF:=evalindets(sys,Vector,freeze);
solve(sysF,{u,v});
thaw(%);
%;

 

Here is a numerical solution if you should be interested. For completeness I bring the whole code:

restart;
Eq2 := diff(Theta(t), t, t) = -omega^2*sin(Theta(t));
ics := Theta(0) = 0, (D(Theta))(0) = Vmax;
SolNum:=dsolve({Eq2,ics},numeric,parameters=[omega,Vmax],output=listprocedure);
TH:=subs(SolNum,Theta(t));
SolNum(parameters=[1,0.5]);
with(plots): with(plottools):
p:=t->display(plot([[0,1],[sin(TH(t)),1-cos(TH(t))]],color=black),disk([sin(TH(t)),1-cos(TH(t))],.03,color=blue),scaling=constrained);
animate(p,[t],t=0..25,frames=500,size=[600,default]);


 

The linear version ought to stop after the first use of dsolve since you use the full ics. The rest doesn't make any sense. It made sense when you used only Theta(0)=0 as in the question you asked a couple of days ago.
For the nonlinear version just use the option implicit in dsolve.
The variable _a is simply an integration variable. As you know, the name of the integration variable in a definite integral is basically irrelevant: int( x^2, x=0..1) is the same integral as int( _a^2, _a=0..1). If you don't like _a being used, then use subs to change it to something else like tau. Just don't use the name t since it obscures things: You will then have t in two meanings, which only physicists can handle.

Eq1 := diff(Theta(t), t, t) = -omega^2*Theta(t);
ics := Theta(0) = 0, (D(Theta))(0) = Vmax;
SHM := dsolve({Eq1, ics});
restart;
Eq2 := diff(Theta(t), t, t) = -omega^2*sin(Theta(t));
ics := Theta(0) = 0, (D(Theta))(0) = Vmax;
Sol:=dsolve({Eq2,ics}, implicit);
## Changing the name of the integration variable if you like.
## Here in Sol[1]:
res:=subs(_a=tau,Sol[1]);

The last integral can be done by Maple. Continue with:

value(res) assuming Vmax>0;

Result: 2*InverseJacobiAM((1/2)*Theta(t), 2*omega/Vmax)/Vmax-t = 0

Use odeplot from the plots package.
Example:

ode:=diff(y(x),x)=y(x);
sol:=dsolve({ode,y(0)=1});
solNum:=dsolve({ode,y(0)=1},numeric);
plots:-odeplot(solNum,[x,y(x)-rhs(sol)],0..3);

 

The error message tells the story. F needs two arguments. It is given only one in

[ seq(seq(Vector([v[],xn]),v in F(n-1)),xn=0..(p-1)) ];

  If you change F(n-1) to F(p,n-1) something is produced. Whether that is what you wanted, you decide.

The syntax for giving the derivatives uses D, not diff. You can only impose 2 conditions on a second order ode. You cannot impose conditions on the nth derivative (or higher) for an nth order ode.
Thus the ics should be written:
ics:=Theta(0) = 0, D(Theta)(0) = Vmax;

Actually, both CAS'es miss (as expected) infinitely many solutions to the initial value problem eq with y(0)=1.
I illustrate how you can construct as many as you like here:

restart;
eq:=(diff(y(x),x))^2=4 * y(x);
ode1:=diff(y(x),x)=2*sqrt(y(x));
ode2:=diff(y(x),x)=-2*sqrt(y(x));;
sol1:=dsolve({ode1,y(0)=1});
plot(rhs(sol1),x=-2..1); 
odetest(sol1,[eq,y(0)=1]);
solp1:=y(x)=piecewise(x<-2,(x+2)^2,x<-1,0,rhs(sol1));
plot(rhs(solp1),x=-3..1,thickness=3);
odetest(solp1,[eq,y(0)=1]);
###
sol2:=dsolve({ode2,y(0)=1});
plot(rhs(sol2),x=-1..2); 
odetest(sol2,[eq,y(0)=1]);
solp2:=y(x)=piecewise(x<1,rhs(sol2),x<2,0,(x-2)^2);
plot(rhs(solp2),x=-1..3,thickness=3);
odetest(solp2,[eq,y(0)=1]);

The graph of the last one:

 

First 33 34 35 36 37 38 39 Last Page 35 of 158