Solve a complex differential equation over different intervals

1 次查看（过去 30 天）
Melvin Balawi
Melvin Balawi 2021-7-23
回答： darova 2021-7-24
Hello,
I am trying to solve this third order differential equation whose constants have different values depending on the intervals considered (3). I wrote a piece of code and tried some things especially with "switch case" but I don't have the necessary knowledge to solve correctly so I gave up, if someone could help me it would be very kind.
Here is the differential equation and the code
clear all
clc
etilde=6/25;
zi=-1;
zf=1;
phii=1;
phif=-1;
n=10^4;
yinit = [phii phif];
init = bvpinit([linspace(zi,zi+etilde,n) linspace(zi+etilde,zf-etilde,n) linspace(zf-etilde,1,n)],yinit);
sol = bvp4c(@(x,y)f(x,y), @bc, init);
x=[linspace(zi,zi+etilde,n) linspace(zi+etilde,zf-etilde,n) linspace(zf-etilde,1,n)];
y = deval(sol, x);
plot(x,y(1,:));
legend('phi(z)')
title('Potentiel électrique en fonction de l épaisseur adimensionnée')
xlabel({'z'})
ylabel('phi(z)')
function dphidz = f(x,y) % equations being solved
epsilon=10^-5;
A1=[51.8 19.8 51.8];
A2=[0 3.35 0];
A3=epsilon*[0.499 1.58 0.499];
A4=[1.18E-1 1.09E-2 1.18E-1];
A5=[2.44 0.617 2.44];
A6=0.6;
Phi2m=[1 0.83 1];
etilde=6/25;
y = zeros(3,1);
dphidz = zeros(3,1);
dphidz(1)=y(2);
dphidz(2)=y(3);
if (-1<z)&&(z<-1+etilde)
% z in [-1 -1+etilde]
j=1;
dphidz(3)=A6/(A3(j)*A1(j)*Phi2m(j))*y(2)*(A1(j)-A2(j)+A3(j)*y(3))/...
(A4(j)+A5(j)/A6*(1-A1(j)*Phi2m(j)/(A1(j)-A2(j)+A3(j)*y(3))));
elseif (-1+etilde<z)&&(z<1-etilde) % z in [-1+etilde 1-etilde]
j=2;
dphidz(3)=A6/(A3(j)*A1(j)*Phi2m(j))*y(2)*(A1(j)-A2(j)+A3(j)*y(3))/...
(A4(j)+A5(j)/A6*(1-A1(j)*Phi2m(j)/(A1(j)-A2(j)+A3(j)*y(3))));
else % z in [1-etilde 1]
j=3;
dphidz(3)=A6/(A3(j)*A1(j)*Phi2m(j))*y(2)*(A1(j)-A2(j)+A3(j)*y(3))/...
(A4(j)+A5(j)/A6*(1-A1(j)*Phi2m(j)/(A1(j)-A2(j)+A3(j)*y(3))));
end
end
function res = bc(phiL,phiR)
epsilon=10^-5;
A3=epsilon*[0.499 1.58 0.499];
res = [phiL(1,1) - 1
phiR(1,1) - phiL(1,2)
A3(1)*phiR(2,1) - A3(2)*phiL(2,2)
A3(1)*phiL(2,1) - A3(2)*phiL(2,3)
A3(1)*phiR(3,1) - A3(2)*phiL(3,2)
phiR(1,2) - phiL(1,3)
A3(2)*phiR(2,2) - A3(3)*phiL(2,3)
A3(2)*phiR(3,2) - A3(3)*phiL(3,3)
phiR(1,3) + 1];
end

请先登录，再回答此问题。

回答（1 个）

darova
darova 2021-7-24
Here are some notes:
  • there is no need of creating special mesh. You have 3d order diff equation so you need 3 yinit
yinit = [phii phif 1];
% init = bvpinit([linspace(zi,zi+etilde,n) linspace(zi+etilde,zf-etilde,n) linspace(zf-etilde,1,n)],yinit);
init = bvpinit(-1:.1:1,yinit);
  • This part of the code can be shorter and more readable
  • too many boundary contions. You need only 3
% res = [phiL(1,1) - 1
% phiR(1,1) - phiL(1,2)
% A3(1)*phiR(2,1) - A3(2)*phiL(2,2)
% A3(1)*phiL(2,1) - A3(2)*phiL(2,3)
% A3(1)*phiR(3,1) - A3(2)*phiL(3,2)
% phiR(1,2) - phiL(1,3)
% A3(2)*phiR(2,2) - A3(3)*phiL(2,3)
% A3(2)*phiR(3,2) - A3(3)*phiL(3,3)
% phiR(1,3) + 1];
res = [phiL(1) - 1
phiR(1) + 1
phiR(2) - phiL(2)];

类别

在 Help Center 和 File Exchange 中查找有关 Boundary Value Problems 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by