How to plot a differential equation?

How to plot the differential equation

 (x-2/3)*f'(x)=6*f(x)-2+[5*(x-1/3)^+]-5*f(min([x+1/3,2/3]))?
  f(0)=0.8, x from 0 to 2/3

[Torsten]

Please insert multiplication signs and clarify what is meant by [5(x-1/3)^+].

Best wishes

Torsten.

[Chong Zhang]

(x-1/3)^+ means if x>1/3,(x-1/3)^+=x-1/3, otherwise,(x-1/3)^+=0. I don't know how to express it.

[Torsten]

And what does

f(min([x+1/3,2/3]))

mean ?

Is min(x+1/3,2/3) really the argument for f ?

[Chong Zhang]

f(min([x+1/3,2/3])) means if x>1/3,f(min([x+1/3,2/3]))=2/3, otherwise,f(min([x+1/3,2/3]))=x+1/3. Sorry for confusing you.

[Torsten]

So your equation reads

(x-2/3)*f'(x)=6*f(x)-2+5*(x-1/3)-5*f(2/3) if 1/3 <= x <= 2/3
(x-2/3)*f'(x)=6*f(x)-2-5*f(x+1/3)         if 0 <= x <= 1/3

and f(2/3), f(x+1/3) really means: f evaluated at 2/3 and f evaluated at x+1/3, respectively ?

[Chong Zhang]

Yes! And f(2/3), f(x+1/3) are what I don't know how to code.

[Torsten]

1. Assume a value for f(1/3) and name it "fmiddle".

2. Solve the differential equation on the interval 1/3 <= x <= 2/3 using bvp4c with f(2/3) as a free parameter.

3. Solve the differential equation on the interval 0 <= x <=1/3 using ODE45 by using the solution from 2 to evaluate f(x+1/3).

4. Compare f(1/3) obtained from the solution in 3. and "fmiddle". If abs(f(1/3)-fmiddle) < tol, accept the solution for f. Otherwise update "fmiddle" and go to 2.

Best wishes

Torsten.

[Chong Zhang]

Will try. Thanks!

Answer 1

[Torsten]

function main

% call root finder
estimate0 = 1.0;
estimate = fzero(@cycle,estimate0);

% call bvp4c for final value for f(1/3)
lambda = 1;
eps = 0.000000001;
solinit = bvpinit(linspace(1/3-eps,2/3-eps,20),@(x)mat4init(x,lambda,estimate),lambda);
sol1 = bvp4c(@mat4ode,@(ya,yb,lambda)mat4bc(ya,yb,lambda,estimate),solinit);

% call ode45 for final value for f(1/3)
x0 = 0.8;
tspan = [0 1/3-eps];
sol2 = ode45(@(x,y)fun_ode45(x,y,sol1),tspan,x0);

% plot entire curve
x1 = linspace(0,1/3-eps,20);
S1 = deval(sol2,x1);
x2 = linspace(1/3-eps,2/3-eps,20);
S2 = deval(sol1,x2);
plot(x1,S1,x2,S2)

end

% function to calculate f(1/3)  
function ret = cycle(estimate)

% call bvp4c
lambda = 1;
eps = 0.000000001;
solinit = bvpinit(linspace(1/3-eps,2/3-eps,20),@(x)mat4init(x,lambda,estimate),lambda);
sol = bvp4c(@mat4ode,@(ya,yb,lambda)mat4bc(ya,yb,lambda,estimate),solinit);  

%call ode45
x0 = 0.8;
tspan = [0 1/3-eps];
[X,Y]=ode45(@(x,y)fun_ode45(x,y,sol),tspan,x0);
ret = Y(end)-estimate;

end

% functions for bvp4c
% ------------------------------------------------------------
function dydx = mat4ode(x,y,lambda)
dydx = (6*y(1)-2+5*(x-1/3)-5*lambda)/(x-2/3); 
end
% ------------------------------------------------------------
function res = mat4bc(ya,yb,lambda,estimate)
res = [ya(1)-estimate ; yb(1)-lambda];
end 
% ------------------------------------------------------------
function yinit = mat4init(x,lambda,estimate)
yinit = estimate+3*(x-1/3)*(lambda-estimate);
end

% function for ode45
function dydx = fun_ode45(x,y,sol)
interpolation = deval(sol,x+1/3);
dydx =(6*y(1)-2-5*interpolation)/(x-2/3);
end

Dirty code, but it works.

Best wishes

Torsten.

[Chong Zhang]

Perfect！