Finite difference method 1st degree
6 views (last 30 days)
Show older comments
I have to solve this problem using a first degree finite difference method. I am really struggling to understand what is not working with my code. If someone could take a look at ir and explain to me what<s wrong I would really appreciate it. Here is the assignement and the code below
function[y]=problimite(N,P,Q,R, a, b, alpha, beta)
h=(b-a)/(N+1);
S=-1+P(1:end-1)*h/2;
D=2+Q(1:end)*h^2;
I=-1-P(2:end)*h/2
B=[(-R(1)*(h^2)+(1+P(1)*h/2))*alpha,-R(2:end-1).*h^2,-R(end)*(h^2)+(1+P(end)*h/2)*beta];
disp(S)
y = tridiagonal(N, S, D, I, B);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
function y = tridiagonal(N, S, D, I, B)
y(1)=B(1)/D(1);
for i=2:N
S(i-1)=(S(i-1))/D(i-1);
D(i)=D(i)-I(i-1)*S(i-1);
y(i)=(B(i)-I(i-1)*y(i-1))/D(i);
end
x(N)=y(N);
for i=N-1:1
x(i)=y(i)-S(i)*x(i+1);
end
%%script%%%%
b=1;
alpha=0;
beta=1;
X=[0.9:0.01:1];
p=@(x) -1/x;
q=@(x) 0;
r=@(x) (1.6)/(x^4);
N=9;
P=arrayfun(p,X);
R=arrayfun(r,X);
Q=arrayfun(q,X);
%h=.01;
%I=[-1-P(2:end)*h/2]
%D=[2+Q(1:end)*h^2]
%S=[-1+P(1:end-1)*h/2]
%display(P);
%display(R);
%display(Q);
plot(problimite(N,P,Q,R, a, b, alpha, beta))
hold on
%% partie 2 %%
X=[0.9:0.005:1];
p=@(x) -1/x;
q=@(x) 0;
r=@(x) (1.6)/(x^4);
N=19;
P=arrayfun(p,X);
R=arrayfun(r,X);
Q=arrayfun(q,X);
%display(P);
%display(R);
%display(Q);
plot(problimite(N,P,Q,R, a, b, alpha, beta))
2 Comments
J. Alex Lee
on 7 Dec 2019
@darova, this looks like a homework assignment requiring hand-coding; it even appears to have an analytical solution for verification.
@Patirick, what do you mean by "first degree" finite difference method? Do you mean first order accuracy? If so I'm not sure there is a first order accurate difference formula for second order derivatives...You can break up your problem into a pair of coupled first order ODEs and use a pair of first order accurate finite differencing schemes. Perhaps you can explain in words the strategy that your code is implementing before jumping into the coding.
Answers (0)
See Also
Categories
Find more on Animation 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!