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
8 个评论
Torsten
2018-5-16
Please insert multiplication signs and clarify what is meant by [5(x-1/3)^+].
Best wishes
Torsten.
Chong Zhang
2018-5-16
Torsten
2018-5-16
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
2018-5-16
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
2018-5-16
编辑：Chong Zhang
2018-5-16
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.
Chong Zhang
2018-5-16
采纳的回答
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 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Boundary Value Problems 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
