False Position method for 2nd order Differential equation BVP

1 次查看(过去 30 天)
I have a good code for the Runge-Kutta method to determine the function for a 2nd order differential equation, but I need to create a code that solves my boundary conditions. My boundaries are y(0)=25 and y(30)=20. I need to determine what y'(0)=? to get y(30)=20 using false position method. When y'(0)=0, y(30)=54.32. When y'(0)=-3, y(30)=-15.3. Basically what I need to do is use the y(30) values [-15.3,54.32] to find what y'(0) would give me y(30)=20 using false position. If anyone has any idea how to code this and can help I would appreciate it! I have put my RK code below. If you need to know more about the problem it is Problem 6 part d in the attached file.
function [KR4] = NewHighway(z1,z2)
C=.053;
f1 = z2 ;
f2 = C*sqrt(1+z2.^2);
KR4 = [f1 ; f2] ;
end
clear all ;
h = 0.01 ; %step size
t = 0:h:30;
a = zeros(1,length(t));
b = zeros(1,length(t));
a(1) = 25 ; %initial values
b(1) = -3 ;
for i = 1:(length(t)-1)
k1 = NewHighway(a(i), b(i)); %steps for RK method
k2 = NewHighway(a(i)+h/2*(k1(1)), b(i) + h/2*(k1(2)));
k3 = NewHighway((a(i)+h/2*k2(1)), b(i) + h/2*(k2(2)));
k4 = NewHighway((a(i)+h*k3(1)), b(i) + h*k3(2));
a(i+1) = a(i) + (h/6)*(k1(1)+ 2*k2(1) + 2*k3(1) + k4(1)); %calculating new values using RK method
b(i+1) = b(i) + (h/6)*(k1(2)+ 2*k2(2) + 2*k3(2) + k4(2));
end
plot(t,a)
a(3001) %the value of the function y at 30 (needs to be 20)

回答(1 个)

Torsten
Torsten 2015-11-20
编辑:Torsten 2015-11-20
b30_target = 20;
b0_high = 0;
b0_low = -3;
b30_high = 54.32;
b30_low = -15.3;
solved = 0;
while solved == 0
b0_test = b0_high + (b30_target-b30_high)/(b30_low-b30_high)*(b0_low-b0_high);
% solve BVP with initial conditions a(1)=25 and b(1)=b0_test
b30_test = b(end);
if abs(b30_test-b30_target) < 0.1
solved = 1;
else
if b30_test < b30_target
b30_low = b30_test;
b0_low = b0_test;
else
b30_high = b30_test;
b0_high = b0_test;
end
end
end
Best wishes
Torsten.
  2 个评论
Kevin Tovson
Kevin Tovson 2015-11-20
Hi Torsten. Thanks so much! Where would this code be put into my current code?
Torsten
Torsten 2015-11-23
编辑:Torsten 2015-11-23
I indicated the position where the BVP has to be solved ...
Best wishes
Torsten.

请先登录,再进行评论。

Community Treasure Hunt

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

Start Hunting!

Translated by