Saltar al contenido

Resolver numéricamente una ecuación integro-diferencial

Siéntete libre de divulgar nuestra web y códigos en tus redes, ayúdanos a ampliar nuestra comunidad.

Solución:

Esta ecuación integro-diferencial se puede resolver con el método mencionado en esta respuesta, es decir, diferenciar la ecuación para convertirla en una puro ODA.

Primero, interprete las ecuaciones para Mathematica código. (Por cierto, si le hubieras dado el Mathematica forma de código de la ecuación en su pregunta, su pregunta habría atraído más atención. )

v = 1; ψ[ζ_] = ζ; f[ζ_, η_] = ζ + η; g[ζ_] = ζ^2;

bc[0] = x[1] == 1;

eq = -v D[x[ζ], ζ] + ψ[ζ] x[ζ] + 
    Integrate[f[ζ, η] x[η], η, 0, ζ] + g[ζ] x[1] == 0 /. Rule @@ bc[0]
(* ζ^2 + Integrate[(ζ + η)*x[η], η, 0, ζ] + ζ*x[ζ] - Derivative[1][x][ζ] == 0 *)

Luego, diferencia la ecuación dos veces.

neweq = D[eq, ζ, ζ]
(* 2 + 3*x[ζ] + 2*Derivative[1][x][ζ] + 2*ζ*Derivative[1][x][ζ] + 
     ζ*Derivative[2][x][ζ] - Derivative[3][x][ζ] == 0 *)

neweq es una EDO de tercer orden, mientras que actualmente solo tenemos 1 bcie bc[0], por lo que debemos deducir otros dos bcs de eq. Esto se puede lograr fácilmente configurando ζ para 0 en eq y D[eq, ζ]:

bc[1] = eq /. ζ -> 0
(* -Derivative[1][x][0] == 0 *)
[email protected] = D[eq, ζ] /. ζ -> 0
(* x[0] - Derivative[2][x][0] == 0 *)

Finalmente, resuelve la ecuación y encuentra $ x (0) $:

sol = NDSolveValue[neweq, bc /@ Range[0, 2], x, ζ, 0, 1]
sol[0]
(* 0.232727 *)

Actualización 1: solución alternativa para $ f ( zeta, eta) = e ^ ( zeta +1) eta $

Como señaló Carl Woll, cuando la forma de $ f $ se vuelve más complicada, puede ser imposible diferenciar la ecuación de una EDO pura. Aún así, existe una solución alternativa al menos para $ f ( zeta, eta) = e ^ ( zeta +1) eta $, es decir, aproximar $ f $ con su expansión de serie.

f[ζ_, η_] = Exp[(ζ + 1) η];

