solving BVP using anonymous function

조회 수: 3 (최근 30일)
Louisa Mezache
Louisa Mezache 2021년 3월 26일
답변: Sulaymon Eshkabilov 2021년 3월 26일
I'm trying to solve a boundary value problem using a linear shooting method, so I've linearized my nonlinear BVP. What I don't understand is how my RK4 function is passing my anonymous function (that has two outputs). Below is my code (including my RK4 function). The first element in my anonymous function is my initial slope value, so I'm struggling to understand how it's moving through my RK4. I'm also struggling to ask my question in a way that can be understood, so please let me know if you need me to try to reword it.
%Given BVP: (d^2 T)/(dx^2 ) = a(T+T_a )^4-b(T-T_b )
%linearized form: (d^2 T)/(dx^2 ) = a*((Ta+Tb)^4 + 4*(Ta+Tb)^3*(X1(1)-Tb)) + b*(X1(1)-Tb)
%define IVP: dT/dx = z ; dz/dx = a*((Ta+Tb)^4 + 4*(Ta+Tb)^3*(X1(1)-Tb)) + b*(X1(1)-Tb)
dx = 0.01; %spatial step
z01 = 5; %first guess
z02 = 100; %second guess
X1 = [T0;z01]; %vector for initial conditions (initial value and initial slope)
%taken from lecture, but how does it work?
diff_eq1 = @(x,X1) [X1(2), a*((Ta+Tb)^4 + 4*(Ta+Tb)^3*(X1(1)-Tb)) + b*(X1(1)-Tb)]; %system of first order initial value problems
[~,Y1] = rk4(diff_eq1,dx,X1,L,x0);
function [time,y] = rk4(diff_eq,h,y0,T,t0)
y(:,1) = y0;
t(1) = t0; %initial value
tf = T; %final value in interval
time=t(1):h:tf;
for i = 1:length(time)-1
k1(:,i) = diff_eq(t(i),y(:,i));
k2(:,i) = diff_eq(t(i)+h/2,y(:,i)+k1*h/2);
k3(:,i) = diff_eq(t(i)+h/2,y(:,i)+k2*h/2);
k4(:,i) = diff_eq(t(i)+h,y(:,i)+k3*h);
y(:,i+1) = y(:,i) + (1/6)*(k1(:,i) + 2*k2(:,i) + 2*k3(:,i) + k4(:,i))*h;
t(i+1) = t(i)+h;
y; %since y is being indexed, y will output values of y
end
end

답변 (1개)

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2021년 3월 26일
Your four variables a, b, Ta and Tb are not defined. Pre-define their values just before your anonymous function is defined.
Good luck.

카테고리

Help CenterFile Exchange에서 Boundary Value Problems에 대해 자세히 알아보기

제품


릴리스

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by