ODE45, differential equation

조회 수: 4 (최근 30일)
Takey Asaad
Takey Asaad 2018년 9월 14일
편집: Stephan 2018년 9월 17일
my function is
dy/dt=k*y*exp(450/y)
k is constant and y(0)=40 and y(15)=95 solve this equation by using ode45 can someone pleaseeeeeeeeeee check the code and make it work .
tspan = [0 300];
y0 = 40;y15=95
[t,y] = ode45(@(t,y) 'k'*y*exp(450/y), tspan, y0,y15);
plot(t,y,'-o')

채택된 답변

Stephan
Stephan 2018년 9월 16일
편집: Stephan 2018년 9월 17일
Hi,
if you search a value for k that complies with the boundary conditions, you could use fzero to solve the problem numerically:
k = fzero(@calculate_k, 0.001);
disp(['k = ', num2str(k)])
[t,y] = ode45(@(t,y) k*y*exp(450/y), [0 300], 40);
plot(t,y,'r')
hold on
scatter(15,95,'ob')
text(23,96,'y(t=15) = 95','Color','b')
hold off
function k_value = calculate_k(x)
tspan = [0 15];
y0 = 40;
[~,y] = ode45(@(t,y) x*y*exp(450/y), tspan, y0);
k_value = y(end) - 95;
end
This will give you as result:
k = 0.00010435
or if you need it more accurate:
k =
1.043495007807761e-04
and the plot which belongs to this result looks like the boundary conditions are met (To be precise, it keeps the boundary conditions. This is because fzero is a powerful tool and also because you can think of a very good initial value for x0 with just a few estimates...):
.
Best regards
Stephan
  댓글 수: 2
Takey Asaad
Takey Asaad 2018년 9월 17일
Thanks so much
Stephan
Stephan 2018년 9월 17일
편집: Stephan 2018년 9월 17일
No problem - please be patient with your questions. Dont ask the same question a second time, if it needs more time than you expected to get an answer.
Better try to give as much as and as good as possible informations that are needed to let the contributers give you a useful answer.

댓글을 달려면 로그인하십시오.

추가 답변 (3개)

James Tursa
James Tursa 2018년 9월 15일
Remove the quotes from 'k', and be sure to define k before you call ode45. Also, ode45 is an initial value problem solver, so the y15 variable is not applicable (remove it from the call).
  댓글 수: 9
Takey Asaad
Takey Asaad 2018년 9월 15일
편집: James Tursa 2018년 9월 17일
function bvp4
xlow=0;xhigh=300;
solinit=bvpinit(linspace(xlow,xhigh,300),[40 95]);
sol=bvp4c(@bvp4ode,@bvp4bc,solinit);
xinit=linspace(xlow,xhigh,20);
sxint=deval(sol,xint);
plot(xint,sxint(1,:);
function dydx=bvp4ode(x,y)
dydx=[y(2) a*y*exp(450/y);
function res=bvp4bc(ya,yb)
res=[ya(1)-1,yb(1)
function bvp4
Error: Function definition not supported in this context. Create functions in code file.
Takey Asaad
Takey Asaad 2018년 9월 15일
편집: James Tursa 2018년 9월 17일
function bvp4
xlow=0;xhigh=300;
solinit=bvpinit(linspace(xlow,xhigh,300),[40 95]);
sol=bvp4c(@bvp4ode,@bvp4bc,solinit);
xinit=linspace(xlow,xhigh,20);
sxinit=deval(sol,xint);
plot(xinit,sxinit(1,:);
function dydx=bvp4ode(x,y)
dydx=[y(2) ,a*y*exp(450/y)];
function res=bvp4bc(ya,yb)
res=[ya(1)-1,yb(1)]

댓글을 달려면 로그인하십시오.


Takey Asaad
Takey Asaad 2018년 9월 15일
if true
function bvp4
xlow=0;xhigh=300;
solinit=bvpinit(linspace(xlow,xhigh,300),[40 95]);
sol=bvp4c(@bvp4ode,@bvp4bc,solinit);
xinit=linspace(xlow,xhigh,20);
sxinit=deval(sol,xint);
plot(xinit,sxinit(1,:);
function dydx=bvp4ode(x,y)
dydx=[y(2) ,a*y*exp(450/y)];
function res=bvp4bc(ya,yb)
res=[ya(1)-1,yb(1)]
end

Takey Asaad
Takey Asaad 2018년 9월 15일
code

카테고리

Help CenterFile 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