Using Shooting Method on Euler Method to Solve Second Order BVP

23 views (last 30 days)
% Project 3 Shooting Method Using Euler Method
%y''+(1-x)y'+xy=x, y(0) =0, y(1) =2
clear; clc
%parameters
x0 = 0;
y0 = 0;
xf = 1;
u0 = 1;%input('yp value at intial x value:\n');
y1 = 2; %input('Input second y condition value:\n');
n = 10;
x(1) = x0;
y(1) = y0;
u(1) = u0;
h = (xf-x0)/n;
t = x0:h:xf ;
for i = 2:length(t)
x(i) = x(i-1)+h;
y(i) = y(i-1)+u(i-1)*h;
u(i) = u(i-1)+h*(x(i-1)-(x(i-1)*y(i-1))-(1-x(i-1)*u(i-1)));
end
A1 = [x;y];
fprintf('\n%8.4f %8.4f\n', A1)
plot(x,y,'r')
This code numerically solves the second ODE using shooting method on Euler Method. However, I'm trying to figure out how to get matlab to automatically test y' values at x0 (I have this variable as u0 in my code) in the iterations until it finds a y' value that satisfies the second boundry condition at y(1). I tried using linspace for u0 to test a string of numbers but it isn't working.

Accepted Answer

Bjorn Gustavsson
Bjorn Gustavsson on 11 May 2022
Edited: Bjorn Gustavsson on 11 May 2022
To do that you should make functions out of this script. That way you will have a much cleaner programme, that is far easier to generalize and extend.
For this example you could wrap your Euler-integrator into one such function:
function [x,y] = your_Euler_odeint(x0,xf,y0,dydx0,n)
x(1) = x0;
y(1) = y0;
u(1) = dydx0;
h = (xf-x0)/n;
t = x0:h:xf ;
for i = 2:length(t)
x(i) = x(i-1)+h;
y(i) = y(i-1)+u(i-1)*h;
u(i) = u(i-1)+h*(x(i-1)-(x(i-1)*y(i-1))-(1-x(i-1)*u(i-1)));
end
end
That function should be possible to call from your script:
x0 = 0;
y0 = 0;
xf = 1;
dydx0 = 1;
y1 = 2;
n = 10;
[x,y] = your_Euler_odeint(x0,xf,y0,dydx0,n);
y_diff_from_target_at_end = y(end)-y1;
Once you've done this you can see the shooting-method as an equation-solving task. Matlab has a large suite of functions to use for that, but from the scope of this task it seems as you should do that yourself from very basic programming. So you could simply do some test and adjustment, Newton-solving style, just test some dydx0 and see what value that gives, test a slightly larger or smaller dydx0 and see if that takes you closer to the solution, and then make the Newton-guess for where the solution might be. Here's a brief description: Newton's method
HTH

More Answers (0)

Community Treasure Hunt

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

Start Hunting!