How to solve a boundary value problem x"+x'=1, x(0)=1, x(1)=2 , with an impulse condition x(i/3+)=2*x(i/3)+1, x'(i/3+)=3*x'(i/3)-1;

3 views (last 30 days)
function impulsbound %x"+x'=1, x(0)=1, x(i/3+)=2*x(i/3)+1, x'(i/3+)=3*x'(i/3)-1, i=1,2. x(1)=2;
format long
clear all
t=linspace(0,1,100);
solinit = bvpinit(t,[1 2]) ;
x1=1;
u=2;
teta(1)=0;
%teta(u+2)=0.2;
for i=1:u
teta(i+1)=i/3;
end
for j=1:1
for k=1:u
%initial_func=[x1,x2];
%[t,x] = ode45(@s,[teta(k):0.0001:teta(k+1)],[x1(j) x2(j)]);
x = bvp4c(@s,@bc2,solinit,[teta(k):0.0001:teta(k+1)],[x1(j) x2(j)]);%
n=length(t);
%disp(size(x));
x1(j)=x(n,1)+x(n,1)+1;
x2(j)=x(n,2)+2*x(n,2)-1;
hold on
%view(30,15);
hold on
figure(1)
subplot(2,1,1);
plot(t,x(:,1),'color','r','Linewidth',1.2);
xlabel('\bf t'); ylabel('$$z$$','interpreter','latex','fontsize',16); zlabel('\bf \psi_2');
grid on
hold on
subplot(2,1,2);
plot(t,x(:,2),'color','r','Linewidth',1.2);
xlabel('\bf t'); ylabel('$$y$$','interpreter','latex','fontsize',16); zlabel('\bf \psi_3');
grid on
hold on
end
end
function dx=s(t,x)
dx=zeros(2,1); % создает нулевой вектор-столбец
dx(1)=x(2);
dx(2)=1-x(1);
function res = bc2(xa, xb)
res = [xa(1)-1
xb(1)-2];
Reference to a cleared variable x2.
Error in impulsbound (line 17)
x = bvp4c(@s,@bc2,solinit,[teta(k):0.0001:teta(k+1)],[x1(j) x2(j)]);%

Answers (1)

Torsten
Torsten on 22 Jun 2023
Edited: Torsten on 22 Jun 2023
syms t x1(t) x2(t)
eqn1 = diff(x1,t,2) + diff(x1,t) == 1;
Dx1 = diff(x1,t);
eqn2 = diff(x2,t,2) + diff(x2,t) == 1;
Dx2 = diff(x2,t);
eqns = [eqn1,eqn2];
conds = [x1(0) == 1, x2(1/3) == 2*x1(1/3)+1, Dx2(1/3) == 3*Dx1(1/3)-1, x2(1) == 2];
sol = dsolve(eqns,conds)
sol = struct with fields:
x2: t + (2*exp(1) + 2*exp(1/3) + exp(2/3) - 12)/(2*exp(1) + exp(2/3) - 3) - (exp(-t)*exp(1/3)*(2*exp(1) - 9*exp(2/3)))/(2*exp(1) + exp(2/3) - 3) x1: t - (4*exp(1) - 3*exp(1/3) - 3*exp(2/3) + 9)/(3*(2*exp(1) + exp(2/3) - 3)) + (exp(-t)*(10*exp(1) - 3*exp(1/3)))/(3*(2*exp(1) + exp(2/3) - 3))
hold on
fplot(sol.x1,[0,1/3])
fplot(sol.x2,[1/3,1])
hold off
grid on
  3 Comments
Torsten
Torsten on 22 Jun 2023
Edited: Torsten on 23 Jun 2023
Then you will have to write a little more: you will have equations and transmission conditions for x1(t),...,x10(t).
And also the above code is set up for only one transmission condition at x=1/3, not both at x=1/3 and x=2/3. So already in this case, you have to use x1(t), x2(t) and x3(t) as follows:
syms t x1(t) x2(t) x3(t)
eqn1 = diff(x1,t,2) + diff(x1,t) == 1;
Dx1 = diff(x1,t);
eqn2 = diff(x2,t,2) + diff(x2,t) == 1;
Dx2 = diff(x2,t);
eqn3 = diff(x3,t,2) + diff(x3,t) == 1;
Dx3 = diff(x3,t);
eqns = [eqn1,eqn2,eqn3];
conds = [x1(0) == 1, ...
x2(1/3) == 2*x1(1/3)+1, Dx2(1/3) == 3*Dx1(1/3)-1,...
x3(2/3) == 2*x2(2/3)+1, Dx3(2/3) == 3*Dx2(2/3)-1,...
x3(1) == 2];
sol = dsolve(eqns,conds)
sol = struct with fields:
x2: t - (10*exp(1) - 57*exp(1/3) + 5*exp(2/3) + 108)/(3*(4*exp(1) + 3*exp(1/3) + 2*exp(2/3) - 9)) - (exp(-t)*exp(1/3)*(4*exp(1) + 3*exp(1/3) - 29*exp(2/3)))/(4*exp(1) + 3*exp(1/3) + 2*exp(2/3)… x1: t - (19*exp(1) - 18*exp(1/3) - 6*exp(2/3) + 27)/(3*(4*exp(1) + 3*exp(1/3) + 2*exp(2/3) - 9)) + (exp(-t)*(31*exp(1) - 9*exp(1/3)))/(3*(4*exp(1) + 3*exp(1/3) + 2*exp(2/3) - 9)) x3: t + (4*exp(1) + 17*exp(1/3) + 6*exp(2/3) - 93)/(4*exp(1) + 3*exp(1/3) + 2*exp(2/3) - 9) - (2*exp(-t)*exp(2/3)*(2*exp(1) - 42*exp(1/3) + 7*exp(2/3)))/(4*exp(1) + 3*exp(1/3) + 2*exp(2/3) - 9)
hold on
fplot(sol.x1,[0,1/3])
fplot(sol.x2,[1/3,2/3])
fplot(sol.x3,[2/3,1])
hold off
grid on

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!