Solving Integro-differential equation with limited integral

Hi,

 How can I solve this equation numerically using matlab

w''''=w''*int(w'^2,0,1)

I tried using the standard form of ODE function, the only problem I faced is how to represent that limited integral Thanks

Answer 1

[Torsten]

Write your integro-differential equation as

w1'=w2
w2'=w3
w3'=w4
w4'=w3*integral_{t=0}^{t=1}w2^2(t') dt'

Then discretize the interval [0:1] in n subintervals 0=t(1)<t(2)<...<t(n)=1.

Compute the derivatives as

wj'(t(i))=(wj(t(i+1))-wj(t(i)))/dt (j=1,2,3,4)

and compute the integral using the trapezoidal rule.

You'll arrive at a polynomial system (order 3) of equations for the unknowns

wj(t(2)),wj(t(3)),...,wj(t(n)) (j=1,2,3,4)

which can be solved by fsolve, e.g.

No chance to use ODE45 in this case.

Another way might be to use ODE45 and iteratively adjust the value of the integral, but I'm not sure whether this method will converge.

Good luck !

Best wishes

Torsten.

Answer 2

[Claudio Gelmi]

Take a look at this solver:

Article "IDSOLVER: A general purpose solver for nth-order integro-differential equations": http://dx.doi.org/10.1016/j.cpc.2013.09.008

MATLAB solver here: http://cpc.cs.qub.ac.uk/summaries/AEQU_v1_0.html

Best wishes,

Claudio

[Fernando Fernandes]

Gelmi help me! How can I use your method to solve this equation?

[Fernando Fernandes]

I've downloaded your paper, but i'm a beginner in Matlab. Do I need the solver in http://cpc.cs.qub.ac.uk/summaries/AEQU_v1_0.html ???

How can I install this?

Answer 3

[ash]

Thanks, sorry for the late reply

I tried to apply the technique you suggested (as much as I understood) using Euler 1st order, also my problem is a BVP (3 initial conditions and 1 BC) kindly find my code bellow The problem is that the code is too slow, and I can only solve small number of points

Is that what you advised me to do in our previous comment?, are there any enhancement for that code?

Thanks

syms a b

w=10;
T=30;
L=500;
F=1;
E=160e3;
I=2500;
A=T*w;
Npnts=11;
x=linspace(0,L/2,Npnts);

q1=sym(zeros(1,length(x)));
q2=sym(zeros(1,length(x)));
q3=sym(zeros(1,length(x)));
q4=sym(zeros(1,length(x)));

h=x(2)-x(1);
q1(1)=0;
q2(1)=0;
q3(1)=b;
q4(1)=-F/2/E/I;

for i=1:length(x)

      q1(i+1)=q1(i)+h*q2(i);
      q2(i+1)=q2(i)+h*q3(i);   
      q3(i+1)=q3(i)+h*q4(i);
      q4(i+1)=q4(i)+h*q3(i)*a*A/2/L/I; 

end
integ_a=sum(q2.^2)*h-a/2;
sol_ab=solve(integ_a==0,q2(i+1)==0,a,b);
sol_a=sol_ab.a;
sol_b=sol_ab.b;
sol_index=1;

q1=subs(q1,a,sol_a(sol_index));
q1=double(subs(q1,b,sol_b(sol_index)));

[ash]

Any help?

[Torsten]

The system to solve is

(w1(t(i+1))-w1(t(i)))/h = w2(t(i))
(w2(t(i+1))-w2(t(i)))/h = w3(t(i))
(w3(t(i+1))-w3(t(i)))/h = w4(t(i))
(w4(t(i+1))-w4(t(i)))/h = [sum_{j=1}^{j=Npnts-1}(w2(t(i+1))+w2(t(i)))*h/2]*w3(t(i))
(i=1,Npnts-1)

These are 4*(Npts-1) equations in which you will have to include the boundary conditions.

You can use fsolve to solve this system of polynomial equations.

Best wishes

Torsten.

[Torsten]

Sorry, should read

(w4(t(i+1))-w4(t(i)))/h = [sum_{j=1}^{j=Npnts-1}(w2(t(j+1))+w2(t(j)))*h/2]*w3(t(i))

Best wishes

Torsten.