# How to plot a differential equation?

32 views (last 30 days)
Chong Zhang on 15 May 2018
Commented: Chong Zhang on 17 May 2018
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

Chong Zhang on 16 May 2018
Yes! And f(2/3), f(x+1/3) are what I don't know how to code.
Torsten on 16 May 2018
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 on 16 May 2018
Will try. Thanks!

Torsten on 17 May 2018
Edited: Torsten on 17 May 2018
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.

#### 1 Comment

Chong Zhang on 17 May 2018
Perfect！