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;

4 ビュー (過去 30 日間)
Nauryzbay
Nauryzbay 2023 年 6 月 22 日
編集済み: Torsten 2023 年 6 月 23 日
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)]);%

回答 (1 件)

Torsten
Torsten 2023 年 6 月 22 日
編集済み: Torsten 2023 年 6 月 22 日
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 件のコメント
Torsten
Torsten 2023 年 6 月 22 日
編集済み: Torsten 2023 年 6 月 23 日
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
Nauryzbay
Nauryzbay 2023 年 6 月 23 日
Thank you very much. You really helped me.

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeCalculus についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by