Energy method with Hermite cubic function (2 essential condition)
Show older comments
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
Accepted Answer
Torsten
on 5 Nov 2023
0 votes
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
Aurora
on 5 Nov 2023
Aurora
on 5 Nov 2023
Hello, the problem is change. The approximatly solution is a littlet different.
In the followinf, you can see the listing and the output.
function [x,u]=ex_enerHermite2(m,nv)
%trasla il problema se necessario con condizione essenziale in b
B=11/14;
C=-6/7;
h=1.5/m;
xc=-2:h:-0.5;
A=zeros(2*m+1);%se non ho condizioni essenziali diventa 2*m+2
b=zeros(2*m+1,1);%se non ho condizioni essenziali diventa 2*m+2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Caso condizione essenziale in a
for l=1:2*m+1 %se non ho condizioni essenziali diventa 2*m+2
[lbl,rbl]=support(l,xc);
for j=1:2*m+1 %se non ho condizioni essenziali diventa 2*m+2
if abs(l-j)<4 %altrimenti l'integrale è zero causa supporto base
[lbj,rbj]=support(j,xc);
lb=max(lbl,lbj);
rb=min(rbl,rbj);
A(l,j)=integral(@(x)fA(x,l,j,xc),lb,rb);
end
end
b(l)=integral(@(x)fb(x,l,xc),lbl,rbl);
end
A(1,1)=2+A(1,1);
%H_m contiene anche il termine fuori dall'integrale
%da aggiungere.
delta=A\b;
h=1.5/nv;
x=-2:h:-0.5;
u=zeros(1,nv+1);
for l=1:2*m+1 %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)
u=u+B*x+C;
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
m = length(xc)-1;
if l==2*m+1
y=fS(x,m,xc);
elseif mod(l,2)==1
y=fS(x,(l-1)/2,xc);
else
y=fH(x,(l-2)/2,xc);
end
end
function y=du(x,l,xc) %derivata di u_l
m = length(xc)-1;
if l==2*m+1
y=dS(x,m,xc);
elseif mod(l,2)==1
y=dS(x,(l-1)/2,xc);
else
y=dH(x,(l-2)/2,xc);
end
end
function [lb,rb]=support(l,xc)
m=length(xc)-1;
%se abbiamo u_0 aggiungo if l==0 lb=xc(1); rb=xc(2);
if l==1 || l==2
lb=xc(1);
rb=xc(2);
elseif l==2*m+1 %cioè il fondo, se non ho u_2m+1 posso scrivere 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)/2);
rb=xc(2+(l-2)/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)+6.*fu(x,l,xc).*fu(x,j,xc);
end
function y=fb(x,l,xc)
B=11/14;
C=-6/7;
%Inserisci la funzione integranda per il vettore b
y=(8.*x+1./x-4*B.*x-6*C).*fu(x,l,xc);
end
Torsten
on 5 Nov 2023
We still cannot run your code and reproduce the plot since we don't know m and nv in your call to "ex_enerHermite2".
Aurora
on 5 Nov 2023
Aurora
on 5 Nov 2023
More Answers (0)
Categories
Find more on Programming in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
