NDSolve with integral equation
Toshiyuki (Toshi) Meshii
NDSolve seems to have difficulties with solving integral equation.
n = 5; NDSolve[{D[\[Sigma]norm[z, t], t] == 3*z*Integrate[\[Sigma]norm[z,
t]^n*z, {z, 0, 1}] - \[Sigma]norm[z, t]^n,
\[Sigma]norm[z, 0] == 1.5*z, \[Sigma]norm[0, t] == 0}*\[Sigma]norm[z,
t], {z, 0.01, 1}, {t, 0.01, 2}]
Mathematica returns a message
NDSolve::"deql": "The first argument must have both an equation and an \
initial condition."
which I cannot understand.
Can anybody tell what's wrong with my attempt?
-Toshi
Alois Steindl
Hello,
I am not really sure that your question isn't a joke.
NDSolve solves differential equations, not integral equations.
Although there are relations between these topics, they are certainly
not the same.
Alois
--
Alois Steindl, Tel.: +43 (1) 58801 / 32558
Inst. for Mechanics II, Fax.: +43 (1) 58801 / 32598
Vienna University of Technology,
A-1040 Wiedner Hauptstr. 8-10 Email: Alois.Ste...@tuwien.ac.at
Andrzej Kozlowski
*(Times) where a , (comma) should be. Presumably your input ought to be:
n = 5; NDSolve[{D[\[Sigma]norm[z, t], t] == 3*z*Integrate[\[Sigma]norm[z,
t]^n*z, {z, 0, 1}] - \[Sigma]norm[z, t]^n,
\[Sigma]norm[z, 0] == 1.5*z, \[Sigma]norm[0, t] ==
0},\[Sigma]norm[z,
t], {z, 0.01, 1}, {t, 0.01, 2}]
However, even in the corrected version the equation can't be solved.
First of all you will get the complaint:
NDSolve::"bcedge": "Boundary conditions must be specified at the edge of
the \
spatial domain."
In other words Mathematica wants a boundary condition for
\[Sigma]norm[z, 0.1] or alternatively you should use {z,0,1} in NDSolve.
But actually I do not think this equation is solvable by any numerical
scheme even if you could provide the initial conditions at the edge of
the boundary that Mathematica requests. To evaluate the integral in your
equation NDSolve needs to know the values of \[Sigma]norm[z, t] for all
z between 0 and 1 and a given t, but this knowledge is not available at
any stage of the evaluation. I am not really an expert, but this seems
to me a clear example of an equation that is not solvable by any
numerical means.
By the way, the fact that you know a solution to a differential
equation, and even the fact that the solution is very simple does not
imply that the equation can be solved by any known method, except of
course guessing, which computer programs generally do not use.
Andrzej Kozlowski
Toyama International University
JAPAN
On Friday, August 23, 2002, at 05:25 AM, Toshiyuki ((Toshi)) Meshii
wrote:
Bill Rowe
Meshii) wrote:
>NDSolve seems to have difficulties with solving integral equation.
>n = 5; NDSolve[{D[\[Sigma]norm[z, t], t] ==
>3*z*Integrate[\[Sigma]norm[z, t]^n*z, {z, 0, 1}] - \[Sigma]norm[z,
>t]^n, \[Sigma]norm[z, 0] == 1.5*z, \[Sigma]norm[0, t] ==
>0}*\[Sigma]norm[z, t], {z, 0.01, 1}, {t, 0.01, 2}]
>Mathematica returns a message
>
>NDSolve::"deql": "The first argument must have both an equation and an
>\ initial condition."
A general solution to a differential equation includes a constant of integration. It isn't possible to arrive a numerical solution without providing sufficient information to determine this constant. Specifying the initial condition provides this information.
Look at the examples using the help browser to see what NDSolve needs.
Daniel Lichtblau
>
> Hello,
>
> NDSolve seems to have difficulties with solving integral equation.
>
> n = 5; NDSolve[{D[\[Sigma]norm[z, t], t] == 3*z*Integrate[\[Sigma]norm[z,
> t]^n*z, {z, 0, 1}] - \[Sigma]norm[z, t]^n,
> \[Sigma]norm[z, 0] == 1.5*z, \[Sigma]norm[0, t] == 0}*\[Sigma]norm[z,
> t], {z, 0.01, 1}, {t, 0.01, 2}]
>
> Mathematica returns a message
>
> NDSolve::"deql": "The first argument must have both an equation and an \
> initial condition."
>
> Can anybody tell what's wrong with my attempt?
>
> -Toshi
I'll show a couple of approaches to this sort of problem. The first is
to iteratively solve for sn[k][z,t] in terms of sn[k-1][z,t], with
sn[0][z,t] initialized to something appropriate (I used 3/2*z). The sn's
are computed as interpolating functions based on results from solving
ODEs in t for many fixed values of z. The code below does this for n=2.
The plots at first appear to flip between two states but eventually
stabilize. Note that this takes many minutes to run to completion.
n = 2;
sn[0][z_,t_] := 3/2*z;
Do [
Do [
sn[k][z,t_] = sn[k][z,t] /.
First[NDSolve[{D[sn[k][z,t], t] ==
3*z*NIntegrate[sn[k-1][w,t]^n*w, {w,0,1}] - sn[k-1][z,t]^n,
sn[k][z,0] == 3/2*z}, sn[k][z,t], {t,0,2}]],
{z,0,1,1/100}
];
snew = Interpolation[Flatten[Table[{z,t,sn[k][z,t]},
{z,0,1,1/100}, {t,0,2,1/100}], 1], InterpolationOrder->7];
sn[k] = snew;
Print["iteration ", k];
Plot3D[sn[k][z,t], {z,0,1}, {t,0,2}],
{k, 1, 15}
]
One can test for convergence as below; it is apparently quite good.
NIntegrate[Abs[sn[15][z,t]-sn[14][z,t]], {z,0,1}, {t,0,2}]
A drawback to this method is that it appears to break down beyond n = 2.
Possibly one simply needs a much better starting function for sn[0], I'm
not sure.
Below is another method that Michael Trott showed me. We expand in a set
of basis functions in z, set up a system of ODEs in t, and solve them. I
tried this using monomials in z for basis functions and ran into some
trouble with the ODE solving, so I will show Michael's attempt using
trig functions for the basis. Note that we now handle the desired case
n=5; another advantage is that this is alot faster than the method
above, though still not exactly blinding in speed. Michael used the
interval {0,2} for z so as to achieve pointwise (not just L^2)
convergence; otherwise there would be a sharp drop-off just before z=1
as all the trigs vanish there. In other words, the basis functions are
of the form Sin[k/2*Pi] rather than Sin[k*Pi] for 1<=k<=deg. Cosines are
excluded due to the vanishing condition at z=0.
integrate[p_Plus, {z_, 0, 1}] := Integrate[#, {z, 0, 1}]& /@ p;
integrate[p_, {z_, 0, 1}] := Integrate[p, {z, 0, 1}];
n = 5;
deg = 6;
vars[t_] = Map[#[t]&,Array[a,deg]];
zFuns = Sin[Range[deg] Pi/2 z];
sn[z_,t_] = vars[t].zFuns;
eqs1 = 3*z*integrate[Expand[sn[w,t]^n*w],{w,0,1}] -
sn[z,t]^n - D[sn[z,t],t];
eqs2 = integrate[Expand[eqs1 #], {z, 0, 1}]& /@ zFuns;
iCs = integrate[Expand[(sn[z,0] - 3/2*z) #], {z, 0, 1}]& /@ zFuns;
fulleqns = # == 0& /@ Join[eqs2, iCs];
nsd = NDSolve[fulleqns, vars[t], {t,0,2}, SolveDelayed->True];
Plot3D[Evaluate[sn[z,t] /. nsd[[1]]],{z,0,1},{t,0,2}];
an initializer for the iteration/interpolation method above in an
attempt to refine the solution.
Daniel Lichtblau
Wolfram Research
Daniel Lichtblau
Steve,
Below is the note that had a line clipped (the one beginning ">From
here..."). I am sending it pretty much in the same way as last time,
with cc to math...@smc.vnet.net (only this time with the intent that
it not appear there).
Daniel
-------------------------------------------
>From here one might try increasing 'deg' or perhaps using the result as