Re: NDSolve with integral equation
- To: mathgroup at smc.vnet.net
- Subject: [mg36250] Re: [mg36144] NDSolve with integral equation
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Wed, 28 Aug 2002 04:16:10 -0400 (EDT)
- References: <200208230425.AAA00142@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
"Toshiyuki (Toshi) Meshii" wrote:
>
> 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."
>
> which I cannot understand.
> 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
- References:
- NDSolve with integral equation
- From: "Toshiyuki \(Toshi\) Meshii" <meshii@mech.fukui-u.ac.jp>
- NDSolve with integral equation