Matlab to solve finite difference for Ode
3 views (last 30 days)
Show older comments
% Finite difference method
% Approximate the solution of y"=(-2/x)y'+(2/x^2)y+ sin(lnx)/x^2
% for 1<=x<=2 with y(1)=1 and y(2)=2.
p = @(x) (-2/x);
q = @(x) (2/x^2);
r = @(x) (sin(log(x))/x^2);
aa = 1; bb = 2; alpha = 1; beta = 2; n=29; % h = (bb-aa)/(n+1); h=0.1
fprintf(' x w \n');
h = (bb-aa)/(n+1);
a = zeros(1,n+1);
b = zeros(1,n+1);
c = zeros(1,n+1);
d = zeros(1,n+1);
l = zeros(1,n+1);
u = zeros(1,n+1);
z = zeros(1,n+1);
w = zeros(1,n+1);
x = aa+h;
a(1) = 2+h^2*q(x);
b(1) = -1+0.5*h*p(x);
d(1) = -h^2*r(x)+(1+0.5*h*p(x))*alpha;
m = n-1;
for i = 2 : m
x = aa+i*h;
a(i) = 2+h^2*q(x);
b(i) = -1+0.5*h*p(x);
c(i) = -1-0.5*h*p(x);
d(i) = -h^2*r(x);
end
x = bb-h;
a(n) = 2+h^2*q(x);
c(n) = -1-0.5*h*p(x);
d(n) = -h^2*r(x)+(1-0.5*h*p(x))*beta;
l(1) = a(1);
u(1) = b(1)/a(1);
z(1) = d(1)/l(1);
for i = 2 : m
l(i) = a(i)-c(i)*u(i-1);
u(i) = b(i)/l(i);
z(i) = (d(i)-c(i)*z(i-1))/l(i);
end
l(n) = a(n)-c(n)*u(n-1);
z(n) = (d(n)-c(n)*z(n-1))/l(n);
w(n) = z(n);
for j = 1 : m
i = n-j;
w(i) = z(i)-u(i)*w(i+1);
end
i = 0;
fprintf('%5.4f %11.8f\n', aa, alpha);
for i = 1 : n
x = aa+i*h;
fprintf('%5.4f %11.8f\n', x, w(i));
end
i = n+1;
fprintf('%5.4f %11.8f\n', bb, beta);am trying to run d code but it keeps writing Error: Invalid text character. The text ' ' contains an unsupported non-ASCII whitespace character. Please where am I doing wrong?
0 Comments
Answers (2)
Viranch Patel
on 9 Jul 2021
I tried running the code on my machine and it's working fine. You can paste this code to your machine and see if its working. I am attaching the output too. You might be using some wrong character for any particular character which is not supported.
% Finite difference method
% Approximate the solution of y"=(-2/x)y'+(2/x^2)y+ sin(lnx)/x^2
% for 1<=x<=2 with y(1)=1 and y(2)=2.
p = @(x) (-2/x);
q = @(x) (2/x^2);
r = @(x) (sin(log(x))/x^2);
aa = 1; bb = 2; alpha = 1; beta = 2; n=29; % h = (bb-aa)/(n+1); h=0.1
fprintf(' x w \n');
h = (bb-aa)/(n+1);
a = zeros(1,n+1);
b = zeros(1,n+1);
c = zeros(1,n+1);
d = zeros(1,n+1);
l = zeros(1,n+1);
u = zeros(1,n+1);
z = zeros(1,n+1);
w = zeros(1,n+1);
x = aa+h;
a(1) = 2+h^2*q(x);
b(1) = -1+0.5*h*p(x);
d(1) = -h^2*r(x)+(1+0.5*h*p(x))*alpha;
m = n-1;
for i = 2 : m
x = aa+i*h;
a(i) = 2+h^2*q(x);
b(i) = -1+0.5*h*p(x);
c(i) = -1-0.5*h*p(x);
d(i) = -h^2*r(x);
end
x = bb-h;
a(n) = 2+h^2*q(x);
c(n) = -1-0.5*h*p(x);
d(n) = -h^2*r(x)+(1-0.5*h*p(x))*beta;
l(1) = a(1);
u(1) = b(1)/a(1);
z(1) = d(1)/l(1);
for i = 2 : m
l(i) = a(i)-c(i)*u(i-1);
u(i) = b(i)/l(i);
z(i) = (d(i)-c(i)*z(i-1))/l(i);
end
l(n) = a(n)-c(n)*u(n-1);
z(n) = (d(n)-c(n)*z(n-1))/l(n);
w(n) = z(n);
for j = 1 : m
i = n-j;
w(i) = z(i)-u(i)*w(i+1);
end
i = 0;
fprintf('%5.4f %11.8f\n', aa, alpha);
for i = 1 : n
x = aa+i*h;
fprintf('%5.4f %11.8f\n', x, w(i));
end
i = n+1;
fprintf('%5.4f %11.8f\n', bb, beta);
Output:
x w
1.0000 1.00000000
1.0333 1.03067949
1.0667 1.06155254
1.1000 1.09262608
1.1333 1.12390423
1.1667 1.15538887
1.2000 1.18708018
1.2333 1.21897697
1.2667 1.25107701
1.3000 1.28337728
1.3333 1.31587416
1.3667 1.34856355
1.4000 1.38144105
1.4333 1.41450201
1.4667 1.44774163
1.5000 1.48115505
1.5333 1.51473732
1.5667 1.54848354
1.6000 1.58238883
1.6333 1.61644834
1.6667 1.65065735
1.7000 1.68501118
1.7333 1.71950528
1.7667 1.75413522
1.8000 1.78889666
1.8333 1.82378540
1.8667 1.85879736
1.9000 1.89392857
1.9333 1.92917520
1.9667 1.96453354
2.0000 2.00000000
0 Comments
Walter Roberson
on 1 Dec 2024
S = 'aa = 1; bb = 2; alpha = 1; beta = 2; n=29; % h = (bb-aa)/(n+1); h=0.1'
idx = find(S>128)
S(idx)
S(idx)+0
Your code contains a number of char(160), which is the Unicode Non-Breaking Whitespace character https://www.compart.com/en/unicode/U+00A0
The first group of NBSP is after the n=29;
S(40:50)
0 Comments
See Also
Categories
Find more on Function Creation in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!