How to plot a differential equation?

1 次查看（过去 30 天）
Chong Zhang
Chong Zhang 2018-5-15
评论： Chong Zhang 2018-5-17
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 个评论
显示 6更早的评论
Torsten
Torsten 2018-5-16
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2018-5-17
编辑：Torsten 2018-5-17
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.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by