You can use the 'known' option, as in the example below:

dsolve({f1(0)=0,diff(f1(x),x)=f1(x)+1},numeric,output=listprocedure);
F1:=subs(%,f1(x));
g:=x->2*F1(x/2);
g(1);
dsolve({f2(0)=1,diff(f2(x),x)=sin(f2(x))+g(x)},numeric,output=listprocedure,known=F1);
F2:=subs(%,f2(x));
plot(F2,0..3);
 

Preben Alsholm

You didn't assign values to the constants, but just wrote equations.

That is one reason why no plot turned up. The other is that you are trying to plot u(x,t) as a function of 2 variables using 'plot'.

I wouldn't assign values to the constants myself, but would use eval as seen below.

plot(eval(u(x, 2), {F = 1, c = 1, l = 1, m = 1, omega = 1}), x = 0 .. 20);
plots:-animate(plot,[eval(u(x, t), {F = 1, c = 1, l = 1, m = 1, omega = 1}), x = 0 .. 20],t=0..5);

One small detail: The line
WaveEqn1 := subs(WaveEqn = u(x, t)):
doesn't do anything. Your assignment of u has already done what you intend to do.

Preben Alsholm

It seems to me that your code works. The warning is just that. A warning, not an error. It is to be expected since the function xx(t-1) at the time is undefined for t > 2.

Using output=piecewise makes it unnecessary to rename, since there is no recursive calling. In fact you made the right hand side of the differential equation into a concrete function of t given by a formula.

So you could change your code to

restart;
xx := t -> 1  :
P := plot(xx, -1 .. 0):
xnum := dsolve({x(0) = xx(0), (D(x))(t) = -xx(t-1)},type = numeric, output = piecewise, range = 0 .. 1):
xx:=unapply(subs(op(xnum),x(t)),t):
P := P, plot(xx, 0 .. 1):
#Just to see that it works
plots:-display(P);
xnum := dsolve({x(1) = xx(1), (D(x))(t) = -xx(t-1)},type = numeric, output = piecewise, range = 1 .. 2):
xx := unapply(subs(op(xnum),x(t)),t):
#Just to see that it still works
P := P, plot(xx, 1 .. 2):
plots:-display(P);
 

Preben Alsholm

The problem may be that the name xx is reused in your second version.

The first of the following corresponds to your first loop and causes no problem.

The second redefines xx to be a procedure that must call a procedure also called xx. This is likely to be the problem.

I wouldn't think that avoiding array or Array would speed up your loop, anyway.

The unproblematic one:

restart;
theta := t -> 1;
xx := theta;
xnum := dsolve({x(0) = xx(0), diff(x(t),t) = -xx(t-1)},
type = numeric, output = listprocedure, 'known' = xx);

xx := rhs(xnum[2]);
xx(.12345);
xnum := dsolve({x(0) = xx(0), diff(x(t),t) = -xx(t-1)},
type = numeric, output = listprocedure, 'known' = xx);

xx1 := rhs(xnum[2]);
xx1(.12345);

The following will cause the loss of connection to the kernel.

restart;
theta := t -> 1;
xx := theta;
xnum := dsolve({x(0) = xx(0), diff(x(t),t) = -xx(t-1)},
type = numeric, output = listprocedure, 'known' = xx);

xx := rhs(xnum[2]);
xx(.12345);
xnum := dsolve({x(0) = xx(0), diff(x(t),t) = -xx(t-1)},
type = numeric, output = listprocedure, 'known' = xx);

xx := rhs(xnum[2]);
xx(.12345);
 

Preben Alsholm

You can do either one of the following. (I have inserted concrete values for a, b and t0).

The first method uses odeplot.

F := diff(y(x),x,x) + x^2*diff(y(x),x) + sin(x)*y(x) ;
EDF := {F = 0, y(1) = -.1, D(y)(1) = 2};
Snum := dsolve(EDF, y(x), type = numeric);
plots:-odeplot(Snum,[x,y(x)],0..10);
plots:-odeplot(Snum,[x,diff(y(x),x)],0..10);

The second method uses output=listprocedure.

Snum2 := dsolve(EDF, y(x), type = numeric,output=listprocedure);
S:=subs(Snum2,y(x));
S1:=subs(Snum2,diff(y(x),x));
plot(S,0..10);
plot(S1,0..10);
 

Preben Alsholm

If you know something about the ranges in which the answers should lie then you could use that information in fsolve.

(What are "the real" alfa values you mention?)

However, fsolve does come up with a solution without being given ranges or starting values:

EQ1:= 0.9868421053+74479.54250*a^4-1.*alfa-33391.41112*P*a^3-.5000000000*C2*a^2;