nmax = 10;
approx[ζ_, η_] = [email protected][f[ζ, η], ζ, 0, nmax];
(* Error Check: *)
Plot3D[f[ζ, x] - approx[ζ, x] // Evaluate, x, 0, 1, ζ, 0, 1, PlotRange -> All]

Gráficos de Mathematica

eq = -v D[x[ζ], ζ] + ψ[ζ] x[ζ] + 
     Integrate[approx[ζ, η] x[η], η, 0, ζ] + g[ζ] x[1] == 0 /. Rule @@ bc[0];

neweq = D[eq, ζ, nmax + 1];

bclst = Table[D[eq, ζ, n] /. ζ -> 0, n, 0, nmax];
sol = NDSolveValue[neweq, [email protected]
, bclst, x, ζ, 0, 1, WorkingPrecision -> 16]; sol[0] (* 0.1498546695665442 *)

Actualización 2: una solución alternativa más simple y probablemente más general

Resulta innecesario aproximar $ f $ con la serie de Taylor, podemos diferenciar la ecuación original suficientes veces y luego simplemente quita el resto Integrate[……]:

order = 6;

eq = -v D[x[ζ], ζ] + ψ[ζ] x[ζ] + 
     Integrate[f[ζ, η] x[η], η, 0, ζ] + g[ζ] x[1] == 
    0 /. Rule @@ bc[0];

rule = [email protected][__] :> 0;

neweq = D[eq, ζ, order + 1] /. rule;

bclst = Table[D[eq, ζ, n] /. ζ -> 0 /. rule, n, 0, order];

solnew = NDSolveValue[neweq, [email protected], bclst, x, ζ, 0, 1, 
    WorkingPrecision -> 16]; // AbsoluteTiming
(* 1.017286, Null *)
solnew[0]
(* 0.1498546688616941 *)

Pero, ¿por qué funciona esto? ¿No es solo una coincidencia?

No, no es.

Una explicación cualitativa es que, cada vez que se diferencia la ecuación, la Integrate[……] término jugará un papel menos importante en la nueva ecuación.

Una explicación cuantitativa es, si aproximamos $ f $ con un n-th orden polinomio por partes (por ejemplo, el polinomio aquí), luego después de diferenciar la ecuación para n+1 veces, el Integrate[…] el término será exactamente 0.

Aunque no lo he probado, creo que este enfoque es más general que el enfoque de expansión de Taylor, porque cuando la forma de $ f $ se vuelve aún más complicada (por ejemplo, por partes), la expansión de Taylor puede no ser adecuada.

Finalmente, una verificación de errores:

help[ζ_?NumericQ, sol_] := 
 NIntegrate[E^((1 + ζ) η) sol[η], η, 0, ζ, 
  Method -> Automatic, "SymbolicProcessing" -> 0]

test[sol_] := 
 Subtract @@ eq /. [email protected][__] -> help[ζ, sol] /. x -> sol
(* Please find cur[8] in Carl's answer *)
Plot[#, ζ, 0, 1] & /@ test /@ solnew, cur[8] // GraphicsRow

Gráficos de Mathematica

Como se puede ver, este enfoque es más preciso que el método iterativo dado por Carl.

Si puede convertir la ecuación integro-diferencial en un IVP o BVP ODE, entonces ese sería el mejor enfoque. Si no puede hacerlo, puede probar el siguiente enfoque iterativo. La idea básica es reemplazar el $ x ( eta) $ desconocido en el integrando con un ansatz (adivinar) y luego resolver la ecuación diferencial resultante. Repita, con la solución de la EDO como la próxima suposición que se utilizará.

Estos son los parámetros sugeridos:

ν = 1;
ψ[ζ_] := ζ
f[ζ_, η_] := ζ + η
g[ζ_] := ζ^2

Dada una conjetura cur[n-1], la siguiente suposición sería:

reset[] := (
    Clear[cur];
    cur[n_] := cur[n] = NDSolveValue[
        
        -ν x'[ζ] + ψ[ζ] x[ζ] + Integrate[f[ζ,η] cur[n-1][η], η, 0, ζ] + g[ζ] == 0,
        x[1] == 1
        ,
        x,
        ζ, 0, 1
    ]
)
reset[]

Comenzando con una suposición inicial:

cur[0] = #&;

Las siguientes gráficas muestran las diferencias en las iteraciones y la respuesta final:

Plot[[email protected][cur[n][t]-cur[n-1][t], n, 1, 8], t, 0, 1, PlotLegends->Range[8]]
Plot[cur[8][t], t,0,1]

ingrese la descripción de la imagen aquíingrese la descripción de la imagen aquí

Con la elección de $ f ( zeta, eta) = zeta + eta $, la ecuación integro-diferencial se puede convertir en una EDO, como lo muestra @xzczd, así que comparemos:

sol = NDSolveValue[
    
    2 + 3 x[ζ] + 2 x'[ζ] + 2 ζ x'[ζ] + ζ x''[ζ] - x'''[ζ] == 0,
    x[1] == 1,
    x'[0] == 0,
    x[0] - x''[0] == 0
    ,
    x,
    ζ, 0, 1
];

Aquí está la solución de la ODE:

Plot[sol[t], t, 0, 1]

ingrese la descripción de la imagen aquí

Y una trama de la diferencia:

Plot[sol[t]-cur[8][t],t,0,1]

ingrese la descripción de la imagen aquí

Como ejemplo de una ecuación integro-diferencial que no se puede convertir en una EDO (o al menos, no creo que pueda), considere:

f[ζ_, η_] := Exp[(ζ+1) η]
reset[] (* to clear out old solution *)
cur[0] = #&;
Plot[[email protected][cur[n][t]-cur[n-1][t], n, 1, 7], t, 0, 1, PlotLegends->Range[7]]
Plot[cur[7][t], t,0,1]

ingrese la descripción de la imagen aquíingrese la descripción de la imagen aquí

La solución iterativa de esta ecuación integro-diferencial converge bastante más lentamente.

Si te ha resultado de utilidad nuestro artículo, agradeceríamos que lo compartas con el resto entusiastas de la programación y nos ayudes a difundir nuestra información.

¡Haz clic para puntuar esta entrada!
(Votos: 0 Promedio: 0)



Utiliza Nuestro Buscador

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *