MATLAB Answers

How to plot a differential equation?

32 views (last 30 days)
Chong Zhang
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

  8 Comments

Show 5 older comments
Chong Zhang
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
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.

Sign in to comment.

Accepted Answer

Torsten
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.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!