EQ2:= 2.9791817*10^5*a^3-0.4044127197e-1-1.001742334*10^5*P*a^2-1.*C2*a+0.4098048893e-1*A;

EQ3:= -2.003484667*10^5*P*a-0.3358800946e-2*A-1.*C2+8.937545100*10^5*a^2;

EQ4:= -2.003484667*10^5*P+0.1376453050e-3*A+0.1113310242e-4*a+0.1358341826e-3;

fsolve(eval({EQ1,EQ2,EQ3,EQ4},a=7),{C2,alfa,A,P});

Preben Alsholm
 

Make the plots (or animations) p1 and p2.

Then

display(Array([p1,p2]));

Preben Alsholm

One way of doing it is to make Y a function of two variables, theta and X:

Y:= (theta,X) -> 1 - ( sin (2*theta)/ (X+cos(2*theta)))^2;
Y(Pi/9,.12345);
evalf(%);
plot(Y(Pi/9,X),X=0..30);
plot(Y(Pi/9,1/X),X=0..30);
plot(Y(Pi/9,1/X^2),X=0..30);

Notice the spelling of Pi (not pi).

Preben Alsholm

Robert Israel beat me to it, but I'll let the following stay.

Assuming that the values x[i] are sorted increasingly, you can do

piecewise(t < x[1], hz[1], t < x[2], hz[2], t< x[3], hz[3], t < x[4], hz[4], t < x[5], hz[5],t < x[6], hz[6], t < x[7], hz[7], t< x[8], hz[8], t < x[9], hz[9], t <= x[10], hz[10]);
 

since piecewise checks the inequalities from the left and stops if the inequality is satisfied.

Even shorter is

piecewise(op(map(op,[seq([t<x[i],h[i]],i=1..9)])),t<=x[10],hz[10]);

where I have singled out the last one since you want <= (if you really do).

Notice that I have used t instead of x. It is not a good idea to use x as a variable name if x[i] is in use.

Preben Alsholm

I don't think you are missing anything.

I tried

F1:=combine(Fder):

solve(F1=0, h);
                                      h
i.e. output is h (!). However, the expression F1 is huge. But certainly this answer is strange.
If you plot F1 you run into numerical problems, which is not surprising.

plot(F1,h=0..Pi,-1.5..1.5);
 

Preben Alsholm

Try executing the two sums

sum(f,i=-infinity..infinity);
sum(f(i),i=-infinity..infinity);

and ask yourself if those sums really are what you want to plot. 

Preben Alsholm

Just one thing.

In your example the first term in the sum is p*ln(0), which prevents solve from doing anything in solve( eqn[21,2], s);

Preben Alsholm

Since

eval(rhs(eqn13),{f(z)=1,e(z)=1,z=0});
evaluates to

-6210.378246*I

i.e. an imaginary number, you may want to rethink the mathematics of the problem. Is there a wrong sign under the squareroot in eqn13?

Preben Alsholm

You are trying to assign to a sequence of sums.

In the help page for assign it says that you can assign to names or functions, where 'function' means an unevalated  function call, like f(8).

Preben Alsholm

I have not changed much of your original code.

restart;
with(plots):
S1:=unapply(u*FresnelS(u)+((1/Pi)*cos((1/2)*Pi*u^2))-(1/Pi),u);
C1:=unapply(u*FresnelC(u)-((1/Pi)*sin((1/2)*Pi*u^2)),u);
R1:=1;
u1:=evalf(sqrt(0.4));
E1:=evalf(Pi*R1*S1(u1)+R1);
D1:=evalf(Pi*R1*C1(u1)+1);
a:=evalf(Pi*R1*(u1));
p1:=plot([1+a*FresnelC(u),a*FresnelS(u)-E1,u=0..u1]):
p2:=plot(-1.064877386,x=0..1):
p3:=plot(1.064877386,x=0..1):
p4:=plot([1+a*FresnelC(u),-a*FresnelS(u)+E1,u=0..u1]):
p5:=plot(-1.064877386,x=-1..0):
p6:=plot(1.064877386,x=-1..0):
p7:=plot([-1+(-a*FresnelC(u)),a*FresnelS(u)-E1,u=0..u1]):
p8:=plot([-1+(-a*FresnelC(u)),-a*FresnelS(u)+E1,u=0..u1]):
X:=1+a*FresnelC(u1);
t1:=fsolve(D1+R1*cos(t)=X,t);
C1:=plot([D1+R1*cos(t),R1*sin(t),t=-t1..t1]):
C2:=plot([-D1-R1*cos(t),R1*sin(t),t=-t1..t1]):
display(p7,p1,p2,p3,p4,p5,p6,p8,C1,C2,scaling=constrained);
 

Preben Alsholm