Energy method with Hermite cubic function (2 essential condition)

1 view (last 30 days)
Aurora
Aurora on 5 Nov 2023
Commented: Aurora on 5 Nov 2023
The different proble is:
-(x^2 y')'+20y=7*x^2 -5
y(-1)=0
y(1)=0
sol: y=(x^2)/2-(x^4)/4-1/4
The next code has the following problem:
Unrecognized function or variable 'm'.
Error in ex_enerHermite>support (line 170)
if l>=2*m
Error in ex_enerHermite (line 10)
[lbl,rbl]=support(l,m);
how can I solve it? Where is the errore?
function [x,u]=ex_enerHermite(m,nv)
%problema con due condizioni essenziali
h=2/m;
xc=-1:h:1;
A=zeros(2*m);%se non ho condizioni essenziali diventa 2*m+2
b=zeros(2*m,1);%se non ho condizioni essenziali diventa 2*m+2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Caso condizione essenziale in a
for l=1:2*m %se non ho condizioni essenziali diventa 2*m+2
[lbl,rbl]=support(l,m);
for j=1:2*m %se non ho condizioni essenziali diventa 2*m+2
if abs(l-j)<4 %altrimenti l'integrale fa zero causa supporto base
[lbj,rbj]=support(j,xc);
lb=max(lbl,lbj);
rb=min(rbl,rbj);
tol=10^-8;
if abs(lb-rb)>=tol
% termine di bordo già inserito alla linea 26
A(l,j)=integral(@(x)fA(x,l,j,xc),lb,rb);
end
end
end
b(l)=integral(@(x)fb(x,l,xc),lbl,rbl);
end
delta=A\b;
h=2/nv;
x=-1:h:1;
u=zeros(1,size(x));
for l=1:2*m %da cambiare in base a dove fai la sommatoria della soluzione u
u=u+delta(l)*fu(x,l,xc);
%se cond. ess. in a allora l=1:2*m+1 e
%u=u+delta(l)*fu(x,l,xc)
%se non ho cond. ess. allora l=0:2*m+1 e
%u=u+delta(l+1)*fu(x,l,xc)
end
%u è la soluzione del problema traslato quindi se necessario aggiungo l(x)
end
function y=fH(x,l,xc) %non toccare
m=length(xc)-1;
y=zeros(size(x));
for j=1:length(x)
if l==0
if xc(1)<=x(j) && x(j)<xc(2)
y(j)=3*(((xc(2)-x(j))/(xc(2)-xc(1)))^2)-...
2*(((xc(2)-x(j))/(xc(2)-xc(1)))^3);
end
elseif l==m
if xc(m)<x(j) && x(j)<=xc(m+1)
y(j)=3*(((x(j)-xc(m))/(xc(m+1)-xc(m)))^2)-...
2*(((x(j)-xc(m))/(xc(m+1)-xc(m)))^3);
end
else
if xc(l)<x(j) && x(j)<=xc(l+1)
y(j)=3*(((x(j)-xc(l))/(xc(l+1)-xc(l)))^2)-...
2*(((x(j)-xc(l))/(xc(l+1)-xc(l)))^3);
elseif xc(l+1)<=x(j) && x(j)<xc(l+2)
y(j)=3*(((xc(l+2)-x(j))/(xc(l+2)-xc(l+1)))^2)-...
2*(((xc(l+2)-x(j))/(xc(l+2)-xc(l+1)))^3);
end
end
end
end
function y=fS(x,l,xc) %non toccare
m=length(xc)-1;
y=zeros(size(x));
for j=1:length(x)
if l==0
if xc(1)<=x(j) && x(j)<xc(2)
y(j)=((xc(2)-x(j))^2/(xc(2)-xc(1)))-...
((xc(2)-x(j))^3/(xc(2)-xc(1))^2);
end
elseif l==m
if xc(m)<x(j) && x(j)<=xc(m+1)
y(j)=-((x(j)-xc(m))^2/(xc(m+1)-xc(m)))+...
((x(j)-xc(m))^3/(xc(m+1)-xc(m))^2);
end
else
if xc(l)<x(j) && x(j)<=xc(l+1)
y(j)=-((x(j)-xc(l))^2/(xc(l+1)-xc(l)))+...
((x(j)-xc(l))^3/(xc(l+1)-xc(l))^2);
elseif xc(l+1)<=x(j) && x(j)<xc(l+2)
y(j)=((xc(l+2)-x(j))^2/(xc(l+2)-xc(l+1)))-...
((xc(l+2)-x(j))^3/(xc(l+2)-xc(l+1))^2);
end
end
end
end
function y=dH(x,l,xc) %non toccare
m=length(xc)-1;
y=zeros(size(x));
for j=1:length(x)
if l==0
if xc(1)<=x(j) && x(j)<xc(2)
y(j)=(6/(xc(2)-xc(1)))*((xc(2)-x(j))/(xc(2)-xc(1)))*...
((xc(2)-x(j))/(xc(2)-xc(1))-1);
end
elseif l==m
if xc(m)<x(j) && x(j)<=xc(m+1)
y(j)=(6/(xc(m+1)-xc(m)))*((x(j)-xc(m))/(xc(m+1)-xc(m)))*...
(1-((x(j)-xc(m))/(xc(m+1)-xc(m))));
end
else
if xc(l)<x(j) && x(j)<=xc(l+1)
y(j)=(6/(xc(l+1)-xc(l)))*((x(j)-xc(l))/(xc(l+1)-xc(l)))*...
(1-((x(j)-xc(l))/(xc(l+1)-xc(l))));
elseif xc(l+1)<=x(j) && x(j)<xc(l+2)
y(j)=(6/(xc(l+2)-xc(l+1)))*((xc(l+2)-x(j))/(xc(l+2)-xc(l+1)))*...
((xc(l+2)-x(j))/(xc(l+2)-xc(l+1))-1);
end
end
end
end
function y=dS(x,l,xc) %non toccare
m=length(xc)-1;
y=zeros(size(x));
for j=1:length(x)
if l==0
if xc(1)<=x(j) && x(j)<xc(2)
y(j)=3*(((xc(2)-x(j))/(xc(2)-xc(1)))^2)-...
2*((xc(2)-x(j))/(xc(2)-xc(1)));
end
elseif l==m
if xc(m)<x(j) && x(j)<=xc(m+1)
y(j)=3*(((x(j)-xc(m))/(xc(m+1)-xc(m)))^2)-...
2*((x(j)-xc(m))/(xc(m+1)-xc(m)));
end
else
if xc(l)<x(j) && x(j)<=xc(l+1)
y(j)=3*(((x(j)-xc(l))/(xc(l+1)-xc(l)))^2)-...
2*((x(j)-xc(l))/(xc(l+1)-xc(l)));
elseif xc(l+1)<=x(j) && x(j)<xc(l+2)
y(j)=3*(((xc(l+2)-x(j))/(xc(l+2)-xc(l+1)))^2)-...
2*((xc(l+2)-x(j))/(xc(l+2)-xc(l+1)));
end
end
end
end
function y=fu(x,l,xc) %scrivi quello che hai trovato nei calcoli della base u_l, attento all'indice l
if l==2*m
y=fS(x,m,xc);
elseif mod(l,2)==1
y=fS(x,(l-1)/2,xc);
else
y=fH(x,l/2,xc);
end
end
function y=du(x,l,xc) %derivata di u_l
if l==2*m
y=dS(x,m,xc);
elseif mod(l,2)==1
y=dS(x,(l-1)/2,xc);
else
y=dH(x,l/2,xc);
end
end
function [lb,rb]=support(l,xc)
%se abbiamo u_0 aggiungo if l==0 lb=xc(1); rb=xc(2);
if l>=2*m
lb=xc(m);
rb=xc(m+1);
elseif mod(l,2)==1 %sono le funzioni S
%In generale alcuni supporti sono più piccoli ma questo va bene per
%tutti, idem nelle funzioni H
lb=xc((l-1)/2);
rb=xc(2+(l-1)/2);
else %sono le funzioni H
lb=xc(l/2);
rb=xc(2+l/2);
end
end
function y=fA(x,l,j,xc)
%Inserisci la funzione integranda per la matrice A
y=(x.^2).*du(x,l,xc).*du(x,j,xc)+20.*fu(x,l,xc).*fu(x,j,xc);
end
function y=fb(x,l,xc)
%Inserisci la funzione integranda per il vettore b
y=(7.*(x.^2)-5).*fu(x,l,xc);
end

Sign in to answer this question.

Accepted Answer

Torsten
Torsten on 5 Nov 2023
m and nv are inputs to the function "ex_enerHermite".
Where do you set these values and where do you call the function ?
  5 Comments
Show 3 older comments
Aurora
Aurora on 5 Nov 2023
You can choose m, nv like you prefer. In particoular, i choose 15,6
Aurora
Aurora on 5 Nov 2023
>> [x,u]=ex_enerHermite2(30,6)
x =
-2.0000 -1.7500 -1.5000 -1.2500 -1.0000 -0.7500 -0.5000
u =
-2.9103 -2.6970 -2.4510 -2.1745 -1.8722 -1.5543 -1.2500
>> y=((x.^2)/3)+(2.*x)+1./(6.*x)
y =
-2.7500 -2.5744 -2.3611 -2.1125 -1.8333 -1.5347 -1.2500
>> plot(x,u,'r*',x,y)

Sign in to comment.

More Answers (0)

Categories

Find more on Robotics in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!