Are you talking about applying fsolve at each step of the integration, or are you talking about using the results of a prior integration, then fsolve, to start a new integration?
You could try something like the following:
> # Build a periodic quadratic over 0..1.2 with a period of 1.2
> forcefun := proc(x) local usex;
> if not type(x,'numeric') then
> 'procname'(x);
> else
> usex := 1.2*frac(x/1.2);
> 0.36-(usex-0.6)^2;
> end if;
> end proc:
> plot(forcefun,0..2.4);
+ AAAAAAA AAAAAA
0.35 AA AA AA AA
+ AA AA AA AA
0.3 A AA AA AA
+ A A A A
+ A A AA AA
0.25 A AA A A
+ AA AA AA AA
+ A A A A
0.2 A A A A
+ A A AA A
+ A A A A
0.15 A A A A
+ AA AA AA AA
0.1 A A A A
+ AA AA AA AA
+ A A A A
0.05A A A A
+A A A A
* AAA A
*--+--+--+--+--+--+---+--+--+--+--+--*--+--+--+--+--+---+--+--+--+--+--+--*
0 0.5 1 1.5 2
> # Now use if to solve a DE:
> dsn := dsolve({diff(y(x),x)=forcefun(x),y(0)=0},numeric,known={forcefun});
dsn := proc(x_rkf45) ... end proc
> plots[odeplot](dsn,0..2.4);
+ AAAAAA
+ AAAAA
+ AAA
0.5+ AA
+ AAA
+ AAA
0.4+ AA
+ AAA
+ AAA
+ AAAA
0.3+ AAAAAAAAA
+ AAAA
+ AAA
0.2+ AAA
+ AA
+ AAA
+ AAA
0.1+ AA
+ AAA
+ AAAAA
-******+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+-
0.5 1 1.5 2
If you have particular values for 'sigma' and 'b' that you are interested in, it may be possible to obtain a numerical answer as follows:
> xpr := (exp(-1/2*(ln(y)-(x*b+2))^2/(si\
> gma^2*(x^2+1)))/(y*sigma*sqrt(2*Pi*(x^2+1))))*(2-5/3*x);
2
(ln(y) - x b - 2) 1/2 / 5 x\
exp(-1/2 ------------------) 2 |2 - ---|
2 2 \ 3 /
sigma (x + 1)
xpr := 1/2 -------------------------------------------
2 1/2
y sigma (Pi (x + 1))
> ii := Int(xpr,[x=0.6..1.2, y=3.6..4.5]);
2
(ln(y) - x b - 2) 1/2 / 5 x\
/ exp(-1/2 ------------------) 2 |2 - ---|
| 2 2 \ 3 /
| sigma (x + 1)
ii := | 1/2 ------------------------------------------- d
| 2 1/2
| y sigma (Pi (x + 1))
|
/
[x = 0.6 .. 1.2, y = 3.6 .. 4.5]
> evalf(eval(ii,{b=1,sigma=1}));
0.01152171054
> evalf(eval(ii,{b=2,sigma=1}));
0.004940803670
I will assume that you are looking for a numerical solution. I will also assume that you have an initial value problem. There are several situations in which this can arise: 1) The problem is complex valued, but not obviously so - in this case you need to use the 'complex'=true option. 2) The ODE system cannot be evaluated, as there is insufficient information. 3) Something else... Could you please post the system, and we can narrow it down further? - Allan
odeplot has the ability to plot expressions in terms of all returned derivatives (those returned when calling the dsolve numeric procedure), but your problem has higher derivatives. The below shows an approach for solving this:
> restart:
interface(warnlevel=0):
with(SolveTools):with( plots ): with( plottools ):
with( DEtools ): with( PDEtools ):
g:=9.81:
x[0][2]:=0:
y[0][2]:=0:
beams:=proc(n,la,ma)
global l,m,x,y:
l[n]:=la: m[n]:=ma:
x[n][2]:=0: y[n][2]:=0:
x[n][1]:=0: y[n][1]:=0:
end proc:
Force:=proc(n,fm,fn)
global fx,fy:
fx[n]:=fm:
fy[n]:=fn:
end proc:
beams(1,1,1):
beams(2,2,1):
beams(3,3,1):
beams(4,4,1):
Force(2,80,100):
> la:=proc(nn)
global x,y,q,vy,vx,T,lg,lm,lp,L,Ll,Lk,ax,ay:
n:=nn:
for i from 1 to n do:
x[i][2]:=l[i]*cos(q[i](t))+x[(i-1)][2]:
y[i][2]:=l[i]*sin(q[i](t))+y[(i-1)][2]:
x[i][1]:=x[i][2]-0.5*l[i]*cos(q[i](t)):
y[i][1]:=y[i][2]-0.5*l[i]*sin(q[i](t)):
end do:
#generate the velocity of the center of the beams
for i from 1 to n do:
vx[i]:=diff(x[i][1],t):
vy[i]:=diff(y[i][1],t):
ax[i]:=diff(x[i][1],t$2):
ay[i]:=diff(y[i][1],t$2):
end do:
T:=0:
for i from 1 to n do:
T:=T+0.5*m[i]*(vx[i]^2+vy[i]^2)+1/24*m[i]*l[i]^2*diff(q[i](t),t)^2:
end do:
for i from 1 to 4 do:
T:=subs(diff(q[i](t),t)=omega[i],T);
T:=subs(q[i](t)=theta[i],T):
end do;
for i from 1 to n do:
lg[i]:=diff(T,omega[i]):
for j from 1 to n do:
lg[i]:=subs({theta[j]=theta[j](t),omega[j]=omega[j](t)},lg[i]):
end do;
end do:
for p from 1 to n do
lp[p]:=diff(lg[p],t):
end do:
for o from 1 to n do:
lm[o]:=diff(T,theta[o]):
for k from 1 to n do:
lm[o]:=subs({theta[k]=theta[k](t),omega[k]=omega[k](t)},lm[o]):
end do:
end do:
for k from 1 to n do:
L[k]:=-lp[k]+lm[k]+(fy[k]+0.5*m[k]*g)*cos(theta[k](t))*l[k]-fx[k]*sin(theta[k](t))*l[k]=0:
#L[k]:=-lp[k]+lm[k]+(fy[k]+0.5*m[k]*g)*cos(theta[k](t))*l[k]-fx[k]*sin(theta[k](t))*l[k]
#-(fy[k]+0.5*m[k]*g)*theta[k](t)*sin(theta[k]##(t))*l[k]-fx[k]*cos(theta[k](t))*l[k]*theta[k](t)=0:
#Ll[k]:=lp[k]*(1/l[k])*cos(theta[k](t))*omega[k](t)/sin(theta[k](t))^2-lm[k]/(l[k]*sin(theta[k](t)))-fx[k]:
#Lk[k]:=lp[k]*(1/l[k])*sin(theta[k](t))*omega[k](t)/cos(theta[k](t))^2+lm[k]/(l[k]*cos(theta[k](t)))-fy[k]-m[#k]*g:
Ll[k]:=fx[k-1]-fx[k]=m[k]*ax[k]:
Lk[k]:=fy[k-1]-fy[k]-0.5*m[k]*g=m[k]*ay[k]:
end do:
for s from 1 to n do:
for r from 1 to n do:
L[s]:=subs(diff(omega[r](t),t)=diff(theta[r](t),t,t),L[s]):
#Lk[s]:=subs(diff(omega[r](t),t)=diff(theta[r](t),t,t),Lk[s]):
Lk[s]:=subs(q[r](t)=theta[r](t),Lk[s]):
Ll[s]:=subs(q[r](t)=theta[r](t),Ll[s]):
#Ll[s]:=subs(diff(omega[r](t),t)=diff(theta[r](t),t,t),Ll[s]):
L[s]:=subs(omega[r](t)=diff(theta[r](t),t),L[s]):
#Lk[s]:=subs(omega[r](t)=diff(theta[r](t),t),Lk[s]):
#Ll[s]:=subs(omega[r](t)=diff(theta[r](t),t),Ll[s]):
L[s]:=combine(L[s],trig):
Ll[s]:=combine(Ll[s],trig):
Lk[s]:=combine(Lk[s],trig):
end do:
end do:
end proc:
> la(2):
solve_la:=proc(nn,inc)
global L,Lk,Ll,dsys:dsys2,ini,dsn:
ini:=inc:
n:=nn:
dsys:={seq(Lk[i],i=1..n),seq(Ll[i],i=1..n)}:
ss:=solve(dsys,{seq(fx[i],i=0..n-1),seq(fy[i],i=0..n-1)}):
dsys2:=eval({seq(L[i],i=1..n)},ss);
dsn := dsolve({op(dsys2)} union ini,numeric,implicit=true);
end proc:
n:=2:
j:=n-1:
#Lk[1];
#Lk[2];
#Ll[1];
#Ll[2];
#L[1];
#L[2];
#dsys:={Lk[2],Lk[1],Ll[1],Ll[2],L[1],L[2]};
dsys:={seq(Lk[i],i=1..n),seq(Ll[i],i=1..n)};
#dsys:={Lk[2],Lk[1],Ll[1],Ll[2]};
> ss:=solve(dsys,{seq(fx[i],i=0..n-1),seq(fy[i],i=0..n-1)});
> dsys2:=eval({seq(L[i],i=1..n)},ss);
> ini := {theta[1](0) = 1/6*Pi, theta[2](0) = 1/3*Pi, D(theta[1])(0) = 0, D(theta[2])(0) = 0};
dsn := dsolve({op(dsys2)} union ini,type=numeric,output=listprocedure);
Solution starts here
odeplot has the ability to plot expressions of all derivatives explicitly returned from the dsolve numeric procedure, but only those.
In your problem, the forces fx[0],fx[1],fy[0],fy[1] also depend upon second derivatives of the theta, so we would like to solve for these, and eliminate them in terms of the derivatives returned from dsolve as follows:
> dsys3 := combine(solve(convert(dsys2,rational),{diff(theta[1](t),t,t),diff(theta[2](t),t,t)}),trig);
Note the convert,rational was just to prevent introduction of extra terms due to round-off.
Now that we have the second derivatives, substitute into the forces to obtain:
> forces := subs(dsys3,subs(ss,[fx[0],fx[1],fy[0],fy[1]]));
> indets(forces,'function');
So we can see that there are at worst first derivatives present, which are returned by dsolve, so we can proceed with plotting:
> plots[odeplot](dsn,[[t,theta[1](t)],[t,theta[2](t)]],0..1);
Now plot fx[0] in red, and fx[1] in green
> plots[odeplot](dsn,[[t,forces[1]],[t,forces[2]]],0..1);
Now plot fy[0] in red, and fy[1] in green:
> plots[odeplot](dsn,[[t,forces[3]],[t,forces[4]]],0..1);
>
This post generated using the online HTML conversion tool
Download the original worksheet
Here's the original equation
> eq1:=diff(v(t),t)+4.646e-4*v(t)^3-7.13e-3*v(t)^2-0.085*v(t)-2.82;
Now relate distance 's(t)' and velocity 'v(t)'
> eq2 := v(t)=diff(s(t),t);
Now since we want the distance travelled from t=3.35 to t=4.77, and have the velocity at t=3.35, we set s(3.35) to zero and then s(4.77) will give us the final displacement from zero:
> ini := {s(3.35)=0,v(3.35)=14.77};
Now solve and look at a plot
> dsn := dsolve({eq1,eq2} union ini,numeric);
plots[odeplot](dsn,[t,s(t)],3.35..4.77);
Now since the displacement is monotonic, it also corresponds to the distance travelled, so we have:
> dsn(4.77);
Or that the total distance travelled is 25.04
This post generated using the online HTML conversion tool
Download the original worksheet
Hello,
You need to make the behavior of the W(i,k,j,l) explicit in the input system if they depend upon either the independent or dependent variables of the problem.
If you can define the function W(i,k,j,l) in advance of forming the equations, then all should work as expected.
Alternatively, for numerical solution, you can define the problem procedurally, defining your system in terms of a procedure, in which case any of the B[i,j] that the W depend upon must be explicitly present in the argument list to W.
Can you provide an example of some of the W(i,k,j,l) you would like to use in practice?
Thanks,
Allan Wittkopf
What initial conditions do you have?
Do you want an exact solution, or numerical solution?
This looks like a duplicate of another thread.
The approach is outlined in http://beta.mapleprimes.com/forum/how-to-solve-implicit-ode-system-like-this-0
You need to do a little manual work to make the input recognizable as a differential system, then use the 'implicit' option of dsolve/numeric to obtain a numerical solution. By the way, equation 3 is a duplicate of equation 1. So what I did was to (after removing equation 3) solve the first 4 equations for fx[0],fy[0],fx[1],fy[1], and substitute into the last two, obtaining equations in the actual differential variables of the problem only (theta[1](t) and theta[2](t)). I then invoked dsolve(...,numeric,implicit=true) to use a simple implicit numerical approach (this approach only works when the ODE system is linear in it's highest order derivatives). You will see in the included worksheet snip (below) that the solution is easily obtained and examined.
Here is the Maple input for the system:
> dsys := [fy[1]-100 = -1.0*sin(theta[2](t))*diff(theta[2](t),t)^2+1.0*cos(theta[2](t))*diff(theta[2](t),`$`(t,2))-sin(theta[1](t))*diff(theta[1](t),t)^2+cos(theta[1](t))*diff(theta[1](t),`$`(t,2)),
fy[0]-fy[1] = -.5*sin(theta[1](t))*diff(theta[1](t),t)^2+.5*cos(theta[1](t))*diff(theta[1](t),`$`(t,2)),
fx[0]-fx[1] = -.5*cos(theta[1](t))*diff(theta[1](t),t)^2-.5*sin(theta[1](t))*diff(theta[1](t),`$`(t,2)),
fx[1]-80 = -1.0*cos(theta[2](t))*diff(theta[2](t),t)^2-1.0*sin(theta[2](t))*diff(theta[2](t),`$`(t,2))-cos(theta[1](t))*diff(theta[1](t),t)^2-sin(theta[1](t))*diff(theta[1](t),`$`(t,2)),
-2.083333333*diff(theta[1](t),`$`(t,2))-2.000000000*diff(theta[2](t),t)^2*sin(theta[1](t)-theta[2](t))-2.000000000*diff(theta[2](t),`$`(t,2))*cos(theta[1](t)-theta[2](t))+cos(theta[1](t))*fy[1]+4.905*cos(theta[1](t))-fx[1]*sin(theta[1](t)) = 0,
-160*sin(theta[2](t))+209.810*cos(theta[2](t))-4.333333333*diff(theta[2](t),`$`(t,2))+2.000000000*diff(theta[1](t),t)^2*sin(theta[1](t)-theta[2](t))-2.000000000*diff(theta[1](t),`$`(t,2))*cos(theta[1](t)-theta[2](t)) = 0];
Now we need an actual differential system, so we explcitly remove the auxiliary variables:
> ss := solve(dsys[1..4],{fx[0],fx[1],fy[0],fy[1]});
> dsys := eval(dsys[5..6],ss);
So now we have 2 equations in 2 differential unknowns, so we can solve the problem.
Here's the ICs:
> ini := {theta[1](0) = 1/6*Pi, theta[2](0) = 1/3*Pi, D(theta[1])(0) = 0, D(theta[2])(0) = 0};
And use this call to solve - note the 'implicit=true' option.
> dsn := dsolve({op(dsys)} union ini,numeric,implicit=true);
> plots[odeplot](dsn,[[t,theta[1](t)],[t,theta[2](t)]],0..1);
>
This post generated using the online HTML conversion tool
Download the original worksheet
Hello Alejandro,
In the setup in the worksheet, the timing measurements for dsolve are actually setting up and applying the numerical problem every time you request the solution.
Given that, it's actually quite fast ;-)
So basically at every requested point the worksheet calls the top-level dsolve/numeric which has to:
1) Process the input problem
2) Build the numerical evaluator
3) Compute the interpolant over 0..2
4) Evaluate the interpolant for the requested value
Ideally you want to set up the problem once for the values of beta,lambda,sigma,P,B (steps 1-3), then use this for multiple 'z' values (step 4).
This can be accomplished by replacing your definition for DLn2:
DLn2:=(z,beta,lambda,sigma,P,B)->
(1+z)*sqrt((1+sigma)^P+B)*rhs(op(2,rr(beta,lambda,sigma,P,B)))(z);
With the following:
DLn2:=proc(z,beta,lambda,sigma,P,B) global DLn2state;
if not assigned(DLn2state) or DLn2state[1]<>[beta,lambda,sigma,P,B] then
DLn2state := [[beta,lambda,sigma,P,B],
rhs(op(2,rr(beta,lambda,sigma,P,B)))];
end if;
(1+z)*sqrt((1+sigma)^P+B)*DLn2state[2](z);
end proc:
Now in this I am using the global variable DLn2state to store the dsolve/solution once computed, then using it until the parameters change, at which point the solution is rebuilt.
This is not the preferred method, as when you set this up in a main program you probably want to obtain the solution procedure once, then use it multiple times, not using a global as I have done above.
With this change the timings become much better, and the dsolve approach is a clear winner.
You may also want to consider the 'output=piecewise' option of dsolve/numeric for this problem, as then you can get a piecewise solution.
And finally, if the z values at which you want to compute the solution are known in advance, the output=Array option (see ?dsolve,numeric) is by far the fastest.
- Allan Wittkopf
First off, the problem as stated is more difficult than it ought to be. This is a reaction process, so this means that since we have a+b=c, we also have: a(0)-a(t)=b(0)-b(t) (reactants must combine equally) a(0)-a(t)+c(0)-c(t)=0 (balance between reactants and result) From this a simplified equation in 'c' alone can be obtained: > s1 := solve({a(0)-a(t)=b(0)-b(t),a(0)-a(t)+c(0)-c(t)=0},{a(t),b(t)}); s1 := {b(t) = c(0) - c(t) + b(0), a(t) = a(0) + c(0) - c(t)} > dsys := { > diff(a(t),t)=-k1*a(t)*b(t)+k2*c(t), > diff(b(t),t)=-k1*a(t)*b(t)+k2*c(t), > diff(c(t),t)= k1*a(t)*b(t)-k2*c(t) > }: > dsys := eval(dsys,s1); /d dsys := {-|-- c(t)| = -k1 (a(0) + c(0) - c(t)) (c(0) - c(t) + b(0)) + k2 c(t), dt / d -- c(t) = k1 (a(0) + c(0) - c(t)) (c(0) - c(t) + b(0)) - k2 c(t)} dt > dsys[1]+dsys[2]; 0 = 0 So we have the one equation: > de := dsys[1]; /d de := -|-- c(t)| = -k1 (a(0) + c(0) - c(t)) (c(0) - c(t) + b(0)) + k2 c(t) dt / Utilizing the initial data gives: > ics := {a(0)=10,b(0)=10,c(0)=0}; ics := {a(0) = 10, b(0) = 10, c(0) = 0} > de := eval(de,ics); /d 2 de := -|-- c(t)| = -k1 (10 - c(t)) + k2 c(t) dt / > de := eval(de,ics); /d 2 de := -|-- c(t)| = -k1 (10 - c(t)) + k2 c(t) dt / From which we can obtain an *exact* solution for the problem: > sol := dsolve({de,c(0)=0}); / | sol := c(t) = 1/2 |20 k1 + k2 - | 1/2 t (k2 (k2 + 40 k1)) 20 k1 + k2 tanh(---------------------- + arctanh(--------------------)) 2 1/2 (k2 (k2 + 40 k1)) 1/2| (k2 (k2 + 40 k1)) |/k1 | / Now from the exact solution, and the data we can do: > X := Vector([seq(i,i=0..25)]): > Y := Vector([0, 4.981119825867072, 6.619863704639184, 7.424230767960948, > 7.895291258366989, 8.200043893374819, 8.410083696006469, 8.561234239001111, > 8.673418582102649, 8.75860469142499, 8.824418203860477, 8.875947188583103, > 8.916714710357871, 8.949234971834077, 8.975347266499856, 8.996425002439704, > 9.01351139331151, 9.027409882146703, 9.038747104014094, 9.048016122443334, > 9.055608432249834, 9.061836841504377, 9.066952814916256, 9.071159284138764, > 9.074620916958418, 9.07747159947345]): > Statistics[NonlinearFit](rhs(sol),X,Y,[t],parameternames=[k1,k2]); Error, (in Statistics:-NonlinearFit) complex value encountered Oops, no go, try: > Statistics[NonlinearFit](Re(rhs(sol)),X,Y,[t],parameternames=[k1,k2]); bytes used=8001568, alloc=4979824, time=0.50 bytes used=12002500, alloc=5307444, time=0.71 bytes used=16003524, alloc=5504016, time=0.93 bytes used=20004628, alloc=5504016, time=1.15 10.0454545440218 - 0.954545439467380 Re(tanh(0.0954545522498140 t + 0.0953101783033397 - 1.57079632679490 I)) And we get the function, and note that the main problem with the complex value is the fact that there is a small complex round-off error causing all the trouble. The comparison of the collected data and the computed solution looks quite good.